此教程的例子讲解当使用 scipy 进行数值计算时如何利用 Mayavi 交互地可视化 numpy 数组。此教程假设你熟悉 Python 的数值工具，然后讲解如何使用 Mayavi 与这些工具协同工作。
让我们来研究电势(potential)中例子的轨迹。这在物理和工程中是很常见的问题，电势和轨迹可视化的关键是理解此问题。
我们所关心的电势是一个周期性栅格(periodic lattice)， 沉浸在抛物线约束(parabolic confinement)中。我们振动此电势然后查看粒子从一个栅格的洞(hole)跳到另一个的情况。抛物线约束在那里限制粒子的偏移 (excursion)。:
import numpy as np
def V(x, y, z):
""" A 3D sinusoidal lattice with a parabolic confinement. """
return np.cos(10*x) + np.cos(10*y) + np.cos(10*z) + 2*(x**2 + y**2 + z**2)
既然我们定义好了电势，我们看看在三维中看起来是什么样子。我们可以创建一个点组成的三维网格，然后在这些点处进行采样:
X, Y, Z = np.mgrid[-2:2:100j, -2:2:100j, -2:2:100j]
V(X, Y, Z)
我们将使用 mlab 模块(请参考simple-scripting-with-mlab)交互地可视化此容量数据(volumetric data)。为此，最好在 Python 的交互命令行输入命令，使用 Mayavi2 内建的命令行或 ipython -wthread。下面我们可视化电势的三维 iso 曲面:
from enthought.mayavi import mlab
mlab.contour3d(X, Y, Z, V)
我们可以与上面的命令产生的可视化结果进行交互，例如旋转视图。但是为了更好的理解电势的结构，变化 iso 曲面是很有用的。我们可以在 Mayavi 管道树(如果你从 ipython 中运行，需要单击场景中的 Mayavi 图标以弹出管道树) 中双击`IsoSurface`。这将打开一个对话框，允许我们选择轮廓的数值。好的电势视图可以这样获得：关闭 auto contour 并选择 -0.5 作为第一个轮廓数值 (例如，在左边的文本输入框输入数值并按 tab 键)。单击蓝色箭头并选择"Add after"可以加入第二个轮廓值。使用数值15可以得到较好的结果。
我们现在可以管道线中单击 Colors and legends 然后选择不同的 LUT (查找表)以改变使用的颜色。让我们选择 'Paired'，因为它可以很好的分开层次。
为了得到更好的电势视图，我们想显示更多的轮廓，但是其问题是封闭的轮廓隐藏了内部情况。一个解决方案是使用剖面。右击 IsoSurface 节点，通过 "Add module" 子菜单添加 ScalarCutPlane。单击并拖动它可以移动此剖面。
为了实现 numpy 数组和可视化之间的联系，我们可以使用相同的菜单添加 Axes 和 Outline。最后，让我们添加一个颜色条。输入下面的命令:
mlab.colorbar(title='Potential', orientation='vertical')
或者使用先前使用的 LUT 中的选项。
我们想研究此电势中粒子的运动。为此我们需要推导电势陡度(gradient)相应的力。我们建立陡度函数:
def gradient(f, x, y, z, d=0.01):
""" Return the gradient of f in (x, y, z). """
fx = f(x+d, y, z)
fx_ = f(x-d, y, z)
fy = f(x, y+d, z)
fy_ = f(x, y-d, z)
fz = f(x, y, z+d)
fz_ = f(x, y, z-d)
return (fx-fx_)/(2*d), (fy-fy_)/(2*d), (fz-fz_)/(2*d)
为了检查陡度函数是否工作正常，让我们可视化它产生的向量场。为了避免显示太多的向量，我们将仅仅评定 X=50 处剖面的陡度，并且只在网格中每三个点评定一个:
Vx, Vy, Vz = gradient(V, X[50, ::3, ::3], Y[50, ::3, ::3], Z[50, ::3, ::3])
mlab.quiver3d(X[50, ::3, ::3], Y[50, ::3, ::3], Z[50, ::3, ::3],
Vx, Vy, Vz, scale_factor=-0.2, color=(1, 1, 1))
现在我们可以使用 scipy 积分轨迹。首先必须定义动态流体，这是一个这些参数和时间的函数，返回不同参数的导数。此流体被每个 ODE (常微分方程) 解算器(solver)所使用。它给出系统动态(dynamic)。我们关心的动态由电势推导出的力造成，而我们在三个方向上随时间振动此电势，作为阻尼力。调节阻尼系数及振动的数量和频率可以给出有趣的动态。 :
def flow(r, t):
""" The dynamical flow of the system """
x, y, z, vx, vy, vz = r
fx, fy, fz = gradient(V, x-.2*np.sin(6*t), y-.2*np.sin(6*t+1), z-.2*np.sin(6*t+2))
return np.array((vx, vy, vz, -fx - 0.3*vx, -fy - 0.3*vy, -fz - 0.3*vz))
现在我们可以积分轨迹:
from scipy.integrate import odeint
# Initial conditions
R0 = (0, 0, 0, 0, 0, 0)
# Times at which we want the integrator to return the positions:
t = np.linspace(0, 50, 500)
R = odeint(flow, R0, t)
右击对应的管道线节点，选择 delete 可除去剖面和向量场。然后我们开始绘制此轨迹。同时关闭相应的 Colors and legends 节点的第一个颜色条。我们用与它相关的标量信息绘制轨迹，通过调色板(colormap)显示时间:
x, y, z, vx, vy, vz = R.T
trajectory = mlab.plot3d(x, y, z, t, colormap='hot',
tube_radius=None)
mlab.colorbar(trajectory, title='Time', orientation='vertical')
没有评论:
发表评论