0

Having a volume implicitly defined by

x*y*z <= 1 

for

-5 <= x <= 5 
-5 <= y <= 5 
-5 <= z <= 5 

how would I go about plotting its outer surface using available Python modules, preferably mayavi?

I am aware of the function mlab.mesh, but I don't understand its input. It requires three 2D arrays, that I don't understand how to create having the above information.

EDIT:

Maybe my problem lies with an unsufficient understanding of the meshgrid()-function or the mgrid-class of numpy. I see that I have to use them in some way, but I do not completely grasp their purpose or what such a grid represents.

EDIT:

I arrived at this:

import numpy as np

from mayavi import mlab
x, y, z = np.ogrid[-5:5:200j, -5:5:200j, -5:5:200j]
s = x*y*z

src = mlab.pipeline.scalar_field(s)

mlab.pipeline.iso_surface(src, contours=[1., ],)
mlab.show()

This results in an isosurface (for x*y*z=1) of a volume though, which is not quite what I was looking for. What I am looking for is basically a method to draw an arbitrary surface, like a "polygon in 3d" if there is such a thing.

I created the following code, which plots a surface (works with mayavi, too). I would need to modify this code to my particular problem, but to do that I need to understand why and how a 3d surface is defined by three 2d-arrays? What do these arrays (x, y and z) represent?

import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d, Axes3D

phi, theta = np.mgrid[0:np.pi:11j, 0:2*np.pi:11j]
x = np.sin(phi) * np.cos(theta)
y = np.sin(phi) * np.sin(theta)
z = np.cos(phi)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(x,y,z)
fig.show()
karlson
  • 5,325
  • 3
  • 30
  • 62
  • 2
    So, what have you tried? You will get much better answers here if you show you have done some research and at least _tried_ before posting. Your question currently reads as 'please do my work for me' (even if that is not your intent) which tends to annoy people which makes them less likely to answer your question. – tacaswell Jul 31 '13 at 12:43
  • I've already read through the examples on the mayavi website, e.g. the demo on [link](http://docs.enthought.com/mayavi/mayavi/mlab.html), but there the surface is defined explicitly. – karlson Jul 31 '13 at 12:50
  • Have you tried solving your implicit definition explicitly, i.e. defining domains within x,y,z, and then plotting an explicit surface for each domain? – Ludo Jul 31 '13 at 13:41
  • ok, so show us the code you have tried (even if it hasn't worked) – tacaswell Jul 31 '13 at 13:46
  • Well, the code I've had would not be of much help here, but I will get back to this (hopefully today) and will present my best attempt. Would anyone happen to know a (or some) good resources to learn about the different use cases of `np.meshgrid` and `np.mgrid`, preferably with examples like this? – karlson Aug 02 '13 at 08:42

2 Answers2

2

The outer surface, implicitly defined by

x*y*z = 1,

cannot be defined explicitly globally. To see this, consider x and y given, then:

z = 1/(x*y),

which is not defined for x = 0 or y = 0. Therefore, you can only define your surface locally for domains that do not include the singularity, e.g. for the domain

0 < x <= 5
0 < y <= 5

z is indeed defined (a hyperbolic surface). Similarly, you need to plot the surfaces for the other domains, until you have patched together

-5 <= x <= 5
-5 <= y <= 5

Note that your surface is not defined for x = 0 and y = 0, i.e. the axis of your coordinate system, so you cannot patch your surfaces together to get a globally defined surface.

Using numpy and matplotlib, you can plot one of these surfaces as follows (adopted from http://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html#surface-plots):

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()
ax = fig.gca(projection='3d')

X = np.arange(0.25, 5, 0.25)
Y = np.arange(0.25, 5, 0.25)
X, Y = np.meshgrid(X, Y)
Z = 1/(X*Y)

surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
        linewidth=0, antialiased=False)
ax.set_zlim(0, 10)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')    
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()

I'm not familiar with mayavi, but I would assume that creating the meshes with numpy would work the same.

Ludo
  • 813
  • 1
  • 9
  • 21
  • So basically, I'd need to solve the equation `xyz=1` for each variable and plot the resulting functions separately? And then maybe patch the whole thing up, because I have poles along all three axes? – karlson Aug 02 '13 at 08:36
  • Thank you @Ludo, that answer helped me a lot in understanding the problem, mathematically. The volume would then be cut off across the axes at -5 and +5. I did not accept the answer though, because it is no solution to my initial question, which was meant more like "How is a surface to be defined, to be plottable with the `mlab.mesh` function?" – karlson Aug 02 '13 at 08:58
  • You don't need to solve for all the variables, just for enough variables until you can patch the whole thing together. In this case, solve for `z` in four domains. I don't have particular knowledge of mayavi, but I guess that it's API wouldn't be much different from say matplotlib. I don't have access to either right now (at work), but can get back later tonight if you're still stuck. – Ludo Aug 02 '13 at 09:16
  • I added an example using matplotlib. I would assume that the mesh would work for mayavi also, although I don't have experience with that. – Ludo Aug 02 '13 at 21:50
1

The test case in the Mayavi docs where the function test_mesh() is defined is capable of producing a sphere. This is done by replacing

r = sin(m0*phi)**m1 + cos(m2*phi)**m3 + sin(m4*theta)**m5 + cos(m6*theta)**m7

with r = 1.0 say.

However, your problem is you need to understand that the equations you are writing define a volume when you want to draw a sphere. You need to reformulate them to give a parametric equation of a sphere. This is essentially what is done in the above example, but it may be worth your while to try it yourself. As a hint consider the equation of a circle and extend it.

sodd
  • 12,482
  • 3
  • 54
  • 62
Greg
  • 11,654
  • 3
  • 44
  • 50
  • I do not want to draw a sphere, but the surface of the volume defined by `x*y*z<=1`, i.e. the surface defined by `x*y*z=1`. My problem is exactly the thing @Ludo addressed in the comment. How would I go about defining said domains? I am not much of a math expert, but I'd like to become one eventually :), so any hints as to how to do this would be much appreciated. – karlson Aug 01 '13 at 08:55
  • 1
    My apologies I had not properly understood the question before answering. @Ludo has provided a much better answer which addresses your problems. I would advise you to first try drawing the bounds in 2D to understand the problem. – Greg Aug 02 '13 at 08:12
  • aestrivex, thanks for your very helpful comment... And thank you @Greg, that is indeed a very good idea I will try out. – karlson Aug 02 '13 at 08:40