FEniCS: Resolución de ecuaciones diferenciales en Python

Introducción

En este artículo os voy a presentar el proyecto FEniCS, una colección de bibliotecas escritas con interfaz en Python para la resolución de ecuaciones diferenciales por el método de los elementos finitos. FEniCS proporciona un método muy inteligente para automatizar los procesos más laboriosos de la solución de este tipo de ecuaciones, de forma que podemos atacar cualquier problema que se nos ocurra. Las posibilidades son inmensas y la documentación es bastante exhaustiva, así que aquí veremos una introducción a cómo resolver ecuaciones en derivadas parciales con FEniCS.

Placa rectangular simplemente apoyada resuelta con FEniCS
Placa rectangular simplemente apoyada resuelta con FEniCS

Continue reading

Ecuaciones de Lotka-Volterra: modelo presa depredador

Introducción

Resulta intuitivo pensar que las poblaciones de un animal depredador y su presa están relacionadas de algún modo en el que si una aumenta, la otra lo hace también. Utilizaremos como ejemplo en este artículo un ecosistema aislado y formado por leones y cebras que viven en armonía, es decir, los leones comiéndose a las cebras. Imaginemos que por cualquier circunstancia, por ejemplo, por disponer de mayor cantidad de alimento, aumenta la población de cebras; los leones dispondrán de más alimento y su población aumentará, pero ¿qué ocurrirá a partir de este momento? Si la población de leones llega a ser demasiado grande para el número de cebras en nuestra sabana, podrían acabar con todas, provocando su propia extinción por inanición. Pero incluso si el festín no es tan grande como para comerse todas las cebras, pero sí para dejar una población muy mermada, probablemente los leones tendrán que pasar hambre una buena temporada y algunos de ellos morirán hasta que las cebras tengan tiempo suficiente para reproducirse y volver a ser pasto de los leones. ¿Cuántas cebras morirán en el atracón? ¿Cuánto tiempo pasarán los leones hambre? ¿Cuántos morirán?

Continue reading

Cómo crear una matriz tridiagonal en Python con NumPy y SciPy

Introducción

En este rápido apunte vamos a ver cómo construir una matriz tridiagonal en Python utilizando NumPy y SciPy. Una matriz tridiagonal es una matriz cuadrada que solamente tiene elementos distintos de cero en su diagonal principal y en las dos diagonales adyacentes a esta (la superdiagonal y la subdiagonal). Las matrices tridiagonales aparecen mucho en cálculo numérico, por ejemplo en la discretización de ecuaciones diferenciales, y tienen la característica de ser matrices dispersas (en lugar de densas) al ser la mayoría de sus elementos cero.

Sin que sirva de precedente, hoy vamos a escribir código que sea compatible tanto con Python 2 como con Python 3. Es un cambio nimio, pero merece la pena ir acostumbrándose a pensar que tarde o temprano habrá que abandonar Python 2 🙂

En esta entrada se ha usado python 2.7.3, numpy 1.6.1 y scipy 0.10.1 y es compatible con python 3.2.3.

Continue reading

Regiones de estabilidad de métodos numéricos en Python

Introducción

Hoy veremos cómo computar con Python la región de estabilidad absoluta de un método numérico para resolver problemas de valores iniciales en ecuaciones diferenciales ordinarias, una herramienta muy importante a la hora de escoger un método numérico adecuado para integrar nuestro problema concreto. Se trata simplemente de otro ejemplo aplicado de lo que publicamos hace unos días sobre cómo pintar curvas de nivel con matplotlib; si quieres ver otro más, puedes leer nuestro ejemplo de uso de Basemap y NetCDF4, donde vimos cómo representar datos climatológicos.

En esta entrada se ha usado python 2.7.3, numpy 1.6.1 y matplotlib 1.1.0.

Continue reading

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\).

Continue reading