Empfehlung für die Finite-Differenz-Methode in Scientific Python

20

Für ein Projekt, an dem ich arbeite (in hyperbolischen PDEs), möchte ich anhand einiger Zahlen einen groben Überblick über das Verhalten erhalten. Ich bin jedoch kein sehr guter Programmierer.

Können Sie einige Ressourcen für das Lernen , wie man effektiv empfehlen Code Finite - Differenzen - Schemata in Scientific Python (andere Sprachen mit kleiner Lernkurve auch willkommen)?

Um Ihnen eine Vorstellung vom Publikum (mir) für diese Empfehlung zu geben:

  • Ich bin ausgebildeter Mathematiker und mit den theoretischen Aspekten von Finite-Differenzen-Schemata einigermaßen vertraut
  • Ich brauche Hilfe dabei, wie ich den Computer dazu bringen kann, das zu berechnen, was ich möchte, insbesondere auf eine Art und Weise, bei der ich nicht zu viel von den Anstrengungen anderer dupliziere (um das Rad nicht neu zu erfinden, wenn ein Paket ist bereits verfügbar). (Eine andere Sache, die ich vermeiden möchte, ist, etwas dumm von Hand zu codieren, wenn es etablierte Datenstrukturen gibt, die zum Zweck passen.)
  • Ich habe einige Programmiererfahrungen gesammelt. Aber ich hatte noch keine in Python (daher macht es mir nichts aus, wenn es gute Ressourcen gibt, um eine andere Sprache zu lernen [zum Beispiel Octave]).
  • Bücher, Dokumentation wäre sowohl nützlich, als auch Sammlungen von Beispielcode.
Willie Wong
quelle
Das Hauptproblem ist, dass ich nicht einmal weiß, wo ich anfangen soll, also wären selbst grundlegende Vorschläge hilfreich.
Willie Wong
Die Einschränkung ist nur, dass ich (noch) nicht mit Methoden für endliche Volumina vertraut bin; also muss ich die methode in verbindung lernen. Eine solche Antwort würde ich natürlich nicht beanstanden.
Willie Wong
PyClaw kann mit nichtlinearen Quellbegriffen umgehen, aber das Schreiben eines eigenen Riemann-Lösers wäre insbesondere in der 2. oder höheren Dimension kompliziert. Wenn Sie ein einfaches Finite-Differenzen-Schema mit strukturierten Gittern ausprobieren möchten , ist Ihre nächste Option, etwas in petsc4py (Offenlegung: Ich bin auch mit diesem Projekt verbunden) auszuprobieren , das allgemeiner und nicht so gut ist. dokumentiert.
Aron Ahmadia
Lasst uns diese Diskussion im Chat fortsetzen
Aron Ahmadia
Hallo Willie (und für Leser, die sich den Chat nicht angesehen haben), ich glaube, Sie wissen das bereits, aber da Sie hyperbolische PDEs erwähnt haben, wären Sie mit einer Methode mit endlichem Volumen wahrscheinlich besser dran.
Matthew Emmett

Antworten:

10

Hier ist ein 97-Zeilen-Beispiel für die Lösung einer einfachen multivariaten PDE mit Finite-Differenzen-Methoden, die von Prof. David Ketcheson aus dem von mir verwalteten py4sci-Repository beigesteuert wurden . Bei komplizierteren Problemen, bei denen Sie mit Erschütterungen oder der Konservierung in einer Diskretisierung mit endlichem Volumen umgehen müssen , empfehle ich pyclaw , ein Softwarepaket, das ich entwickle.

"""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
quelle
8

Sie können sich Fenics ansehen , ein Python / C-Framework, mit dem sich ganz allgemeine Gleichungen mit einer speziellen Auszeichnungssprache lösen lassen. Meistens werden Finite Elemente verwendet, die jedoch einen Blick wert sind. Das Tutorial soll Ihnen einen Eindruck vermitteln, wie einfach es sein kann, Probleme zu lösen.

Moyner
quelle
3

Diese Referenz könnte für Sie sehr nützlich sein. Dies ist ein offenes Buch im Internet. Ich habe Python aus diesem Buch gelernt (lerne noch). Ich fand es in der Tat eine sehr gute Ressource.

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

Für die numerische Berechnung sollte man auf jeden Fall 'numpy' wählen. (Vergewissern Sie sich nur, dass Sie das 'Array' und 'Matrix' und 'Liste' richtig verstanden haben.)

Subodh
quelle