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

20

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.
Willie Wong
la source
Le problème principal est que je ne sais même pas par où commencer: donc même des suggestions de base seraient utiles.
Willie Wong
La restriction est seulement que je ne suis pas (encore) familier avec les méthodes de volumes finis; donc je vais devoir apprendre la méthode en même temps. Je ne m'opposerais pas à une telle réponse, bien sûr.
Willie Wong
PyClaw peut gérer des termes sources non linéaires, mais l'écriture de votre propre solveur Riemann serait compliquée, en particulier dans les dimensions 2 ou supérieures. Si vous voulez essayer un schéma de différenciation finie simple avec des grilles structurées, votre prochaine option serait d'essayer quelque chose dans petsc4py , (Divulgation: je suis également affilié à ce projet), qui est plus général et pas aussi bien- documenté.
Aron Ahmadia
Salut Willie (et pour les lecteurs qui n'ont pas regardé le chat), je pense que vous le savez déjà, mais puisque vous avez mentionné les PDE hyperboliques, vous seriez probablement mieux avec une méthode de volume fini.
Matthew Emmett

Réponses:

10

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

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.

moyner
la source
3

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)

Subodh
la source