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

20

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.
Willie Wong
źródło
Głównym problemem jest to, że nie wiem nawet, od czego zacząć: więc pomocne byłyby nawet podstawowe sugestie.
Willie Wong,
Ograniczeniem jest tylko to, że nie jestem (jeszcze) zaznajomiony z metodami o skończonej objętości; więc będę musiał nauczyć się tej metody w połączeniu. Oczywiście nie sprzeciwiłbym się takiej odpowiedzi.
Willie Wong
PyClaw może obsługiwać nieliniowe terminy źródłowe, ale napisanie własnego solwera Riemanna byłoby skomplikowane, szczególnie w 2 lub wyższych wymiarach. Jeśli chcesz wypróbować prosty schemat różnicowania za pomocą strukturalnych siatek, następną opcją byłoby wypróbowanie czegoś w Petsc4py (Ujawnienie: Jestem również związany z tym projektem), który ma bardziej ogólny cel, a nie równie dobrze- udokumentowane.
Aron Ahmadia
niech nam kontynuować tę dyskusję w czacie
Aron Ahmadia
Cześć Willie (i dla czytelników, którzy nie patrzyli na czat), myślę, że już o tym wiesz, ale skoro wspomniałeś o hiperbolicznych PDE, prawdopodobnie lepiej byłoby, gdybyś zastosował metodę skończonej objętości.
Matthew Emmett

Odpowiedzi:

10

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()
Aron Ahmadia
źródło
8

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.

Moyner
źródło
3

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)

Subodh
źródło