Trajectoire particule chargée et ODEINT:Excess work done on this call (perhaps wrong Dfun type)

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)

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58

# -*- 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

Code:

---

1 2 3 4 5 6 7

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)

---

J'ai ensuite simplifié le code : ce qui ne change rien d'ailleurs au problème. Les appels aux fonctions n'ont pas l'air d'être le problème.

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

# -*- 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=1;m=1;alpha=q/m   # le champ E=np.array([0,0,0]);B=np.array([0,0,1]);#composantes sur x,y,z   #les CI X0=np.array([0.0,0.0,0.0,100.0,0.0,0.0]);v0=X0[3:6]#X0([x0,y0,z0,vx0,vy0,vz0])   # la durée t=np.linspace(0,15,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 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]; dX[3:6]=alpha*(E+vectoriel(v,B)); # print('a= ',a) return dX     X=odeint(Xprim,X0,t); pypl.plot(X[:,1],X[:,0]) pypl.grid() pypl.show()

---

Et voilà que jusqu'à tmax=15 no pb et à 16 boum, ça revient à l'origine... C'est vraiment un problème d'info qui dépasse le débutant que je suis.