Méthode d'Euler vs Odeint

Bonjour,

Je cherche à résoudre numériquement une équation différentielle du second ordre non linéaire que je ne vais pas écrire ne sachant pas comment le faire proprement :roll:. J'ai d'abord utilisé la méthode d'Euler mais je n'ai pas obtenu la courbe voulue(celle d'un oscillateur amorti). J'ai donc essayé par curiosité la fonction odeint de scipy.integrate et là j'obtiens la courbe voulue. Je ne comprends pas d'où vient une différence aussi énorme alors que je devrais obtenir à peu près la même chose pour les deux... Voilà le code :

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

import numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint     rho = 1.204 C = 0.44 v = 10 S = 0.0028 s = 0.0005 Mb = 0.002 Ma = 0.003 l = 0.02     CI = [np.pi, 28] n=1000 t = np.linspace(0,5,n) h=5/n   # Avec Odeint def fun(z,t): return np.array([z[1],(-1/2)*rho*C*v*S*((z[1]/(Mb*(1+(Mb/Ma)))+(np.sin(z[0])*v/(l*(Ma+Mb)))))]) sol=odeint(fun,CI,t)[:,0]   #Il s'agit d'un angle donc je le ramène dans [-pi,pi] for i in range(len(sol)): while sol[i]<-np.pi or sol[i] > np.pi: if sol[i]>np.pi: sol[i]-=2*np.pi else: sol[i]+=2*np.pi   #Avec Euler   def fun2(y,yd): return (-1/2)*rho*C*v*S*((yd/(Mb*(1+(Mb/Ma)))+(np.sin(y)*v/(l*(Ma+Mb)))))   L1 = [-np.pi] L2=[28] for i in range(1,n): L1.append(L1[i-1] + h*L2[i-1]) L2.append(L2[i-1] + h*fun2(L1[i-1],L2[i-1]))   #idem à ci-dessus for i in range(n): while L1[i]< -np.pi or L1[i] > np.pi: if L1[i]>np.pi: L1[i]-=2*np.pi else: L1[i]+=2*np.pi   plt.plot(t,sol) plt.plot(t,L1) plt.show()

---

Je vous remercie d'avance pour vos réponses

Difficile de corriger ta méthode d'Euler sans l'équation différentielle.

Un exemple :

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 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90

import numpy as np import matplotlib.pyplot as plt # Outils graphiques   R=1e3 # Resistance du Pont de Wien C=100e-9 # Capacite du Pont de Wien K=3.2 # Gain de l'etage d'amplification tfin=15e-3 # instant final de la simulation dt=0.01e-3 # pas vsat=15 # tension de saturation de l'AO   # Conditions Initiales t=[0.0] e=[0.01] s=[0.0] ds=[0.0] de=[0.01]     while t[-1]<tfin: de.append((1/(R*C)*ds[-1]-3/(R*C)*de[-1]-e[-1]/(R*R*C*C))*dt+de[-1]) e.append(de[-1]*dt+e[-1]) if e[-1]<-vsat/K: s.append(-vsat) elif e[-1]>vsat/K: s.append(vsat) else: s.append(K*e[-1]) ds.append((s[-1]-s[-2])/dt) t.append(t[-1]+dt)   plt.clf() plt.plot(t,e) plt.axis([-0.0005,0.015,-6,6]) plt.xlabel("temps t (s)") plt.ylabel("Entrée $v_e$ (V)") plt.title("Demarrage des oscillations") plt.plot(0,0.01,'ro') plt.savefig("wien-demarrage-ve.eps") plt.show()   plt.clf() plt.plot(e,de) plt.xlabel("Entrée $v_e$ (V)") plt.ylabel("Dérivée d$v_e$/dt (V/s)") plt.title("Portrait de phase") plt.plot(0.01,0.01,'ro') plt.savefig("wien-portrait-ve.eps") plt.show()   de_ar=np.array(de,float) de_red=de_ar*(R*C) plt.axis('equal') plt.plot(e,de_red) plt.xlabel("Entrée $v_e$ (V)") plt.ylabel("RC * d$v_e$/dt (V/s)") plt.title("Portrait de phase réduit") plt.plot(0.01,0.01,'ro') plt.savefig("wien-portrait-red-ve.eps") plt.show()   plt.clf() plt.plot(t,s) plt.axis([-0.0005,0.015,-16,16]) plt.xlabel("temps t (s)") plt.ylabel("Sortie $v_s$ (V)") plt.title("Demarrage des oscillations") plt.plot(0,0.01,'ro') plt.savefig("wien-demarrage-vs.eps") plt.show()   plt.clf() plt.plot(s,ds) plt.axis([-16,16,-200000,200000]) plt.xlabel("Sortie $v_s$ (V)") plt.ylabel("Dérivée d$v_s$/dt (V/s)") plt.title("Portrait de phase") plt.plot(0.01,0.01,'ro') plt.savefig("wien-portrait-vs.eps") plt.show()   ds_ar=np.array(ds,float) ds_red=ds_ar*(R*C) plt.axis('equal') plt.plot(s,ds_red) plt.xlabel("Sortie $v_s$ (V)") plt.ylabel("RC * d$v_s$/dt (V/s)") plt.title("Portrait de phase réduit") plt.plot(0.01,0.01,'ro') plt.savefig("wien-portrait-red-vs.eps") plt.show()

---

Voilà l'équation en image :

Pièce jointe 209273

Pas trop le temps de creuser en ce moment mais je suis sceptique concernant ces lignes là :

Code:

---

1 2 3

for i in range(1,n): L1.append(L1[i-1] + h*L2[i-1]) L2.append(L2[i-1] + h*fun2(L1[i-1],L2[i-1]))

---

J'ai repris la syntaxe utilisée dans votre code pour remplacer les lignes qui vous semblaient douteuses et ça fonctionne. Merci beaucoup !
Pour les intéressés j'ai remplacé

Code:

---

1 2 3 4 5 6

L1 = [-np.pi] L2=[28] for i in range(1,n): L1.append(L1[i-1] + h*L2[i-1]) L2.append(L2[i-1] + h*fun2(L1[i-1],L2[i-1]))

---

par

Code:

---

1 2 3 4 5 6 7 8 9 10

tps = [ti] L1 = [-np.pi] L2=[28] ti = 0 tf = 5 while tps[-1] < tf: L1.append(L1[-1] + h*L2[-1]) L2.append(L2[-1] + h*fun2(L1[-1],L2[-1])) tps.append(tps[-1] + h )

---