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