PythonでBloch方程式を可視化する②〜odeint編〜

プログラム Python

前回（https://kousukejapan.hatenablog.jp/entry/2019/03/16/222053）はEuler法でBloch方程式を解いたが，odeintを使うことを学んだのでそっちでもやってみた．

 

import numpy as np
# from math import sin, cos
from scipy.integrate import odeint
from matplotlib import cm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

M0 = 1
def func(v, t, ga, om, H0, H1, T1, T2):
    return [ga * (v[1] * H0 - v[2] * H1 * np.sin(om * t)) - v[0] / T2,
    ga * (v[2] * H1 *np.cos(om * t) - v[0] * H0) - v[1] / T2,
    ga * H1 * (v[0] *np.sin(om * t) - v[1] *np.cos(om * t)) - (v[2] - M0) / T1]
ga, om, T1, T2 = 10 ** 11, 10 ** 9, 5 * 10 ** -10, 10 ** -10
H0, H1 = 10 ** -2, 5 * 10 ** -1
v0 = [0, 0, 1]
t = np.arange(0, 10 ** -8, 10 ** -13)
v = odeint(func, v0, t, args=(ga, om, H0, H1, T1, T2))

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
U = np.linspace(0, 2 * np.pi, 100)
V = np.linspace(0, np.pi, 100)
X = np.outer(np.cos(U), np.sin(V))
Y = np.outer(np.sin(U), np.sin(V))
Z = np.outer(np.ones(np.size(U)), np.cos(V))
ax.plot_surface(X, Y, Z, cmap=cm.coolwarm, alpha = 0.5)

ax = fig.gca(projection='3d')
ax.plot(v[:, 0], v[:, 1], v[:, 2])
ax.scatter3D(0, 0, 1,color='b')
ax.scatter3D(v[-1,0], v[-1,1], v[-1,2],color='b')
ax.set_xlabel('x(t)')
ax.set_ylabel('y(t)')
ax.set_zlabel('z(t)')

 

forループがなくなった分スッキリした気がする．

精度はどうなんだろう？

 

2019-03-25追記
numpyに三角関数が入っているから，mathモジュールは要らなかった．