Discussion :

[Herr_Vey]

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 . 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 :

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

[marco056]

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

Un exemple :

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
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()

[Herr_Vey]

Voilà l'équation en image :

[marco056]

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

Code :

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

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]))

[Herr_Vey]

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 :

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

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 :

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

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 )