Rekomendasi untuk Metode Perbedaan Hingga dalam Python Ilmiah

20

Untuk proyek yang sedang saya kerjakan (dalam PDE hiperbolik), saya ingin menangani perilaku tersebut dengan melihat beberapa angka. Namun, saya bukan programmer yang sangat baik.

Bisakah Anda merekomendasikan beberapa sumber untuk mempelajari cara kode skema perbedaan hingga secara efektif dalam Scientific Python (bahasa lain dengan kurva belajar kecil juga diterima)?

Untuk memberi Anda gambaran tentang audiens (saya) untuk rekomendasi ini:

  • Saya seorang ahli matematika murni dengan pelatihan, dan saya agak akrab dengan aspek teoritis dari skema beda hingga
  • Apa yang saya perlu bantuan adalah bagaimana membuat komputer menghitung apa yang saya inginkan, terutama dengan cara yang saya tidak menduplikasi terlalu banyak upaya yang sudah dilakukan oleh orang lain (agar tidak menciptakan kembali roda saat paket sudah tersedia). (Hal lain yang ingin saya hindari adalah dengan bodohnya mengkodekan sesuatu dengan tangan ketika ada struktur data yang sesuai dengan tujuannya.)
  • Saya telah memiliki beberapa pengalaman pengkodean; tapi saya tidak punya di Python (maka saya tidak keberatan jika ada sumber daya yang baik untuk belajar bahasa yang berbeda [katakanlah, Oktaf misalnya]).
  • Buku-buku, dokumentasi keduanya akan berguna, seperti koleksi contoh kode.
Willie Wong
sumber
Masalah utama adalah bahwa saya bahkan tidak tahu harus mulai dari mana: bahkan saran dasar pun akan sangat membantu.
Willie Wong
Pembatasannya hanya bahwa saya belum (belum) terbiasa dengan metode volume yang terbatas; jadi saya harus belajar metode dalam hubungannya. Saya tidak akan keberatan dengan jawaban seperti itu, tentu saja.
Willie Wong
PyClaw dapat menangani istilah sumber nonlinier, tetapi menulis pemecah Riemann Anda sendiri akan rumit, terutama dalam dimensi ke-2 atau lebih tinggi. Jika Anda ingin mencoba skema pembedaan terbatas sederhana dengan grid terstruktur, opsi Anda berikutnya adalah mencoba sesuatu di petsc4py , (Pengungkapan: Saya juga berafiliasi dengan proyek ini), yang merupakan tujuan yang lebih umum dan tidak juga didokumentasikan.
Aron Ahmadia
mari kita lanjutkan diskusi ini dalam obrolan
Aron Ahmadia
Hai Willie (dan untuk pembaca yang belum melihat obrolan), saya pikir Anda sudah tahu ini, tetapi karena Anda menyebutkan PDE hiperbolik Anda mungkin akan lebih baik dengan metode volume yang terbatas.
Matthew Emmett

Jawaban:

10

Berikut adalah contoh 97-garis penyelesaian PDE multivariat sederhana menggunakan metode beda hingga, disumbangkan oleh Prof. David Ketcheson , dari repositori py4sci yang saya pelihara . Untuk masalah yang lebih rumit di mana Anda perlu menangani guncangan atau konservasi dalam diskritisasi volume terbatas, saya sarankan melihat pyclaw , paket perangkat lunak yang saya bantu kembangkan.

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

Anda bisa melihat Fenics , yang merupakan kerangka python / C yang memungkinkan persamaan yang cukup umum untuk diselesaikan menggunakan bahasa markup khusus. Ini sebagian besar menggunakan elemen hingga, tetapi patut dilihat. The tutorial harus memberikan kesan betapa mudahnya dapat untuk memecahkan masalah.

moyner
sumber
3

Referensi ini mungkin sangat berguna bagi Anda. Ini adalah buku terbuka di internet. Saya belajar (masih belajar), python dari buku ini. Saya menemukan sumber daya yang sangat bagus.

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

Untuk perhitungan numerik, seseorang harus menggunakan 'numpy'. (hanya pastikan bahwa Anda telah memahami 'array' dan 'matrix' dan 'daftar' dengan benar) (lihat dokumentasi numpy untuk itu)

Subodh
sumber