Zalecenie dotyczące metody różnic skończonych w Pythonie naukowym

W przypadku projektu, nad którym pracuję (w hiperbolicznych PDE), chciałbym nieco zorientować się w zachowaniu, patrząc na niektóre wartości liczbowe. Nie jestem jednak zbyt dobrym programistą.

Czy możesz polecić niektóre zasoby do nauki skutecznego kodowania schematów różnic skończonych w Pythonie naukowym (mile widziane są również inne języki o małej krzywej uczenia się)?

Aby przedstawić publiczności (mnie) tę rekomendację:

• Z wykształcenia jestem matematykiem i jestem obeznany z teoretycznymi aspektami schematów różnic skończonych
• Potrzebuję pomocy w tym, jak sprawić, by komputer obliczył to, co chcę, aby obliczył, szczególnie w taki sposób, aby nie powielać zbyt dużego wysiłku już włożonego przez innych (aby nie wymyślać koła ponownie, gdy pakiet jest już dostępny). (Inną rzeczą, której chciałbym uniknąć, jest głupie ręczne kodowanie, gdy istnieją ustalone struktury danych pasujące do celu.)
• Mam trochę doświadczenia w programowaniu; ale nie miałem żadnego w Pythonie (stąd nie mam nic przeciwko, czy istnieją dobre zasoby do nauki innego języka [na przykład Octave na przykład]).
• Przydałyby się zarówno książki, dokumentacja, jak i zbiory przykładowego kodu.

Oto 97-liniowy przykład rozwiązania prostego wielowymiarowego PDE przy użyciu metod różnic skończonych, napisanych przez prof. Davida Ketchesona z repozytorium py4sci, które utrzymuję. W przypadku bardziej skomplikowanych problemów, w których musisz poradzić sobie z wyładowaniami lub konserwacją w skończonej objętości dyskrecji, polecam przyjrzeć się pyclaw , pakietowi oprogramowania, który pomagam opracować.

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

Możesz spojrzeć na Fenics , który jest szkieletem Pythona / C, który pozwala na rozwiązywanie dość ogólnych równań przy użyciu specjalnego języka znaczników. Wykorzystuje jednak głównie elementy skończone, ale warto zajrzeć. Poradnik powinien dać wrażenie, jak łatwo można rozwiązać problemy.

To odniesienie może być dla ciebie bardzo przydatne. To jest otwarta książka w Internecie. Nauczyłem się (wciąż się uczę), python z tej książki. Uważam, że to bardzo dobry zasób.

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

Do obliczeń numerycznych zdecydowanie należy wybrać „numpy”. (po prostu upewnij się, że dobrze zrozumiałeś „tablicę”, „macierz” i „listę”) (zapoznaj się z dokumentacją numpy)