Integrando ecuaciones diferenciales: método leapfrog en Python

Introducción

En Python tenemos numerosas herramientas listas para que podamos integrar ecuaciones diferenciales ordinarias sin tener que preocuparnos en implementar un esquema numérico. Sin ir más lejos, en el módulo integrate de SciPy existen varias funciones a tal efecto. La función odeint es una interfraz en Python a la biblioteca ODEPACK, escrita en Fortran. Sin embargo, con un propósito didáctico vamos a estudiar cómo programaríamos la solución a un problema de ecuaciones diferenciales, en este caso utilizando la regla del punto medio o método leapfrog[1].

Para esta entrada se ha utilizado python 2.7.2, numpy 1.6.1 y matplotlib 1.1.0.

Enunciado

Problema de Cauchy: integrar la siguiente ecuación diferencial
\( \frac{d^2 x}{d t^2} + x = 0\)
Con las condiciones iniciales \( x(0) = 1, \dot{x}(0) = 0\).

Tratamiento matemático

Sabemos que la solución de este sistema es:

\( x(t) = \cos{t}\)

Siento fastidiar la sorpresa 🙂

En primer lugar, hemos de darle a este problema el tratamiento matemático necesario para poder resolverlo numéricamente con facilidad. Para ello, transformaremos la EDO (Ecuación Diferencial Ordinaria) de 2º orden en un sistema de dos EDOs de 1º orden. Definiendo el vector

\( U = \begin{pmatrix} x \\ \dot{x} \end{pmatrix}\)

tendremos

\( \frac{d U}{d t} = \begin{pmatrix} \dot{x} \\ -x \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} x \\ \dot{x} \end{pmatrix} = A U = F(U)\)

donde \( A\) es la matriz del sistema. Utilizaremos el esquema leapfrog para resolver el sistema. Su EDD (Ecuación en Diferencias) es

\( U^{n + 1} = U^{n – 1} + 2 \Delta t F(U^n)\)

Como se puede ver, no podemos aplicar directamente el esquema para el paso 1 porque no tenemos \( U^{-1}\). Debemos arrancar con un esquema de un paso, como por ejemplo un Euler explícito:

\( U^{n + 1} = U^n + \Delta t F(U^n)\)

Puesto que tenemos

\( U^0 = \begin{pmatrix} x(0) \\ \dot{x}(0) \end{pmatrix} = \begin{pmatrix} 1.0 \\ 0.0 \end{pmatrix}\)

ya podemos resolver el problema, sin más que escoger un valor del paso.

Implementación

Aunque no sería estrictamente necesario, como la expresión de la matriz del sistema es muy sencilla vamos a utilizar las capacidades de NumPy para manejar matrices. La función \( F(U)\) quedará, en Python,

Aquí hemos utilizado la función np.dot(a, b), que para el caso de que a sea un array bidimensional y b un array unidimensional es el producto de matriz por vector al que estamos acostumbrados en Álgebra Lineal.

Nota: Nótese que no podíamos escribir directamente A * u. ¡Esto no sería el producto de matriz por vector, sería el producto de dos arrays!

Ahora llegó el momento de implementar los esquemas Euler y leapfrog. Como queremos que sea lo más general posible, aceptaremos que la función F sea cualquiera y pueda depender del tiempo, que el paso dt se pueda también escoger y que cada vez que llamemos a las funciones nos den solamente el paso siguiente.

Para el caso del Euler, los datos de entrada serán el instante t_n, el vector u_n, la función F(t, u) y el paso dt, y el método nos dará el vector u_n1.

Traducción literal de la EDD del esquema a Python. Más sencillo imposible 🙂
Hacemos lo mismo con el esquema leapfrog:

Ya sólo nos queda implementar la lógica del programa.

Solución numérica

En primer lugar decidiremos el número n de pasos que queremos dar, o hasta dónde queremos hallar la solución, y guardaremos los sucesivos valores de x en un array de dimensión n. Para ello utilizamos la función empty(shape), que nos inicializa un array con la forma dada por el primer argumento, y linspace(a, b, n), que discretizará el intervalo \( [a, b]\) con n puntos. El código correspondiente a esto y a dar las condiciones iniciales quedará:

Ya podemos empezar a integrar. El primer paso lo daremos con el euler, y los que queden hasta n con el método leapfrog. Después de cada paso guardamos el valor x hallado, y vamos avanzando. Como necesitaremos guardar dos pasos del vector U para poder aplicar el método leapfrog, tendremos que escribir a continuación

A partir de ahora, todos los pasos son iguales. Iremos sobreescribiendo en u_n0 y u_n1 los valores de los vectores que necesitemos para cada paso del leapfrog, y en u_n2 escribiremos la solución dada por el método. El código del bucle será, finalmente,

¡Fácil, rápido y para toda la familia! 🙂

Representación gráfica

Y ya, para dar el toque de gracia, representemos gráficamente la solución con estas sencillas líneas:

¿Esperabas que fuese más difícil? 😀

Conclusión

El código final es este:

Como se puede ver, en 60 líneas de código incluyendo abundantes comentarios y espacios en blanco hemos implementado un esquema numérico para resolver un problema de ecuaciones diferenciales y hemos representado la solución, sin salirnos de Python.

Espero que te haya resultado útil, nos vemos en Pybonacci 🙂

Referencias

2 pensamientos sobre “Integrando ecuaciones diferenciales: método leapfrog en Python”

  1. El leapfrog no es estable para ecuaciones disipativas como la de tu ejemplo. Si integras a tiempo mayor verás que deja de converger al cabo de un rato. Lo ves enseguida si te pintas la región de convergencia. El leapfrog funciona sólo con ecuaciones puramente conservativas como las que te encuentras en mecánica molecular.

    1. Tu comentario me produce unas cuantas dudas. Pensaba coger mañana un par de referencias sobre análisis numérico para solucionarlas, pero así para empezar lo que más me llama la atención es que digas que esta ecuación es disipativa: lo que me dice el sentido común es que lo sería si tuviese un término de fricción. ¿En qué me estoy equivocando?

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

fifteen − = 12