scipyでMayaviを使う

このチュートリアルの例では， scipy を使用して数値計算を行いながら，Mayaviをインタラクティブに使用して numpy を表示する方法を示します．ここでは，Pythonの数値ツールに精通していることを前提に，これらのツールと組み合わせてMayaviを使用する方法について説明します．

ポテンシャル中の粒子の軌道を調べてみましょう．これは物理学と工学において非常に一般的な問題であり，ポテンシャルと軌跡の可視化は問題の理解を発展させる鍵です．

興味があるポテンシャルは,放物線状閉込めに浸した周期的格子である．このポテンシャルを振って，パーティクルがラティスの穴から別の穴にジャンプする様子を確認します．放物線状の閉じ込めは，パーティクルのエクスカーションを制限するために存在します．

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)

ポテンシャルを定義したので，3Dでどのように見えるかを確認します．これを行うには，点の3Dグリッドを作成し，次の点でサンプルします．

X, Y, Z = np.mgrid[-2:2:100j, -2:2:100j, -2:2:100j]
V(X, Y, Z)

このボリュームデータをインタラクティブに視覚化するには，( mlab: 3Dプロット用のPythonスクリプト を参照してください．)モジュールを使用します．そのためには，Mayavi2アプリケーションの組み込みシェルを使用するか，または ipython --gui=qt でコマンドをインタラクティブなPythonシェルに入力することをお勧めします．ポテンシャルの3D等値面を表示します．

from mayavi import mlab
mlab.contour3d(X, Y, Z, V)

上記のコマンドで作成されたビジュアライゼーションは，ビューを回転することで操作できますが，ポテンシャルの構造を十分に理解するには，等値面を変更すると便利です．これを行うには，Mayaviパイプラインツリーの IsoSurface をダブルクリックします．(ipythonから実行している場合は，シーン上のMayaviアイコンをクリックしてパイプラインをポップアップする必要があります．)ダイアログボックスが開き，使用する等高線の値を選択できます．自動コンターをオフにして，コンターの最初の値として-0.5を選択すると，電位を適切に表示できます(例えば，右側のテキストボックスに入力して tab キーを押します．)．青い矢印をクリックして "Add after" を選択すると，別の輪郭を追加できます．値15を使用すると良い結果が得られます．

警告

"Add after" UIは以前，バージョン3.2 .0(既知のバグは 既知のバグと問題 を参照してください．)までのMayaviでwxPythonバックエンドを持つLinux上でクラッシュしていました．

パイプラインの Colors and legends をクリックし，別のLUT(テーブルの検索)を選択して使用するカラーを変更します．レベルを分ける 'Paired' を選択してみましょう．

_images/example_potential_ipython.jpg

可能性をよりよく把握するために，より多くの等高線を表示したいと考えていますが，この方法の問題は，閉じた等高線によって内部が隠されることです．解決策としては，切断面を使用します． IsoSurface ノードを右クリックし， "Add module" サブメニューから ScalarCutPlane を追加します．切断面を移動するには，その切断面をクリックしてドラッグします．

numpy配列と視覚化の間のリンクを作成するには，同じメニューを使用して軸とアウトラインを追加します．最後に，カラーバーを追加します．次のように入力します．

mlab.colorbar(title='Potential', orientation='vertical')

または，前述のLUTダイアログボックスのオプションを使用します．

_images/example_potential.jpg

このポテンシャルにおける粒子の運動を研究します．このためには，ポテンシャルの勾配によって与えられる，対応する力を導く必要があります．グラデーション関数を作成します．

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))
_images/example_gradient.jpg

これで scipy を使って軌道を統合することができます．最初に動的フローを定義する必要があります．これは，異なるパラメータの導関数を，これらのパラメータと時間の関数として返す関数です．このフローはすべての ODE (常微分方程式)ソルバーで使用され，システムの動的特性を提供します．我々が興味を持っている動力学は,減衰力と同様に,3方向で時間とともに振動するポテンシャルから生じる力で構成されています．減衰係数と振動の量と頻度は興味深い動的を与えるように調整しました．

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)

対応するパイプラインノードを右クリックして削除を選択することで，切断平面とベクトルフィールドを削除した後で軌道をプロットできます．また，対応する Colors and legends ノードの最初のカラーバーをオフにします．追加のスカラー情報をアタッチして軌道をプロットし，カラーマップで時間を表示します．

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')
_images/example_trajectories.jpg