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()
Ne devrait-il pas s'agir de «u_rx = sin (thêta) * cos (phi)»? – LutzL