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

20

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.
Willie Wong
fuente
El principal problema es que ni siquiera sé por dónde empezar a buscar, por lo que incluso las sugerencias básicas serían útiles.
Willie Wong
La restricción es solo que no estoy (todavía) familiarizado con los métodos de volumen finito; así que tendré que aprender el método en conjunto. No me opondría a tal respuesta, por supuesto.
Willie Wong
PyClaw puede manejar términos de fuente no lineales, pero escribir su propio solucionador de Riemann sería complicado, particularmente en 2da o más dimensiones. Si desea probar un esquema simple de diferencia finita con cuadrículas estructuradas, su próxima opción sería probar algo en petsc4py , (Divulgación: también estoy afiliado a este proyecto), que tiene un propósito más general y no tan bueno. documentado
Aron Ahmadia
Hola Willie (y para los lectores que no han mirado el chat), creo que ya lo sabes, pero dado que mencionaste PDE hiperbólicos, probablemente estarías mejor con un método de volumen finito.
Matthew Emmett

Respuestas:

10

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()
Aron Ahmadia
fuente
8

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.

moyner
fuente
3

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)

Subodh
fuente