Einführung

Eine normale Differentialgleichung ist eine Gleichung, die die unabhängigen Variablen t, die von t abhängige Funktion y (t) und ihre n-te Ableitung (n = 0,1,2, ..., N) enthält. Mit anderen Worten

g(t,y,y',y'',y^{(3)},...,y^{(N)}) = 0 \hspace{50pt}(1)

Es ist eine Gleichung, die in Form von beschrieben werden kann.

Solche Gleichungen können durch eine numerische Lösung angenähert werden, die als Eular-Methode bezeichnet wird. Ich werde es anhand eines konkreten Beispiels erklären.

Vorbereitung auf die Eular-Methode

Die obige Differentialgleichung $ (1) $ hat tatsächlich eine Form, die für die Verwendung der Euler-Methode nicht geeignet ist. Sie müssen dies irgendwie in etwas wie $ (2) $ unten umwandeln.

\frac{dx(t)}{dt} = f(t,x(t)) \hspace{50pt}(2)
\\jedoch, x(t)=\begin{pmatrix}y\\y'\\y''\\...\\y^{(N-1)}\end{pmatrix}

Aktuelles Beispiel

Stellen Sie sich ein System mit einer Masse von $ m $ an der Spitze einer Feder mit einer Federkonstante von $ k $ vor. Wenn eine externe Schwingungskraft auf das System ausgeübt wird und der Luftwiderstand berücksichtigt wird, lautet die Differentialgleichung wie folgt.

Es ist möglich, es analytisch zu lösen, aber wir werden in Betracht ziehen, es durch eine numerische Berechnungsmethode zu lösen.

m\frac{d^2y(t)}{dt^2}+2\gamma \frac{dy(t)}{dt}+ky(t)=a\sin(\omega t)\hspace{50pt}(3)

Wenn $ (3) $ transformiert wird,

\frac{d^2y(t)}{dt^2}=-\frac{2\gamma}{m} \frac{dy(t)}{dt}-\frac{k}{m}y(t)+\frac{a}{m} \sin(\omega t) \hspace{50pt}(3-1)

Ebenfalls,

\frac{dy(t)}{dt}=\frac{dy(t)}{dt}\hspace{50pt}(3-2)

Wo der Vektor $ x (t) $

x(t)=\begin{pmatrix} x_1(t)  \\ x_2(t) \end{pmatrix} = \begin{pmatrix} y(t)  \\ \frac{dy(t)}{dt} \end{pmatrix}

Wenn Sie von $ (3-1) $, $ (3-2) $,

\frac{dx(t)}{dt}
=\frac{d}{dt} \begin{pmatrix} y(t) \\ \frac{dy(t)}{dt}  \end{pmatrix}
=\begin{pmatrix} \frac{dy(t)}{dt} \\ \frac{d^2y(t)}{dt^2}  \end{pmatrix}
=\begin{pmatrix} \frac{dy(t)}{dt} \\ -\frac{2\gamma}{m} \frac{dy(t)}{dt}-\frac{k}{m}y(t)+\frac{a}{m} \sin(\omega t)  \end{pmatrix}
=\begin{pmatrix} x_2(t) \\ -\frac{2\gamma}{m} x_2(t)-\frac{k}{m}x_1(t)+\frac{a}{m} \sin(\omega t)  \end{pmatrix}

Daher kann es wie folgt ausgedrückt werden.

\frac{dx(t)}{dt} = f(t,x(t))\hspace{50pt}(4)
\\where : f(t,x)=\begin{pmatrix} x_2 \\ -\frac{2\gamma}{m} x_2-\frac{k}{m}x_1+\frac{a}{m} \sin(\omega t)  \end{pmatrix}

Eular Methode

Mit der Arbeit bis zu diesem Punkt konnten wir eine gut aussehende Differentialgleichung erstellen. Lassen Sie uns nun über die Eular-Methode selbst sprechen.

Zunächst aus der Definition der Differenzierung

\frac{dx(t)}{dt} = f(t,x(t))

⇒\lim_{\Delta t \to 0} \frac{x(t+\Delta t)-x(t)}{\Delta t} = f(t,x(t))

Wenn Sie hier $ \ Delta t $ verwenden, das klein genug ist,

\frac{x(t+\Delta t)-x(t)}{\Delta t} \fallingdotseq f(t,x(t))

⇒x(t+\Delta t) \fallingdotseq x(t)+f(t,x(t))\Delta t\hspace{50pt}(5)

Wenn Sie mit dieser Formel das erste $ x (t) $ kennen, können Sie $ x (t + \ Delta t) $ nach $ \ Delta t $ Sekunden berechnen.

Programm (Python)

Die Eular-Methode wird wie folgt in Python implementiert.

`f.py`


import numpy as np

a = 10           #a(Externe Kraftamplitude)
(m,gm2,k) = (30,5,1)#m,2Γ,k
omg = np.pi     #ω

def f(t,x):
    f0 = x[1]
    f1 = -(gm2/m)*x[1] -(k/m)*x[0] + (a/m)*np.sin(omg*t)
    return np.array([f0,f1])

`f_D.py`


import numpy as np
from f import f as f

dt = 0.01

def f_D_eular(t,x):
    return (t + dt, x + f(t,x)*dt)

`main.py`


import numpy as np
from f_D import f_D_eular as f_D

def main():
    (t0,x0)=(0,np.array([1,0]))
    t1 = 100
    (t,x)=(t0,x0)
    while(t < t1):
       print(t,x[0])
       (t,x)=f_D(t,x)

main()

Es wird so sein. Übrigens ist der Graph des Ausführungsergebnisses dieses Programms wie folgt.