# Recommendation for Finite Difference Method in Scientific Python

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. May 16, 2012 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. May 16, 2012 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. May 16, 2012 at 13:51
• May 16, 2012 at 13:53
• 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. May 16, 2012 at 20:43

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


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.

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)