2016-04-05 15 views
0

私は、パラメータに依存する偏微分にして、これを別の関数で使用してから、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() 
+0

'u_rx = sin(theta)* cos(phi)'ではないはずですか? – LutzL

答えて

1

あなたのエラーがanizotropy_energy.pyです。関数のパラメータと同じ名前のグローバル変数があります。 この場合、thetaとphiという名前のグローバルなthetaとphi AND関数のパラメータがあります。それらのペアの名前を変更する必要があります。 関数パラメータthetaとphiをin_thetaとin_phiにラベルを付けました。あなたは本当にそれらの名前を変更する必要があります。いくつかの場所で使用すると、ANI otropy、あなたはANI Z otropyを言ういくつかを言うこと

また

注意

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(in_theta, in_phi, theta_k, phi_k): 

    u_rx = sin(in_theta)*sin(in_phi) 
    u_ry = sin(in_theta)*sin(in_phi) 
    u_rz = cos(in_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 

    return diff((u_kx*u_rx + u_ky*u_ry + u_kz*u_rz)**2, theta) 
関連する問題