Recommandation pour la méthode des différences finies en Python scientifique

Pour un projet sur lequel je travaille (dans les PDE hyperboliques), j'aimerais avoir une idée approximative du comportement en regardant quelques chiffres. Je ne suis cependant pas un très bon programmeur.

Pouvez-vous recommander des ressources pour apprendre à coder efficacement des schémas de différences finies en Python scientifique (d'autres langages à petite courbe d'apprentissage sont également les bienvenus)?

Pour vous donner une idée du public (moi) pour cette recommandation:

• Je suis un pur mathématicien de formation et je connais un peu les aspects théoriques des schémas de différences finies
• Ce dont j'ai besoin d'aide, c'est de savoir comment faire pour que l'ordinateur calcule ce que je veux qu'il calcule, en particulier de manière à ne pas dupliquer trop les efforts déjà déployés par d'autres (afin de ne pas réinventer la roue lorsque un package est déjà disponible). (Une autre chose que j'aimerais éviter est de coder stupidement quelque chose à la main lorsqu'il existe des structures de données établies correspondant à l'objectif.)
• J'ai eu une certaine expérience de codage; mais je n'en ai pas eu en Python (donc cela ne me dérange pas s'il existe de bonnes ressources pour apprendre une langue différente [disons, Octave par exemple]).
• Livres, documentation seraient tous deux utiles, tout comme les collections d'exemples de code.

Voici un exemple de 97 lignes de résolution d'un PDE multivarié simple utilisant des méthodes de différences finies, contribué par le professeur David Ketcheson , du référentiel py4sci que je maintiens. Pour les problèmes plus compliqués où vous devez gérer des chocs ou la conservation dans une discrétisation à volume fini, je recommande de regarder pyclaw , un logiciel que j'aide à développer.

"""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()

Vous pouvez jeter un œil à Fenics , qui est un cadre en python / C qui permet de résoudre des équations assez générales à l'aide d'un langage de balisage spécial. Il utilise principalement des éléments finis, mais mérite le détour. Le didacticiel devrait vous donner une idée de la facilité avec laquelle il peut être possible de résoudre des problèmes.

Cette référence pourrait vous être très utile. Ceci est un livre ouvert sur Internet. J'ai appris (j'apprends toujours) le python de ce livre. Je l'ai trouvé en effet très bonne ressource.

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

Pour le calcul numérique, il faut absolument opter pour «numpy». (assurez-vous simplement que vous avez bien compris le 'tableau' et la 'matrice' et la 'liste') (reportez-vous à la documentation numpy pour cela)