Python - Differentialgleichung Numerische Lösung Euler-Methode & Zentrale Differenzmethode & Rungekutta-Methode

Die Lösungen der folgenden gewöhnlichen Differentialgleichungen werden nach der Euler-Methode, der zentralen Differenzmethode und der Lungekutter-Methode berechnet.

\frac{dx}{dt}=f(t)=cos(t)

Allgemeine Lösung

\begin{eqnarray}
dx&=&cos(t)dt\\
\int{dx}&=&\int{cos(t)}dt\\
x&=&sin(t)+x(t=0)
\end{eqnarray}

Euler-Methode

Taylor-Expansion von $ x (t + Δt) $

x(t+Δt)=x(t)+\frac{dx(t)}{dt}Δt+O(Δt^2)

Wenn $ O (Δt ^ 2) $ als Kürzungsfehler abgeschnitten wird

\begin{eqnarray}
x(t+Δt)&≃&x(t)+\frac{dx(t)}{dt}Δt\\
&=&x(t)+cos(t)Δt...①
\end{eqnarray}

Dann die Methode der sequentiellen Berechnung von $ x $ bei $ t + Δt $ unter Verwendung von $ x $ bei $ t $.

Der Algorithmus ist

1. Geben Sie $ t und x $ Anfangswerte.
2. Finden Sie $ x $ bei $ t + Δt $ nach Gleichung (1).
3. Wiederholen Sie den zu berechnenden Bereich 2.

·Quellcode

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 Anfangswert
x = 0.0 # t=X bei 0

x_hist = [x]
t_hist = [t]

#Sequenzielle Berechnung
while t < MAX_T:
    x += f(t)*DELTA_T
    t += DELTA_T
    x_hist.append(x)
    t_hist.append(t)

#Numerisches Lösungsdiagramm
plt.plot(t_hist, x_hist)

#Genaue Lösung(sin(t))Handlung
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()

·Ergebnis

Δt=0.01

Δt=0.001

Zentrale Differenzmethode

Wenn $ x (t + Δt) $ und $ x (t-Δt) $ von Taylor erweitert werden,

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

①-② und Abschneiden von $ 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}

Wenn Sie es durch $ t + Δt $ ⇒ $ t $ ersetzen

\begin{eqnarray}
x(t)&≃&x(t-2Δt)+2cos(t-Δt)Δt...③\\
\end{eqnarray}

Wenn $ x $ bei $ t-2Δt $ bekannt ist, wird das Verfahren zum sequentiellen Berechnen von $ x $ bei $ t $ durch Gleichung (3) durchgeführt.

·Quellcode

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 Anfangswert
x = 0.0 # t=X bei 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)

#Numerisches Lösungsdiagramm
plt.plot(t_hist, x_hist)

#Genaue Lösung(sin(t))Handlung
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()

·Ergebnis

Δt=0.01

Δt=0.001

Rungekutta-Methode 4. Ordnung

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

Mit wird $ x $ bei $ t + Δt $ nacheinander unter Verwendung von $ x $ bei $ t $ berechnet.

Wenn $ f (t, x) = f (t) $, von $ k_2 = k_3 $

\begin{eqnarray}
x(t+Δt)&=&x(t)+\frac{1}{6}(k_1+4k_2+k_4)
\end{eqnarray}

Der in $ x (t + Δt) $ enthaltene Fehler ist $ O (Δt ^ 5) $

·Quellcode

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 Anfangswert
x = 0.0 # t=X bei 0

x_hist = [x]
t_hist = [t]

#Sequenzielle Berechnung
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)
    
#Numerisches Lösungsdiagramm
plt.plot(t_hist, x_hist)

#Genaue Lösung(sin(t))Handlung
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()

·Ergebnis

Δt=0.001

Δt=0.001