PythonでBloch方程式を可視化する - kousuke

磁化ベクトルの微分方程式であるBloch方程式をPythonで解いてみた．

微分方程式の数値解は以下のpdfを参考に，Euler法で実装した．

http://www.isc.meiji.ac.jp/~mizutani/python/print/diffeq/differential_eq.pdf

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

ga, om, T1, T2 = 10 ** 11, 10 ** 9, 5 * 10 ** -10, 10 ** -10
H0, H1 = 10 ** -2, 5 * 10 ** -1
M0 = 1
def dx(t, x, y, z):
return(ga * (y*H0 - z*H1*np.sin(om*t)) - x/T2)
def dy(t, x, y, z):
return(ga * (z*H1*np.cos(om*t) - x*H0) - y/T2)
def dz(t, x, y, z):
return(ga * H1 * (x*np.sin(om*t) - y*np.cos(om*t)) - (z - M0)/T1)
def euler_orbit(x0, T, dt, vectorfield):
width = len(x0)
x = x0
t = 0
orbit = [] orbit.append(list(x))
while t <= T:
vx = list(map(lambda v: v(t, *x), vectorfield))
for i in range(width):
x[i] += dt * vx[i]
orbit.append(list(x)) t += dt
return(orbit)
vector = [dx, dy, dz]
x0 = [0, 0, 1]
dt = 10 ** -13
T = 10 ** -7
orbit = euler_orbit(x0, T, dt, vector)
nporbit = np.array(orbit)
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(nporbit[:,0], nporbit[:,1], nporbit[:,2])
ax.scatter3D(0, 0, 1,color='b')
ax.scatter3D(nporbit[-1,0], nporbit[-1,1], nporbit[-1,2],color='b')
ax.set_xlabel('x(t)')
ax.set_ylabel('y(t)')
ax.set_zlabel('z(t)')
plt.show()

自分なりの工夫としては，単位球を描かせて解ベクトルがその中を動いているかを確認できるようにしたことと，初期値（[0,0,1]）と終値を青点で表示したこと．前者については，MathematicaのSphereみたいに一発で書くことはできないのかな．

これは実験室座標系（固定されている座標系）で，本来は回転座標系（電磁波の周波数で回転する座標系）に移らないといけない．方程式をいじるだけなのだが，解を見るとどうも上手くいっていない．さて，どうしたものか．

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