私は、パラメータに依存する偏微分にして、これを別の関数で使用してから、ODEシステムを解く関数を持っています。誰が私が派生する必要があり、私はsympyでそれを作ったシータとphiPythonの偏微分、数値では使用できません
import sympy
import numpy as np
from sympy import Symbol, diff, sin, cos
theta = Symbol('theta')
phi = Symbol('phi')
theta_k = Symbol('theta_k')
phi_k = Symbol('phi_k')
def anizotropy_energy(theta, phi, theta_k, phi_k):
u_rx = sin(theta)*sin(phi)
u_ry = sin(theta)*sin(phi)
u_rz = cos(theta)
u_kx = sin(theta_k)*cos(phi_k)
u_ky = sin(theta_k)*sin(phi_k)
u_kz = cos(theta_k)
u_kx*u_rx + u_ky*u_ry + u_kz*u_rz
diff((u_kx*u_rx + u_ky*u_ry + u_kz*u_rz)**2, theta)
のに関してanizotropy_energyですが、私は私が私が追加する必要が中間の機能を持っているodeint でこれらのderrivatesを使用することはできませんthetaとphiを持つ構造体が完成し、この関数をメインプログラムで使用します。
中級nction:
import math
import numpy as np
from anisotropy_energy import anizotropy_energy
def thetafunc_anisotropy(alpha,theta,phi, hx, hy, hz, theta_k, phi_k):
return alpha*(hx*np.cos(theta)*np.cos(phi) + hy*np.cos(theta)*np.sin(phi)- \
hz*np.sin(theta) + \
anizotropy_energy(theta, phi, theta_k, phi_k) + \
(-hx*np.sin(phi) + hy*np.cos(phi))
# diff(anizotropy_energy(theta, phi, theta_k, phi_k), phi))
メインプログラム:
import matplotlib as mpl
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import thetafunc_anizotropy
from thetafunc_anizotropy import thetafunc_anisotropy
import phifunc_anizotropy
from phifunc_anizotropy import phifunc_anisotropy
import sympy
def LLG(y, t, alpha, hx, hy, hz, theta_k, phi_k):
theta, phi = y
dydt = [thetafunc_anisotropy(alpha,theta,phi, hx, hy, hz, theta_k, phi_k), thetafunc_anisotropy(alpha,theta,phi, hx, hy, hz, theta_k, phi_k)]
return dydt
alpha = 0.1
H = 1.0
t0 = 60.0*np.pi/180.0
p0 = 0.0*np.pi/180.0
hx = H*np.sin(t0)*np.cos(p0)
hy = H*np.sin(t0)*np.sin(p0)
hz = H*np.cos(t0)
theta_k = 60.0*np.pi/180.0
phi_k = 60.0*np.pi/180.0
y0 = [np.pi, -0*np.pi]
t = np.linspace(0, 1000, 10000)
sol = odeint(LLG, y0, t, args=(alpha,hx, hy, hz, theta_k, phi_k))
print sol
mpl.rcParams['legend.fontsize'] = 10
#
fig = plt.figure()
ax = fig.gca(projection='3d')
#x = np.sol[:, 0]
#y = sol[:, 1]
#ax.plot(x, y, label='parametric curve')
#ax.legend()
x = np.sin(sol[:, 0])*np.cos(sol[:, 1])
y = np.sin(sol[:, 0])*np.sin(sol[:, 1])
z = np.cos(sol[:, 0])
ax.plot(x, y, z, label='parametric curve')
ax.legend()
#plt.show()
#plt.axes(projection='3d')
#plt.plot(t, sol[:, 0], 'b', label='$\\theta(t)$')
#plt.plot(t,sol[:,1], 'r', label='$\\varphi(t)$')
#
#plt.legend(loc='best')
#plt.xlabel('t')
#
#plt.grid()
#plt.show()
'u_rx = sin(theta)* cos(phi)'ではないはずですか? – LutzL