2016-04-05 2 views
0

J'ai une fonction que je dois lui faire des dérivées partielles dépendant d'un paramètre et de l'utiliser dans une autre fonction, puis de résoudre un système ODE, La fonction qui je dois anizotropy_energy est dérivé par rapport thêta et phipython dérivées partielles, je ne peux pas l'utiliser avec

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) 

Je l'ai fait avec sympy, mais je ne peux pas utiliser ces derrivates dans odeint j'ai une fonction intermédiaire où je dois ajouter un à un structure avec thêta et phi cela dérive et ensuite utiliser cette fonction dans le programme principal avec Ode

Intermédiaire fu 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)) 

Le programme principal:

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

Ne devrait-il pas s'agir de «u_rx = sin (thêta) * cos (phi)»? – LutzL

Répondre

1

Votre erreur est anizotropy_energy.py. Vous avez des variables globales avec le même nom que les paramètres de fonction. Dans ce cas, vous avez un thêta global et les paramètres de fonction phi ET nommés theta et phi. Vous devez renommer une paire d'entre eux. J'ai étiqueté les paramètres de fonction thêta et phi à in_theta et in_phi. Vous devriez vraiment les renommer.

En outre, note que dans certains endroits que vous dites ani s otropy et certains vous disent 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)