Recomendação para o método das diferenças finitas em Python científico

20

Para um projeto em que estou trabalhando (em PDEs hiperbólicas), gostaria de ter uma idéia mais detalhada do comportamento, analisando alguns números. No entanto, não sou um programador muito bom.

Você pode recomendar alguns recursos para aprender a codificar efetivamente esquemas de diferenças finitas no Scientific Python (outros idiomas com pequena curva de aprendizado também são bem-vindos)?

Para lhe dar uma idéia do público (eu) para esta recomendação:

  • Sou um matemático puro por treinamento e estou um pouco familiarizado com os aspectos teóricos dos esquemas de diferenças finitas
  • O que eu preciso de ajuda é como fazer o computador calcular o que eu quero, especialmente de uma maneira que não duplique muito o esforço já realizado por outras pessoas (para não reinventar a roda quando um pacote já está disponível). (Outra coisa que eu gostaria de evitar é codificar algo estupidamente à mão quando houver estruturas de dados estabelecidas adequadas ao objetivo.)
  • Eu tive alguma experiência em codificação; mas eu não tive nenhum em Python (portanto, não me importo se existem bons recursos para aprender um idioma diferente [digamos, Octave por exemplo]).
  • Livros, a documentação seria útil, assim como coleções de código de exemplo.
Willie Wong
fonte
O principal problema é que eu nem sei por onde começar a procurar: portanto, mesmo sugestões básicas seriam úteis.
Willie Wong
A restrição é apenas que eu ainda não estou familiarizado com métodos de volume finito; então terei que aprender o método em conjunto. Eu não contestaria essa resposta, é claro.
Willie Wong
O PyClaw pode lidar com termos de fontes não lineares, mas escrever seu próprio solucionador Riemann seria complicado, principalmente em 2ª ou mais dimensões. Se você quiser tentar um esquema simples de diferenciação finita com grades estruturadas, sua próxima opção seria tentar algo em petsc4py (Divulgação: também sou afiliado a esse projeto), que é de uso geral e não tão adequado. documentado.
Aron Ahmadia
Olá Willie (e para leitores que não viram o bate-papo), acho que você já sabe disso, mas desde que você mencionou PDEs hiperbólicos, provavelmente se sairia melhor com um método de volume finito.
Matthew Emmett

Respostas:

10

Aqui está um exemplo de 97 linhas de solução de um PDE multivariado simples usando métodos de diferenças finitas, contribuído pelo Prof. David Ketcheson , do repositório py4sci que mantenho. Para problemas mais complicados nos quais você precisa lidar com choques ou conservação com uma discretização de volume finito, recomendo examinar o pyclaw , um pacote de software que eu ajudo a desenvolver.

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

Você pode dar uma olhada no Fenics , que é um framework em python / C que permite que equações bastante gerais sejam resolvidas usando uma linguagem de marcação especial. No entanto, usa principalmente elementos finitos, mas vale a pena dar uma olhada. O tutorial deve lhe dar uma impressão de como pode ser fácil resolver problemas.

moyner
fonte
3

Esta referência pode ser muito útil para você. Este é um livro aberto na internet. Eu aprendi (ainda estou aprendendo), python deste livro. Achei muito bom recurso.

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

Para o cálculo numérico, deve-se definitivamente usar 'numpy'. (apenas certifique-se de ter entendido as 'matriz' e 'matriz' e 'lista' corretamente) (consulte a documentação numpy para isso)

Subodh
fonte