Pybonacci

Computación Científica con Python en castellano

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

with 3 comments

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.

Región de estabilidad absoluta

A la hora de integrar numéricamente un problema diferencial del tipo

\frac{d U}{d t} = A U

nos interesa que la solución numérica tenga el mismo carácter de estabilidad que la solución analítica y saber para qué valores del paso de integración podemos conseguir esto [Hernández]. Mediante los autovalores \lambda de la matriz del sistema o los de la matriz jacobiana que resulta de linealizar el sistema [LeVeque], podemos definir la

Región de estabilidad absoluta: región del plano complejo \mathcal{R} \subset \mathbb{C} tal que el método es estable \forall \, \omega = \lambda \Delta t \in \mathcal{R}.

Si la solución de la ecuación en diferencias resultante de discretizar el problema diferencial es de la forma

u^n = {\displaystyle \sum_k c_k r_k^n},

entonces los r_k deben ser raíces del polinomio característico de estabilidad del método numérico:

\pi(r) |_{\omega} = {\displaystyle \sum_{j = 0}^p (\alpha_j - \omega \beta_j) r^{p - j}}

que no depende del problema que estamos integrando. Si se cumple que el mayor de los valores absolutos de las raíces es menor o igual que 1 entonces la solución numérica será estable (asintóticamente estable en el caso estrictamente menor), e inestable si hay alguna raíz mayor que 1 en valor absoluto. Con esta información, vamos a computar la región de estabilidad absoluta siguiendo el algoritmo que propone [Hernández]:

Algoritmo

  1. Discretizar una región del plano complejo (x_i, y_j).
  2. Sea \rho una matriz. Para cada \omega_{ij} \leftarrow x_i + i y_j,
    1. Obtener las raíces r_k \quad (k = 1, \, \dots, p) de \pi(r) |_{\omega_{ij}}.
    2. \rho_{ij} \leftarrow \max(|r_k|)
  3. Representar las curvas de nivel de \rho. La \rho = 1 será la frontera de la región de estabilidad absoluta y las \rho < 1 serán el interior de la misma.

Código y representación

Vamos a incluir el código en un módulo para poder acceder a esta funcionalidad más fácilmente. El código completo, añadiendo documentación y ejemplos con el estándar de documentación de NumPy/SciPy, quedaría así:

# -*- coding: utf-8 -*-
#
# Región de estabilidad absoluta
# Juan Luis Cano Rodríguez

import numpy as np

def region_estabilidad(p, X, Y):
    """Región de estabilidad absoluta

    Computa la región de estabilidad absoluta de un método numérico, dados
    los coeficientes de su polinomio característico de estabilidad.

    Argumentos
    ----------
    p : function
        Acepta un argumento w y devuelve una lista de coeficientes
    X, Y : numpy.ndarray
        Rejilla en la que evaluar el polinomio de estabilidad generada por
        numpy.meshgrid

    Devuelve
    --------
    Z : numpy.ndarray
        Para cada punto de la malla, máximo de los valores absolutos de las
        raíces del polinomio de estabilidad

    Ejemplos
    --------
    >>> import numpy as np
    >>> x = np.linspace(-3.0, 1.5)
    >>> y = np.linspace(-3.0, 3.0)
    >>> X, Y = np.meshgrid(x, y)
    >>> Z = region_estabilidad(lambda w: [1,
    ... -1 - w - w ** 2 / 2 - w ** 3 / 6 - w ** 4 / 24], X, Y)  # RK4
    >>> import matplotlib.pyplot as plt
    >>> cs = plt.contour(X, Y, Z, np.linspace(0.05, 1.0, 9))
    >>> plt.clabel(cs, inline=1, fontsize=10)  # Para etiquetar los contornos
    >>> plt.show()

    """

    Z = np.zeros_like(X)
    w = X + Y * 1j
    
    for j in range(len(X)):
        for i in range(len(Y)):
            r = np.roots(p(w[i, j]))
            Z[i, j] = np.max(abs(r if np.any(r) else 0))

    return Z

Como sucede casi siempre en Python, el código es casi una traducción literal del algoritmo. Vamos a explicar un poco el código:

  1. Obtenemos los coeficientes del polinomio para w[i, j].
  2. Calculamos las raíces de dicho polinomio con np.roots, que acepta un array de rango 1 o equivalente (por esto no se puede vectorizar el bucle) y devuelve un array con las raíces del polinomio.
  3. La expresión r if np.any(r) else 0 es un condicional ternario, y tiene la misión de devolver un 0 si el polinomio no tiene raíces (porque todos los coeficientes se han hecho nulos, por ejemplo) para que la función np.max no reciba un array vacío. Es equivalente a
    if np.any(r):
        return r
    else:
        return 0
    
  4. Calculamos el valor absoluto de estas raíces.
  5. Con la función np.max calculamos el mayor de estos valores.
  6. Lo asignamos a Z[i, j].

Resultados

Una vez tenemos esto en un archivo region_estabilidad.py, si abrimos un intérprete de IPython en el mismo directorio será tan sencillo como

In [1]: import numpy as np

In [2]: x = np.linspace(-3.0, 1.5)

In [3]: y = np.linspace(-3.0, 3.0)

In [4]: X, Y = np.meshgrid(x, y)

In [5]: def p(w):
   ...:     """Polinomio de estabilidad del método RK4."""
   ...:     return [1, -1 - w - w ** 2 / 2 - w ** 3 / 6 - w ** 4 / 24]
   ...:

In [6]: from region_estabilidad import region_estabilidad

In [7]: Z = region_estabilidad(p, X, Y)

In [8]: import matplotlib.pyplot as plt

In [9]: plt.contour(X, Y, Z, np.linspace(0.0, 1.0, 9))
Out[9]:

In [10]: plt.show()

El segundo argumento de plt.contour es para especificar los niveles que queremos pintar (normalmente sólo nos interesará hasta el 1.0). Podemos obtener resultados similares a estos:

Región de estabilidad absoluta de un método Adams-Bashforth 3

Región de estabilidad absoluta de un método Runge-Kutta 4

Y hasta aquí llega el artículo de hoy. Espero que os haya resultado interesante y útil, no olvidéis difundir y comentar :)

¡Un saludo!

Referencias

  • HERNÁNDEZ, Juan A. Cálculo Numérico en Ecuaciones Diferenciales Ordinarias. ADI, 2001.
  • LEVEQUE, Randall J. Finite Difference Methods for Ordinary and Partial Differential Equations. Society for Industrial and Applied Mathematics, 2007.
About these ads

Written by Juanlu001

27 de abril de 2012 a 9:36

3 comentarios

Suscríbete a los comentarios mediante RSS.

  1. Yo también estoy en 3º de aeronáutica y estaba comprobando mis resultados con los tuyos, y creo que debe de haber algún error en la región de estabilidad del AB2, pues a mi me sale una región que está en el plano Re(omega)<0. El único punto del eje imaginario que entra dentro de la región es el (0,0). Lo miré en el libro de Hernández y le sale igual así que te lo aviso para que no lo entregues así. Gracias por el artículo.

    Manuel

    15 de mayo de 2012 at 2:50

    • Gracias Manuel, pero en el artículo está pintada la región del AB3! :P

      Un saludo

      Juanlu001

      15 de mayo de 2012 at 8:19

      • Mis disculpas.

        Manuel

        15 de mayo de 2012 at 15:30


¡Deja un comentario!

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s

Seguir

Recibe cada nueva publicación en tu buzón de correo electrónico.

Únete a otros 1.103 seguidores

%d personas les gusta esto: