Discussion :

[sebfoy]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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.