Imitez la fonction ode45() en Python

Les équations différentielles ordinaires sont utilisées dans MatLab pour résoudre de nombreux problèmes scientifiques. Le ode45() est utilisé dans MatLab pour résoudre des équations différentielles.

Cet article va montrer comment on peut imiter la fonction ode45() en Python.

Imitez la fonction ode45() en Python

Pour imiter la fonction ode45() en python, on peut utiliser la méthode solve_ivp() définie dans le module scipy. La méthode solve_ivp() intègre un système d’équations différentielles ordinaires (ODEs).

• La méthode solve_ivp() prend une fonction comme premier argument d’entrée. La fonction donnée en argument d’entrée doit renvoyer un tableau contenant les coefficients de l’équation différentielle.
• Dans le second argument d’entrée, la méthode solve_ivp() prend un tuple ou une liste contenant deux valeurs numériques. Les valeurs représentent l’intervalle d’intégration, où la première valeur du tuple représente le début de l’intervalle et la deuxième valeur du tuple représente la valeur la plus élevée de l’intervalle.
• Dans le troisième argument d’entrée, la méthode solve_ivp() prend un tableau représentant les valeurs initiales.
• Après exécution, la méthode solve_ivp() renvoie un objet bouquet avec divers attributs.
  1. L’attribut t contient un tableau numpy contenant des points de temps.
  2. L’attribut y contient un tableau numpy avec des valeurs et des points de temps dans t.
  3. L’attribut sol contient un objet Odesolution contenant la solution de l’équation différentielle. Si l’argument dense_output est défini sur false dans la méthode solve_ivp(), l’attribut sol contient None.

Pour mieux comprendre cela, voir l’exemple suivant.

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)

Production:

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]]

Dans l’exemple ci-dessus, nous avons d’abord défini une fonction nommée fonction qui prend t et y comme argument d’entrée et renvoie une valeur basée sur y.

Ensuite, nous avons défini un intervalle et des valeurs initiales pour l’ODE en utilisant respectivement les variables interval et initial_values. Nous passons function, interval et initial_values comme arguments d’entrée à la fonction solve_ivp(), et enfin, nous obtenons la sortie dans la variable solution.

Dans la sortie, vous pouvez observer que les valeurs de temps sont réparties sur l’intervalle de 0 à 10. De même, la sortie contient une valeur y correspondant à chaque valeur de temps.

On peut aussi spécifier explicitement les points temporels dans l’attribut t de la solution. Pour cela, nous devons passer un tableau contenant les valeurs de temps souhaitées pour lesquelles nous avons besoin des valeurs y à l’argument t_eval de la méthode solve_ivp(), comme indiqué ci-dessous.

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)

Production:

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]]

Vous pouvez voir que les valeurs de temps ne contiennent que les valeurs qui sont passées comme arguments d’entrée au paramètre t_eval. De même, l’attribut y contient des valeurs uniquement pour les valeurs t spécifiées.

Cette approche peut vous aider à obtenir des valeurs pour certains points de l’intervalle.