Recommendation for Finite Difference Method in Scientific Python

Question

For a project I am working on (in hyperbolic PDEs) I would like to get some rough handle on the behavior by looking at some numerics. I am, however, not a very good programmer.

Can you recommend some resources for learning how to effectively code finite difference schemes in Scientific Python (other languages with small learning curve also welcome)?

To give you an idea of the audience (me) for this recommendation:

I am a pure mathematician by training, and am somewhat familiar with the theoretical aspects of finite difference schemes
What I need help with is how to make the computer compute what I want it to compute, especially in a way that I don't duplicate too much of the effort already put in by others (so as to not re-invent the wheel when a package is already available). (Another thing I would like to avoid is to stupidly code something by hand when there are established data structures fitting the purpose.)
I have had some coding experience; but I have had none in Python (hence I don't mind if there are good resources for learning a different language [say, Octave for example]).
Books, documentation would both be useful, as would collections of example code.

The main problem is that I don't even know where to start looking: so even basic suggestions would be helpful. — Willie Wong, May 16 '12 at 12:53
The restriction is only that I am not (yet) familiar with finite volume methods; so I will have to learn the method in conjunction. I would not object to such an answer, of course. — Willie Wong, May 16 '12 at 13:12
PyClaw can handle nonlinear source terms, but writing your own Riemann solver would be complicated, particularly in 2nd or higher dimensions. If you want to try a simple finite-differencing scheme with structured grids, your next option would be to try something in petsc4py, (Disclosure: I'm affiliated with this project as well), which is more general purpose and not as well-documented. — Aron Ahmadia, May 16 '12 at 13:51
Hi Willie (and for readers that haven't looked at the chat), I think you already know this, but since you mentioned hyperbolic PDEs you'd probably be better off with a finite volume method. — Matthew Emmett, May 16 '12 at 20:43
Hi, Willie, finite differences doesn't work quite well for hyperbolic PDEs especially when it is nonlinear, even it couldn't capture the true solution of a simple transport equation where discontinuities are present, I wonder by any chance if you could share with us a minimal example of the equation you want to solve. — Shuhao Cao, Sep 18 '12 at 15:30
@Shuhao: discontinuities will not be a problem (that I need to worry about), see my discussion with Araon Ahmadia linked above. — Willie Wong, Sep 24 '12 at 09:00
See New technique for solving time fractional wave equation: Python — KRISHNA GHODE, Aug 05 '21 at 11:03

Aron Ahmadia · Accepted Answer · 2012-05-16T14:40:51.957

Here is a 97-line example of solving a simple multivariate PDE using finite difference methods, contributed by Prof. David Ketcheson, from the py4sci repository I maintain. For more complicated problems where you need to handle shocks or conservation in a finite-volume discretization, I recommend looking at pyclaw, a software package that I help develop.

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

score 10 · Answer 2 · answered May 16 '12 at 14:18

You could have a look at Fenics, which is a python/C framwork which allows quite general equations to be solved using a special markup language. It mostly uses finite elements though, but worth a look. The tutorial should give you an impression of how easy it can be to solve problems.

score 3 · Answer 3 · answered Jun 03 '12 at 09:08

This reference might be very useful for you. This is an open book on internet. I learned (still am learning), python from this book. I found it very good resource indeed.

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

For numerical calculation, One should definitely go for 'numpy'. (just make sure that you have understood the 'array' and 'matrix' and 'list' properly) (refer numpy documentation for that)

Recommendation for Finite Difference Method in Scientific Python

3 Answers3

Linked