10.6: Integrando ODEs con Scipy

Excepto con fines educativos, casi siempre es una mala idea implementar tu propio solucionador de ODE; en cambio, debes usar un solucionador preescrito.

10.6.1 El solucionador scipy.integrate.odeint

En Scipy, el solucionador de ODE más simple de usar es la función scipy.integrate.odeint, que se encuentra en el módulo scipy.integrate. Esto es en realidad una envoltura alrededor de una biblioteca numérica de bajo nivel conocida como LSODE (el L ivermore S olver para ODE s”), que es parte de una biblioteca de solucionador ODE ampliamente utilizada conocida como ODEPACK. La característica más importante de este solucionador es que es “adaptativo”: puede averiguar automáticamente (i) qué esquema de integración usar (eligiendo entre un método Adams-Moulton de alto orden, u otro método implícito conocido como la Fórmula de Diferenciación Hacia Atrás que no hemos descrito), y ii) el tamaño de los pasos de tiempo discretos, con base en el comportamiento de las soluciones a medida que se están elaborando. En otras palabras, el usuario solo necesita especificar la función derivada, el estado inicial y los tiempos de salida deseados, sin tener que preocuparse por los detalles internos del método de solución.

La función toma varias entradas, de las cuales las más importantes son:

1. func, una función correspondiente a la función derivada\(\vec{F}(\vec{y}, t)\).
2. y0, ya sea un array numérico o 1D, correspondiente al estado inicial\(\vec{y}(t_0)\).
3. t, una matriz de tiempos en los que se genera la solución ODE. El primer elemento correspondiente al tiempo inicial\(t_{0}\). Tenga en cuenta que estos son los tiempos de “salida” solamente, no especifican los pasos de tiempo reales que el solucionador utiliza para encontrar las soluciones; esos son determinados automáticamente por el solucionador.
4. (opcional) args, una tupla de entradas adicionales para pasar a la función derivada func. Por ejemplo, si args =( 2,3), entonces func debería aceptar cuatro entradas, y se le pasará 2 y 3 como las dos últimas entradas.

La función luego devuelve una matriz y, donde y [n] contiene la solución en el tiempo t [n]. Tenga en cuenta que y [0] será exactamente el mismo que la entrada y0, el estado inicial que especificó.

Aquí hay un ejemplo del uso de odeint para resolver el problema del oscilador armónico amortiguado\(m \ddot{x} = - \lambda \dot{x} - k x(t)\), usando el truco de vectorización mencionado anteriormente para convertirlo en una ODE de primer orden:

from scipy import *
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def ydot(y, t, m, lambd, k):
    x, v = y[0], y[1]
    return array([v, -(lambd/m) * v - k * x / m])

m, lambd, k = 1.0, 0.1, 1.0   # Oscillator parameters
y0 = array([1.0, 5.0])        # Initial conditions [x, v]
t  = linspace(0.0, 50.0, 100) # Output times

y = odeint(ydot, y0, t, args=(m, lambd, k))

## Plot x versus t
plt.plot(t, y[:,0], 'b-')
plt.xlabel('t')
plt.ylabel('x')
plt.show()

Hay una limitación importante de odeint: no maneja ODE complejas, y siempre asume que\(\vec{y}\) y\(\vec{F}\) son reales. Sin embargo, esto no es un problema en la práctica, porque siempre se puede convertir una ODE compleja de primer orden en una real, reemplazando los vectores complejos\(\vec{y}\) y\(\vec{F}\) por vectores reales de doble longitud:

\[\vec{y}' \equiv \begin{bmatrix}\mathrm{Re}(\vec{y})\\ \mathrm{Im}(\vec{y})\end{bmatrix}, \;\; \vec{F}' \equiv \begin{bmatrix}\mathrm{Re}(\vec{F})\\ \mathrm{Im}(\vec{F})\end{bmatrix}.\]

10.6.2 Los solucionadores scipy.integrate.ode

Aparte de odeint, Scipy proporciona una interfaz más general a una variedad de solucionadores de ODE, en la forma de la clase scipy.integrate.ode. Esta es una interfaz mucho más bajo; en lugar de llamar a una sola función, hay que crear un “objeto” ODE, luego usar los métodos de este objeto para especificar el tipo de solucionador de ODE a usar, las condiciones iniciales, etc.; luego hay que llamar repetidamente al método integrate del objeto ODE, para integrar la solución hasta cada paso de tiempo de salida deseado.

El es una inconsistencia extremadamente agravante entre la función odeint y esta clase ode: ¡se invierte el orden esperado de entradas para las funciones derivadas! La función odeint asume que la función derivada tiene la forma F (y, t), pero la clase ode asume que tiene la forma F (t, y). ¡Cuidado con esto!

Aquí hay un ejemplo del uso de la clase ode con el problema del oscilador armónico amortiguado\(m \ddot{x} = - \lambda \dot{x} - k x(t)\), usando un solucionador Runge-Kutta:

from scipy import *
import matplotlib.pyplot as plt
from scipy.integrate import ode

## Note the order of inputs (different from odeint)!
def ydot(t, y, m, lambd, k):
    x, v = y[0], y[1]
    return array([v, -(lambd/m) * v - k * x / m])

m, lambd, k = 1.0, 0.1, 1.0   # Oscillator parameters
y0 = array([1.0, 5.0])        # Initial conditions [x, v]
t  = linspace(0.0, 50.0, 100) # Output times

## Set up the ODE object
r = ode(ydot)
r.set_integrator('dopri5')    # A Runge-Kutta solver
r.set_initial_value(y0)
r.set_f_params(m, lambd, k)

## Perform the integration.  Note that the "integrate" method only integrates
## up to one single final time point, rather than an array of times.
x    = zeros(len(t))
x[0] = y0[0]
for n in range(1,len(t)):
    r.integrate(t[n])
    assert r.successful()
    x[n] = (r.y)[0]

## Plot x versus t
plt.plot(t, x, 'b-')
plt.xlabel('t')
plt.ylabel('x')
plt.show()

Consulte la documentación para obtener una lista más detallada de opciones, incluida la lista de solucionadores de ODE entre los que puede elegir.