Imitieren die Funktion ode45() in Python

Gewöhnliche Differentialgleichungen werden in MatLab verwendet, um viele wissenschaftliche Probleme zu lösen. Die ode45() wird in MatLab verwendet, um Differentialgleichungen zu lösen.

Dieser Artikel zeigt, wie wir die Funktion ode45() in Python nachahmen können.

Imitieren Sie die Funktion ode45() in Python

Um die Funktion ode45() in Python zu imitieren, können wir die im Modul scipy definierte Methode solve_ivp() verwenden. Die Methode solve_ivp() integriert ein System gewöhnlicher Differentialgleichungen (ODEs).

• Die Methode solve_ivp() nimmt eine Funktion als erstes Eingabeargument. Die im Eingabeargument angegebene Funktion muss ein Array zurückgeben, das die Koeffizienten der Differentialgleichung enthält.
• Im zweiten Eingabeargument nimmt die Methode solve_ivp() ein Tupel oder eine Liste mit zwei numerischen Werten. Die Werte stellen das Integrationsintervall dar, wobei der erste Wert im Tupel den Beginn des Intervalls und der zweite Wert des Tupels den höchsten Wert im Intervall darstellt.
• Im dritten Eingabeargument nimmt die Methode solve_ivp() ein Array, das die Anfangswerte darstellt.
• Nach der Ausführung liefert die Methode solve_ivp() ein Bunch-Objekt mit verschiedenen Attributen zurück.
  1. Das Attribut t enthält ein numpy-Array mit Zeitpunkten.
  2. Das Attribut y enthält ein numpy-Array mit Werten und Zeitpunkten in t.
  3. Das Attribut sol enthält ein Objekt Odesolution, das die Lösung der Differentialgleichung enthält. Wenn in der Methode solve_ivp() das Argument dense_output auf false gesetzt ist, enthält das Attribut sol None.

Um dies besser zu verstehen, sehen Sie sich das folgende Beispiel an.

from scipy.integrate import solve_ivp


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


interval = [0, 10]
initial_values = [10, 15, 25]
solution = solve_ivp(function, interval, initial_values)
print("Time:", solution.t)
print("Y:", solution.y)

Ausgabe:

Time: [ 0.          0.07578687  0.56581063  1.18741382  1.85887096  2.55035821
  3.25007544  3.95320486  4.65775424  5.36289544  6.06828346  6.77377445
  7.47930839  8.18486026  8.89041961  9.59598208 10.        ]
Y: [[1.00000000e+01 1.16366412e+01 3.10073783e+01 1.07492109e+02
  4.11689241e+02 1.64114780e+03 6.65071446e+03 2.71362627e+04
  1.11036049e+05 4.54874443e+05 1.86437495e+06 7.64300835e+06
  3.13352156e+07 1.28474398e+08 5.26752964e+08 2.15973314e+09
  4.84541488e+09]
 [1.50000000e+01 1.74549617e+01 4.65110674e+01 1.61238163e+02
  6.17533861e+02 2.46172171e+03 9.97607169e+03 4.07043941e+04
  1.66554074e+05 6.82311665e+05 2.79656243e+06 1.14645125e+07
  4.70028233e+07 1.92711598e+08 7.90129446e+08 3.23959970e+09
  7.26812231e+09]
 [2.50000000e+01 2.90916029e+01 7.75184457e+01 2.68730272e+02
  1.02922310e+03 4.10286951e+03 1.66267862e+04 6.78406569e+04
  2.77590123e+05 1.13718611e+06 4.66093739e+06 1.91075209e+07
  7.83380389e+07 3.21185996e+08 1.31688241e+09 5.39933284e+09
  1.21135372e+10]]

Im obigen Beispiel haben wir zuerst eine Funktion namens Funktion definiert, die t und y als Eingabeargument verwendet und einen Wert basierend auf y zurückgibt.

Dann haben wir ein Intervall und Anfangswerte für die ODE mit den Variablen interval bzw. initial_values definiert. Wir übergeben function, interval und initial_values als Eingabeargumente an die solve_ivp()-Funktion, und zuletzt erhalten wir die Ausgabe in der Variablen Lösung.

In der Ausgabe können Sie beobachten, dass die Zeitwerte über das Intervall 0 bis 10 verteilt sind. Ebenso enthält die Ausgabe einen y-Wert, der jedem Zeitwert entspricht.

Wir können die Zeitpunkte auch explizit im Attribut t der Lösung angeben. Dazu müssen wir ein Array mit den gewünschten Zeitwerten, für die wir die y-Werte benötigen, an das t_eval-Argument der solve_ivp()-Methode übergeben, wie unten gezeigt.

from scipy.integrate import solve_ivp


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


interval = [0, 10]
initial_values = [10, 15, 25]
time_values = [1, 2, 3, 6, 7, 8]
solution = solve_ivp(function, interval, initial_values, t_eval=time_values)
print("Time:", solution.t)
print("Y:", solution.y)

Ausgabe:

Time: [1 2 3 6 7 8]
Y: [[7.38683416e+01 5.46053271e+02 4.03089733e+03 1.62618365e+06
  1.20160156e+07 8.87210156e+07]
 [1.10802512e+02 8.19079906e+02 6.04634600e+03 2.43927547e+06
  1.80240234e+07 1.33081523e+08]
 [1.84670854e+02 1.36513318e+03 1.00772433e+04 4.06545912e+06
  3.00400390e+07 2.21802539e+08]]

Sie sehen, dass die Zeitwerte nur die Werte enthalten, die als Eingabeargumente an den Parameter t_eval übergeben werden. Ebenso enthält das Attribut y nur Werte für die angegebenen t-Werte.

Dieser Ansatz kann Ihnen helfen, Werte für bestimmte Punkte im Intervall zu erhalten.