[PYTHON] Finden Sie die numerische Lösung der gewöhnlichen Differentialgleichung zweiter Ordnung mit scipy

Einführung

Im Allgemeinen ist die Bewegungsgleichung eine Differentialgleichung zweiter Ordnung für die Zeit des Positionsvektors. Sie müssen diese also lösen, um die Bewegung zu simulieren. Es gibt jedoch Fälle, in denen eine allgemeine Lösung nicht einfach zu erhalten ist. In diesem Fall werden numerische Berechnungen durchgeführt. In einem solchen Fall ist die Methode zum Finden der numerischen Lösung der Differentialgleichung in einem Modul namens ode in SciPy implementiert. Berechnen wir also damit.

scipy.integrate.ode http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html

Problemstellung

Betrachten wir der Einfachheit halber das Problem des freien Falls aus der Höhe "h" für die Qualitätspunkte, die durch Masse "m" und Position "x" beschrieben werden. Hier ist der Luftwiderstand vernachlässigbar und die auf den Qualitätspunkt ausgeübte Kraft ist nur die Schwerkraft.

Bewegungsgleichung

Aus der obigen Annahme ergibt sich die Bewegungsgleichung

m \ddot{x} = - g

Es wird sein. "G" ist jedoch eine Schwerkraftkonstante. In der Physikkonvention repräsentieren die Punkte über der Variablen auch die Zeitdifferenzierung. Wenn zwei Punkte vorhanden sind, handelt es sich um ein Differential zweiter Ordnung.

Nun, dies kann eindeutig eine analytische Lösung finden,

x = - \frac{1}{2} g t^2 + h

Wenn jedoch die Zeit in diese eingefügt wird, gibt es keine Quelle oder untergeordnetes Element, sodass sie nur zum Vergleich mit der numerischen Lösung verwendet wird.

Formelumwandlung, damit es mit scipy gelöst werden kann

Nun, das ist das Hauptthema. Da das Ode-Modul von SciPy (wahrscheinlich wie jedes andere) nur normale Differentialgleichungen erster Ordnung lösen kann, gibt es die folgenden Variablentransformationen:

v := \dot{x}

Daher ist die Bewegungsgleichung

\dot{x} = v
\dot{v} = - \frac{g}{m}

Es wird sein. Hier wird die Differenzierung nach links verschoben, dies ist jedoch keine Frage des Zufalls oder des Aussehens, sondern eine Einschränkung des Odenmoduls. Die lineare Differentialgleichung sollte im Allgemeinen die folgende Form haben.

\frac{\mathrm{d}}{\mathrm{d}t}
\begin{bmatrix}
x_1 \\
\vdots \\
x_n
\end{bmatrix} =
\begin{bmatrix}
a_{11} & a_{12} & \dots & a_{1n} \\
\vdots & \ddots &  & \vdots \\
a_{n1} & \dots &  & a_{nn}
\end{bmatrix}
\begin{bmatrix}
x_1 \\
\vdots \\
x_n
\end{bmatrix} +
\begin{bmatrix}
b_1 \\
\vdots \\
b_n
\end{bmatrix}

In dieser Angelegenheit,

\frac{\mathrm{d}}{\mathrm{d}t}
\begin{bmatrix}
x_1 \\
x_2
\end{bmatrix} =
\begin{bmatrix}
0 & 1 \\
0 & 0
\end{bmatrix}
\begin{bmatrix}
x_1 \\
x_2
\end{bmatrix} +
\begin{bmatrix}
0 \\
- \frac{g}{m}
\end{bmatrix}

ist. ** In dem unten gezeigten Code erscheinen "erste Formel" und "zweite Formel" in der erweiterten Formel **. Hier ergeben sich die folgenden Anfangsbedingungen aus den Beschreibungen "aus der Höhe fallen" h "und" frei fallen ".

x_1(0) = h
x_2(0) = 0

Nachdem beide die Form von normalen Differentialgleichungen erster Ordnung haben, können sie schließlich mit dem Ode-Modul gelöst werden.

Löse normale Differentialgleichungen mit scipy

Lassen Sie uns dies nun im Python-Code tun. Dies kann jedoch fast mechanisch erfolgen. Hier ist der vollständige Code.

#-*- coding:utf-8 -*-

import numpy as np
from scipy.integrate import odeint

g = 9.8     #Schwerkraftkonstante
m = 1.0     #Masse
h = 10      #Ausgangsposition

def f(x, t):
    ret = [
        x[1],      #Rechte Seite von Formel 1
        -g / m     #Rechte Seite von Formel 2
    ]
    return ret


def main():
    #Ausgangszustand
    x0 = [
        h,    #Anfangsbedingungen der ersten Gleichung
        0     #Anfangsbedingung der zweiten Gleichung
    ]

    #Intervall zu berechnen
    #Die Argumente sind in der Reihenfolge "Startzeit", "Endzeit" und "Inkrement".
    t = np.arange(0, 10, 0.1)

    #Integrieren
    x = odeint(f, x0, t)

    #Zeigen Sie das Ergebnis an (drucken Sie es vorerst aus)
    print(x)


if __name__ == '__main__':
    main()

Wenn Sie dies tun, erhalten Sie eine Ausgabe ähnlich der folgenden:

$ python free_fall_sample.py
[[  1.00000000e+01   0.00000000e+00]
 [  9.95100000e+00  -9.80000000e-01]
 [  9.80400000e+00  -1.96000000e+00]
 [  9.55900000e+00  -2.94000000e+00]
...
 [ -4.60596000e+02  -9.60400000e+01]
 [ -4.70249000e+02  -9.70200000e+01]]

Es wird in Form eines Doppelarrays angegeben und speichert für jedes Mal eine Reihe von Lösungen. Die Menge der Lösungen besteht darin, dass das erste Element die Lösung der ersten Gleichung und das zweite Element die Lösung der zweiten Gleichung ist (kurz gesagt, wie Sie es sehen).

Als Test wird die numerische Lösung bei "t = 10" mit dem aus der analytischen Lösung berechneten Wert verglichen.

Numerische Lösung= -4.70249000e+02 = -470.2 ...
Analytische Lösung= -9.8 / 2 * (10 * 10) + 10 = -480.0 ...

Der Fehler beträgt ca. 10 [m](* Herr Koryor wies darauf hin, dass er derzeit korrigiert wird. Weitere Informationen finden Sie im Kommentarbereich). Es gibt verschiedene Möglichkeiten, diese Genauigkeit zu verbessern. Am einfachsten ist es jedoch, die Schrittweite zu verringern. Im Code ist dies der Teil:

    #Intervall zu berechnen
    #Die Argumente sind in der Reihenfolge "Startzeit", "Endzeit" und "Inkrement".
    t = np.arange(0, 10, 0.01) # <=Veränderung

Mit dieser Änderung beträgt die numerische Lösung 479 und der Fehler etwa 1,0.

Numerische Lösung= -4.79020490e+02 = -479.0 ...
Analytische Lösung= -9.8 / 2 * (10 * 10) + 10 = -480.0 ...

Wenn das Reduzieren der Schrittgröße keine ausreichende Genauigkeit bietet oder zu lange dauert, ziehen Sie einen anderen Algorithmus in Betracht.

Erwägung

Übrigens, wenn Sie sich die Lösung ansehen, die ernsthaft herauskam, können Sie sehen, dass die Position der Plakette negativ ist. Dies ist ein Zustand, in dem Studenten, die nach Physik streben, einmal Erfahrungen machen und die Zusagen in den Boden versenkt werden. Da die Randbedingungen und die Bodenabstoßung ignoriert werden, durchdringt die Plakette den Boden und beschleunigt das Eindringen von der anderen Seite der Erde. Um dies jedoch zu verhindern, treten Luftwiderstand und Bodenabstoßung auf Muss definiert werden. Dies ist ein einfacher Fall, sodass Sie ihn sofort bemerken werden. Wenn es sich jedoch um einen komplizierten Fall handelt, können Sie ihn leicht übersehen. Vergessen Sie also nicht, die Ergebnisse ruhig zu betrachten, sobald Sie ihn gelöst haben.

Zusammenfassung

Ich fand die numerische Lösung der gewöhnlichen Differentialgleichung mit scipy. Sobald Sie gelernt haben, wie man es benutzt, ist es sehr einfach und unkompliziert. Denken Sie also daran, wenn Sie die Bewegungsgleichung lösen möchten.

Recommended Posts

Finden Sie die numerische Lösung der gewöhnlichen Differentialgleichung zweiter Ordnung mit scipy
Finden Sie die Lösung der Gleichung n-ter Ordnung mit Python
[Wissenschaftlich-technische Berechnung mit Python] Numerische Lösung der gewöhnlichen Differentialgleichung zweiter Ordnung, Anfangswertproblem, numerische Berechnung
Numerische Analyse gewöhnlicher Differentialgleichungen mit Scipys Odeint und Ode
Lösung des Anfangswertproblems gewöhnlicher Differentialgleichungen mit JModelica
Die Geschichte der numerischen Berechnung von Differentialgleichungen mit TensorFlow 2.0
Finden Sie den Tag nach Datum / Uhrzeit heraus
[Wissenschaftlich-technische Berechnung von Python] Numerische Berechnung zur Ermittlung des Ableitungswerts (Differential)
[Wissenschaftlich-technische Berechnung mit Python] Analytische Lösungssympathie zur Lösung von Gleichungen
Berechnen Sie die Summe der eindeutigen Werte durch Pandas-Kreuztabellen
Finden Sie den Speicherort der mit pip installierten Pakete heraus
[Wissenschaftlich-technische Berechnung mit Python] Numerische Lösung gewöhnlicher Differentialgleichungen erster Ordnung, Anfangswertproblem, numerische Berechnung
Ich habe versucht, die Entropie des Bildes mit Python zu finden
Python - Differentialgleichung Numerische Lösung Euler-Methode & Zentrale Differenzmethode & Rungekutta-Methode
Sehen Sie, wie schnell Sie mit NumPy / SciPy beschleunigen können
Bewegungsgleichung mit Sympy
Finden Sie den optimalen Wert der Funktion mit einem genetischen Algorithmus (Teil 2)
Numerische Berechnung partieller Differentialgleichungen mit Singularität (am Beispiel der thermischen Gleichungen vom Hardy-Hénon-Typ)
Finden Sie die Definition des Wertes von errno
Finden Sie die Bearbeitungsentfernung (Levenshtein-Entfernung) mit Python
Ich habe die numerische Berechnung von Python durch Rust ersetzt und die Geschwindigkeit verglichen
Ermitteln Sie mit NumPy die Trägheitsspindel und das Hauptträgheitsmoment aus dem Trägheitstensor
Finden Sie die allgemeinen Begriffe der Tribonacci-Sequenz in linearer Algebra und Python
[Wissenschaftlich-technische Berechnung durch Python] Lösung des Randwertproblems gewöhnlicher Differentialgleichungen im Matrixformat, numerische Berechnung