Imitieren die Funktion ode45() in Python

Die Funktion ode45() in MATLAB löst Anfangswertprobleme für gewöhnliche Differentialgleichungen (ODEs) mit einem adaptiven Runge-Kutta-Verfahren. Python hat keine Funktion, die wörtlich ode45() heißt, aber SciPys solve_ivp() deckt denselben Arbeitsablauf ab: Sie beschreiben die Gleichung, geben dem Solver ein Zeitintervall und einen Anfangszustand vor und lesen die Lösung aus dem zurückgegebenen Objekt.

Verwenden Sie SciPys solve_ivp() für MATLABs ode45()

Die Funktion solve_ivp() in scipy.integrate ist die nächstliegende praktische Entsprechung zu ode45(). Ihre Standardmethode RK45 ist ein explizites Runge-Kutta-Verfahren der Ordnung 5(4) mit automatischer Schrittweitensteuerung – dieselbe Algorithmenfamilie, die ode45() in MATLAB verwendet. Der Solver wählt daher größere Schritte, wo sich die Lösung langsam ändert, und kleinere Schritte, wo sie sich schnell ändert.

Es ist kein Drop-in-Ersatz: Die Argumentnamen unterscheiden sich, und das Ergebnis wird als einzelnes Objekt zurückgegeben statt als Paar aus t- und y-Arrays. Das allgemeine Muster – eine Ableitungsfunktion, ein Zeitintervall und ein Anfangszustand – lässt sich jedoch direkt zwischen beiden übertragen. Installieren Sie SciPy mit pip install scipy, falls es noch nicht verfügbar ist.

Definieren Sie die ODE und rufen Sie solve_ivp() auf

Der wesentliche Aufruf ist solve_ivp(fun, t_span, y0, method="RK45", t_eval=None, args=None, ...). Die wichtigsten Eingaben sind:

• fun: ein aufrufbares Objekt fun(t, y), das die Ableitung dy/dt für den aktuellen Zeitpunkt t und den Zustand y zurückgibt. Bei einer einzelnen Gleichung gibt es einen Wert zurück; bei einem System ein Array mit einer Ableitung pro Zustandskomponente.
• t_span: ein Paar (t_start, t_end), das Anfang und Ende des Integrationsintervalls angibt.
• y0: der Zustand bei t_start; seine Länge legt die Anzahl der Gleichungen im System fest.
• method (optional): das Integrationsverfahren, standardmäßig 'RK45'.
• args (optional): zusätzliche Werte, die der Funktion als fun(t, y, *args) übergeben werden.
• rtol und atol (optional): relative und absolute Toleranzen, die die adaptive Schrittweite steuern.

Dieses Beispiel löst dy/dt = 2y ausgehend von y(0) = 1 über dem Intervall [0, 2]. Der Aufruf gibt ein Ergebnisobjekt zurück; das Ausgeben seiner Attribute t und y zeigt die vom Solver gewählten Zeitpunkte und den Zustand zu jedem dieser Zeitpunkte. Code und Ausgabe bleiben unverändert, einschließlich der englischen Beschriftungen t: und y::

Der Code und die angezeigten Bezeichnungen oder Ausgaben in diesem Abschnitt bleiben auf Englisch.

from scipy.integrate import solve_ivp


def growth(t, y):
    return 2 * y


result = solve_ivp(growth, (0, 2), [1])
print("t:", result.t)
print("y:", result.y[0])

t: [0.         0.07580098 0.56590702 1.18757739 1.85907272 2.        ]
y: [ 1.          1.16369695  3.10133569 10.75272756 41.18553629 54.59497095]

Die Abstände der ausgegebenen Zeitpunkte sind ungleichmäßig, weil solve_ivp() jeden Schritt an das lokale Verhalten der Lösung anpasst, statt mit einem festen Inkrement fortzuschreiten.

Lesen Sie das Solver-Ergebnis

solve_ivp() gibt ein scipy.integrate.OdeResult-Objekt zurück. Die Attribute, die Sie am häufigsten verwenden werden, sind:

• t: ein NumPy-Array der ausgegebenen Zeitpunkte. Ohne t_eval sind dies die internen adaptiven Schritte des Solvers.
• y: ein 2-D-NumPy-Array der Form (n, len(t)); Zeile i enthält Komponente i des Zustands zu jedem ausgegebenen Zeitpunkt, sodass y[:, k] der Zustand bei t[k] ist. Bei einer einzelnen Gleichung ist y[0] die Lösungskurve.
• sol: ein stetiges OdeSolution-Objekt, das Sie als result.sol(t) aufrufen können, um die Lösung zu beliebigen Zeitpunkten innerhalb des Intervalls auszuwerten. Es wird nur befüllt, wenn Sie dense_output=True übergeben; andernfalls ist es None.
• success und message: ein boolescher Wert dafür, ob die Integration abgeschlossen wurde, und eine menschenlesbare Beschreibung, wie sie endete.
• status: 0, wenn das Intervall vollständig durchlaufen wurde, 1, wenn ein Terminierungsereignis eingetreten ist, und -1, wenn der Solver keine Fortschritte gemacht hat.

Es ist eine gute Angewohnheit, success (oder status) zu prüfen, bevor Sie y verwenden, da ein fehlgeschlagener Lauf dennoch partielle Arrays zurückgeben kann.

Wählen Sie Ausgabezeitpunkte mit t_eval

Die adaptiven Schrittzeitpunkte stimmen selten mit den Zeitpunkten überein, die Sie ausgeben möchten. Übergeben Sie ein sortiertes Array von Zeitpunkten an t_eval; solve_ivp() gibt dann die Lösung genau an diesen Zeitpunkten zurück statt zu seinen internen Schritten. Code und Ausgabe bleiben unverändert, einschließlich der englischen Beschriftungen t: und y::

Der Code und die angezeigten Bezeichnungen oder Ausgaben in diesem Abschnitt bleiben auf Englisch.

from scipy.integrate import solve_ivp


def growth(t, y):
    return 2 * y


times = [0.5, 1.0, 1.5, 2.0]
result = solve_ivp(growth, (0, 2), [1], t_eval=times)
print("t:", result.t)
print("y:", result.y[0])

t: [0.5 1.  1.5 2. ]
y: [ 2.71823524  7.38683161 20.07910026 54.59497095]

Nun enthält t nur die angeforderten Zeitpunkte, und y enthält jeweils den interpolierten Zustand dazu. Jeder Wert in t_eval muss innerhalb von t_span liegen; ein Zeitpunkt außerhalb des Integrationsintervalls bricht den Aufruf mit einem Fehler ab, der besagt, dass die Werte innerhalb von t_span liegen müssen. Endpunkte sind erlaubt, sodass t_eval sowohl t_start als auch t_end enthalten darf.