Python - Solution numérique d'équation différentielle Méthode d'Euler et méthode de différence centrale et méthode Rungekutta
Les solutions des équations différentielles ordinaires suivantes sont calculées par la méthode d'Euler, la méthode des différences centrales et la méthode Lungekutter.
\frac{dx}{dt}=f(t)=cos(t)
Solution générale
\begin{eqnarray}
dx&=&cos(t)dt\\
\int{dx}&=&\int{cos(t)}dt\\
x&=&sin(t)+x(t=0)
\end{eqnarray}
Méthode d'Euler
Expansion de Taylor de $ x (t + Δt) $
x(t+Δt)=x(t)+\frac{dx(t)}{dt}Δt+O(Δt^2)
Si $ O (Δt ^ 2) $ est tronqué comme une erreur de troncature
\begin{eqnarray}
x(t+Δt)&≃&x(t)+\frac{dx(t)}{dt}Δt\\
&=&x(t)+cos(t)Δt...①
\end{eqnarray}
Ensuite, la méthode de calcul séquentiel de $ x $ à $ t + Δt $ en utilisant $ x $ à $ t $.
L'algorithme est
- Donnez les valeurs initiales à $ t et x $.
- Trouvez $ x $ à $ t + Δt $ par l'équation (1).
- Répétez la plage 2 que vous souhaitez calculer.
·Code source
import math
import matplotlib.pyplot as plt
import numpy as np
def f(t):
return math.cos(t)
DELTA_T = 0.001
MAX_T = 100.0
t = 0.0 #t valeur initiale
x = 0.0 # t=X à 0
x_hist = [x]
t_hist = [t]
#Calcul séquentiel
while t < MAX_T:
x += f(t)*DELTA_T
t += DELTA_T
x_hist.append(x)
t_hist.append(t)
#Tracé de solution numérique
plt.plot(t_hist, x_hist)
#Solution exacte(sin(t))Terrain
t = np.linspace(0, MAX_T, 1/DELTA_T)
x = np.sin(t)
plt.plot(t, x)
plt.xlim(0, MAX_T)
plt.ylim(-1.3, 1.3)
plt.show()
·résultat
Δt=0.01
Δt=0.001
Méthode de différence centrale
Lorsque $ x (t + Δt) $ et $ x (t-Δt) $ sont développés par Taylor,
\begin{eqnarray}
x(t+Δt)&=&x(t)+\frac{dx(t)}{dt}Δt+\frac{1}{2}\frac{d^2x}{dt^2}Δt^2+O(Δt^3)...①\\
x(t-Δt)&=&x(t)-\frac{dx(t)}{dt}Δt+\frac{1}{2}\frac{d^2x}{dt^2}Δt^2+O(Δt^3)...②
\end{eqnarray}
①-② et tronquer $ O (Δt ^ 3) $
\begin{eqnarray}
x(t+Δt)-x(t-Δt)&≃&2\frac{dx(t)}{dt}Δt\\
x(t+Δt)&≃&x(t-Δt)+2cos(t)Δt...③\\
\end{eqnarray}
Si vous le remplacez par $ t + Δt $ ⇒ $ t $
\begin{eqnarray}
x(t)&≃&x(t-2Δt)+2cos(t-Δt)Δt...③\\
\end{eqnarray}
Lorsque $ x $ à $ t-2Δt $ est connu, la méthode de calcul séquentiel de $ x $ à $ t $ par l'équation (3).
·Code source
import math
import matplotlib.pyplot as plt
import numpy as np
def f(t):
return math.cos(t)
DELTA_T = 0.001
MAX_T = 100.0
t = 0.0 #t valeur initiale
x = 0.0 # t=X à 0
x_hist = [x]
t_hist = [t]
while t < MAX_T:
x += 2*f(t-DELTA_T)*DELTA_T
t += 2*DELTA_T
x_hist.append(x)
t_hist.append(t)
#Tracé de solution numérique
plt.plot(t_hist, x_hist)
#Solution exacte(sin(t))Terrain
t = np.linspace(0, MAX_T, 1/DELTA_T)
x = np.sin(t)
plt.plot(t, x)
plt.xlim(0, MAX_T)
plt.ylim(-1.3, 1.3)
plt.show()
·résultat
Méthode Rungekutta de 4ème ordre
\begin{eqnarray}
k_1&=&Δtf(t,x)\\
k_2&=&Δtf(t+\frac{Δt}{2}, x(t)+\frac{k_1}{2})\\
k_3&=&Δtf(t+\frac{Δt}{2}, x(t)+\frac{k_2}{2})\\
k_4&=&Δtf(t+Δt, x(t)+k_3)\\
x(t+Δt)&=&x(t)+\frac{1}{6}(k_1+2k_2+2k_3+k_4)
\end{eqnarray}
Par, $ x $ à $ t + Δt $ est calculé séquentiellement en utilisant $ x $ à $ t $.
Lorsque $ f (t, x) = f (t) $, de $ k_2 = k_3 $
\begin{eqnarray}
x(t+Δt)&=&x(t)+\frac{1}{6}(k_1+4k_2+k_4)
\end{eqnarray}
L'erreur contenue dans $ x (t + Δt) $ est $ O (Δt ^ 5) $
·Code source
import math
import matplotlib.pyplot as plt
import numpy as np
def f(t):
return math.cos(t)
DELTA_T = 0.001
MAX_T = 100.0
t = 0.0 #t valeur initiale
x = 0.0 # t=X à 0
x_hist = [x]
t_hist = [t]
#Calcul séquentiel
while t < MAX_T:
k1 = DELTA_T*f(t)
k2 = DELTA_T*f(t+DELTA_T/2)
k3 = DELTA_T*f(t+DELTA_T/2)
k4 = DELTA_T*f(t+DELTA_T)
x += (k1 + 2*k2 + 2*k3 + k4)/6
t += DELTA_T
x_hist.append(x)
t_hist.append(t)
#Tracé de solution numérique
plt.plot(t_hist, x_hist)
#Solution exacte(sin(t))Terrain
t = np.linspace(0, MAX_T, 1/DELTA_T)
x = np.sin(t)
plt.plot(t, x)
plt.xlim(0, MAX_T)
plt.ylim(-1.3, 1.3)
plt.show()
·résultat
Recommended Posts