Imitar la función ode45() en Python

Las ecuaciones diferenciales ordinarias se utilizan en MatLab para resolver muchos problemas científicos. El ode45() se usa en MatLab para resolver ecuaciones diferenciales.

Este artículo mostrará cómo podemos imitar la función ode45() en Python.

Imita la función ode45() en Python

Para imitar la función ode45() en python, podemos usar el método solve_ivp() definido en el módulo scipy. El método solve_ivp() integra un sistema de ecuaciones diferenciales ordinarias (EDO).

• El método solve_ivp() toma una función como primer argumento de entrada. La función dada en el argumento de entrada debe devolver un array que contenga los coeficientes de la ecuación diferencial.
• En el segundo argumento de entrada, el método solve_ivp() toma una tupla o lista que contiene dos valores numéricos. Los valores representan el intervalo de integración, donde el primer valor de la tupla representa el inicio del intervalo y el segundo valor de la tupla representa el valor más alto del intervalo.
• En el tercer argumento de entrada, el método solve_ivp() toma un array que representa los valores iniciales.
• Después de la ejecución, el método solve_ivp() devuelve un objeto montón con varios atributos.
  1. El atributo t contiene un array de numpy que contiene los puntos de tiempo.
  2. El atributo y contiene un array de numpy con valores y puntos de tiempo en t.
  3. El atributo sol contiene un objeto Odesolution que contiene la solución de la ecuación diferencial. Si el argumento dense_output se establece en falso en el método solve_ivp(), el atributo sol contiene None.

Para entender esto mejor, vea el siguiente ejemplo.

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)

Producción :

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

En el ejemplo anterior, primero definimos una función llamada function que toma t e y como argumento de entrada y devuelve un valor basado en y.

Luego, hemos definido un intervalo y valores iniciales para la EDO usando las variables intervalo y valores_inicial, respectivamente. Pasamos función, intervalo y valores_iniciales como argumentos de entrada a la función solve_ivp() y, por último, obtenemos el resultado en la solución variable.

En la salida, puede observar que los valores de tiempo se distribuyen a lo largo del intervalo de 0 a 10. De manera similar, la salida contiene un valor y correspondiente a cada valor de tiempo.

También podemos especificar explícitamente los puntos de tiempo en el atributo t de la solución. Para esto, necesitamos pasar un array que contenga los valores de tiempo deseados para los cuales necesitamos los valores y al argumento t_eval del método solve_ivp(), como se muestra a continuación.

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)

Producción :

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

Puede ver que los valores de tiempo solo contienen aquellos valores que se pasan como argumentos de entrada al parámetro t_eval. De manera similar, el atributo y contiene valores solo para los valores t especificados.

Este enfoque puede ayudarlo a obtener valores para ciertos puntos en el intervalo.