Je débute Python que je souhaite intégrer dans mes TD de physique.
Notamment l'usage de ODEINT.
J'ai réussi à intégrer celui-ci pour les oscillateurs amortis.
J'ai tenté aujourd'hui de faire de même pour le mouvement d'une particule chargée dans un champ électromagnétique.
C'est une catastrophe .
Le mouvement circulaire n'est correctement décrit que pour certaines valeurs de m et de q ainsi que pour l'intervalle de temps d'étude.
Alors, évidemment, le mouvement en 3D n'est pas pour tout de suite.
Je me demande si ce n'est pas l'usage des tableaux qui est inutile : je croyais qu'ils augmentaient la vitesse de calcul.
Ou si c'est mon appel de fonctions l'une dans l'autre à chaque itération qui pose problème : Func dans OE appelle force qui appelle produit vectoriel…
Ou un usage malheureux des variables.
Ou une stupidité énorme que je ne vois pas…
D'avance un grand merci. Je vais bien entendu continuer à chercher seul en attendant ;o)
# -- coding: utf-8 --
"""
Created on Wed Jan 4 10:33:52 2017
@author: seb
"""
import matplotlib.pyplot as pypl
import numpy as np
from scipy . integrate import odeint
#constantes de la particule chargée
q=0.1;m=0.1;
# le champ
E=np.array([0,0,0]);B=np.array([0,0,1]);
#les CI
X0=np.array([0.0,0.0,0.0,10.0,0.0,0.0]);v0=X0[3:6]#X0([x0,y0,z0,vx0,vy0,vz0])
# la durée
t=np.linspace(0,10,100)
def vectoriel(v,w):
z=np.array([0,0,0]) #définition de z
z[0]=v[1]*w[2]-v[2]*w[1];
z[1]=v[2]*w[0]-v[0]*w[2];
z[2]=v[0]*w[1]-v[1]*w[0];
# produit=np.array([v,w,z]) ssi on veut tout dans un tableau
return z
def force(q,E,v,B):
f=np.array([0,0,0]);
w=vectoriel(v,B);
#print ('v^B= ',w);#pour vérifier le produit vectoriel
f=q*(E+w);#on économise une multiplication ici
return f
def Xprim(X,t):
dX=np.array([0.0,0.0,0.0,0.0,0.0,0.0]);F=np.array([0,0,0]);
v=X[3:6]
#dX[dx_dt,dy_dt,dz_dt,d²x_dt²,d²y_dt²,d²z_dt²]
dX[0]=X[3];dX[1]=X[4];dX[2]=X[5];
a=force(q,E,v,B)/m;
dX[3:6]=a[:];
# print('a= ',a)
return dX
#print('force initiale=',force(q,E,v0,B))#pour vérifier le calcul de la force
#X=2*X0
#dX=Xprim(X,t)
#print('a ='+str(dX[3:6]))#pour vérifier le calcul de l'accélération
X=odeint(Xprim,X0,t);
pypl.plot(X[:,1],X[:,0])
pypl.grid()
pypl.show()
Ce qui me renvoie en changeant q en 1,6e-19 et m=9.1e-31 l'erreur suivante
seb@seb-desktop:~$ python qvB.py
lsoda-- at current t (=r1), mxstep (=i1) steps
taken on this call before reaching tout
in above message, i1 = 500
in above message, r1 = 0.4604894560110D-11
/usr/lib/python2.7/dist-packages/scipy/integrate/odepack.py:218: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
warnings.warn(warning_msg, ODEintWarning)
Suivre le flux des commentaires
Note : les commentaires appartiennent à celles et ceux qui les ont postés. Nous n’en sommes pas responsables.