Recomendación para el método de diferencias finitas en Python científico

Para un proyecto en el que estoy trabajando (en PDE hiperbólicas), me gustaría obtener una idea aproximada del comportamiento observando algunos números. Sin embargo, no soy un muy buen programador.

¿Puede recomendarme algunos recursos para aprender a codificar efectivamente esquemas de diferencias finitas en Scientific Python (también son bienvenidos otros idiomas con una curva de aprendizaje pequeña)?

Para darle una idea de la audiencia (yo) para esta recomendación:

• Soy un matemático puro por formación, y estoy algo familiarizado con los aspectos teóricos de los esquemas de diferencias finitas.
• Lo que necesito ayuda es cómo hacer que la computadora calcule lo que quiero que calcule, especialmente de una manera que no duplique demasiado el esfuerzo ya realizado por otros (para no reinventar la rueda cuando un paquete ya está disponible). (Otra cosa que me gustaría evitar es codificar estúpidamente algo a mano cuando hay estructuras de datos establecidas que se ajustan al propósito).
• He tenido algo de experiencia en codificación; pero no he tenido ninguno en Python (por lo tanto, no me importa si hay buenos recursos para aprender un idioma diferente [por ejemplo, Octave]).
• Los libros, la documentación serían útiles, al igual que las colecciones de código de ejemplo.

Aquí hay un ejemplo de 97 líneas para resolver una PDE multivariada simple usando métodos de diferencias finitas, aportados por el profesor David Ketcheson , del repositorio py4sci que mantengo. Para problemas más complicados en los que necesita manejar choques o conservación en una discretización de volumen finito, le recomiendo que busque Pyclaw , un paquete de software que ayudo a desarrollar.

"""Pattern formation code

    Solves the pair of PDEs:
       u_t = D_1 \nabla^2 u + f(u,v)
       v_t = D_2 \nabla^2 v + g(u,v)
"""

import matplotlib
matplotlib.use('TkAgg')
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import spdiags,linalg,eye
from time import sleep

#Parameter values
Du=0.500; Dv=1;
delta=0.0045; tau1=0.02; tau2=0.2; alpha=0.899; beta=-0.91; gamma=-alpha;
#delta=0.0045; tau1=0.02; tau2=0.2; alpha=1.9; beta=-0.91; gamma=-alpha;
#delta=0.0045; tau1=2.02; tau2=0.; alpha=2.0; beta=-0.91; gamma=-alpha;
#delta=0.0021; tau1=3.5; tau2=0; alpha=0.899; beta=-0.91; gamma=-alpha;
#delta=0.0045; tau1=0.02; tau2=0.2; alpha=1.9; beta=-0.85; gamma=-alpha;
#delta=0.0001; tau1=0.02; tau2=0.2; alpha=0.899; beta=-0.91; gamma=-alpha;
#delta=0.0005; tau1=2.02; tau2=0.; alpha=2.0; beta=-0.91; gamma=-alpha; nx=150;

#Define the reaction functions
def f(u,v):
    return alpha*u*(1-tau1*v**2) + v*(1-tau2*u);

def g(u,v):
    return beta*v*(1+alpha*tau1/beta*u*v) + u*(gamma+tau2*v);


def five_pt_laplacian(m,a,b):
    """Construct a matrix that applies the 5-point laplacian discretization"""
    e=np.ones(m**2)
    e2=([0]+[1]*(m-1))*m
    h=(b-a)/(m+1)
    A=np.diag(-4*e,0)+np.diag(e2[1:],-1)+np.diag(e2[1:],1)+np.diag(e[m:],m)+np.diag(e[m:],-m)
    A/=h**2
    return A

def five_pt_laplacian_sparse(m,a,b):
    """Construct a sparse matrix that applies the 5-point laplacian discretization"""
    e=np.ones(m**2)
    e2=([1]*(m-1)+[0])*m
    e3=([0]+[1]*(m-1))*m
    h=(b-a)/(m+1)
    A=spdiags([-4*e,e2,e3,e,e],[0,-1,1,-m,m],m**2,m**2)
    A/=h**2
    return A

# Set up the grid
a=-1.; b=1.
m=100; h=(b-a)/m; 
x = np.linspace(-1,1,m)
y = np.linspace(-1,1,m)
Y,X = np.meshgrid(y,x)

# Initial data
u=np.random.randn(m,m)/2.;
v=np.random.randn(m,m)/2.;
plt.hold(False)
plt.pcolormesh(x,y,u)
plt.colorbar; plt.axis('image'); 
plt.draw()
u=u.reshape(-1)
v=v.reshape(-1)

A=five_pt_laplacian_sparse(m,-1.,1.);
II=eye(m*m,m*m)

t=0.
dt=h/delta/5.;
plt.ion()

#Now step forward in time
for k in range(120):
    #Simple (1st-order) operator splitting:
    u = linalg.spsolve(II-dt*delta*Du*A,u)
    v = linalg.spsolve(II-dt*delta*Dv*A,v)

    unew=u+dt*f(u,v);
    v   =v+dt*g(u,v);
    u=unew;
    t=t+dt;

    #Plot every 3rd frame
    if k/3==float(k)/3:
        U=u.reshape((m,m))
        plt.pcolormesh(x,y,U)
        plt.colorbar
        plt.axis('image')
        plt.title(str(t))
        plt.draw()

plt.ioff()

Podría echar un vistazo a Fenics , que es un marco de Python / C que permite resolver ecuaciones bastante generales utilizando un lenguaje de marcado especial. Sin embargo, utiliza principalmente elementos finitos, pero vale la pena echarle un vistazo. El tutorial debería darle una impresión de lo fácil que puede ser resolver problemas.

Esta referencia puede ser muy útil para usted. Este es un libro abierto en internet. Aprendí (todavía estoy aprendiendo), Python de este libro. Me pareció muy buen recurso de hecho.

http://www.openbookproject.net/thinkcs/python/english2e/

Para el cálculo numérico, uno definitivamente debe ir por 'numpy'. (solo asegúrese de haber entendido la 'matriz' y 'matriz' y 'lista' correctamente) (consulte la documentación numpy para eso)