A recursive finite difference implementation

In this entry I will show how to program recursive functions on Python, using as an example the determination of finite difference expressions.

 Theory

The book Computational fluid dynamics, Fourth edition, by Klaus A. Hoffmann and Steve T. Chiang provides recursive definitions for finite difference approximations of derivatives. First they define the forward difference and backward difference operators:

Higher order operators are defined recursively:

Finally they define higher order derivatives around a node i in terms of these higher order operators.

Forward difference

Backward difference

Central difference (n even)

Central difference (n odd)

Implementation

 A recursive definition is a perfect example to be programmed. A Python module finiteDifference.py was implemented (https://drive.google.com/drive/folders/0B7_nEWx0hv56UE9SLWVQUmQyTjQ?usp=sharing). The module contains the basic class definitions needed to obtain higher order finite difference approximations. The basic class is GridPoint, which allows to define a function evaluation at some specified grid location. Here it is __init__ method:

class GridPoint(object):
       def __init__(self, coeff = 1, ind = 0):
....

This allows you to define a gridPoint instance and inspect it using print:

>>> import finiteDifference as fd
>>> a=fd.GridPoint(coeff=2,ind=1)
>>> print(a)
  2.0*f(i+1)

The next class is Den, which is the denominator of a finite difference approximation, it is a very simple class:

class Den(object):
 def __init__(self, power = 1):
...
 

>>> b=fd.Den(power=2)
>>> print(b)
(dh)^2

Next is the Num class,  which is the numerator of a finite difference approximation. It is constructed from a list of instances of the GridPoint class, as the other classes, it can be inspected by using print on it:

class Num(object):
 def __init__(self, exp = None):
...

>>> c=fd.GridPoint(coeff=2,ind=2)
>>> d=fd.GridPoint(coeff=1.0,ind=-1)
>>> e=fd.Num([c,d])
>>> print(e)
  2.0*f(i+2)+  1.0*f(i-1)

The Num by default groups GridPoint instances terms with equal ind parameter:

>>> f=fd.GridPoint(coeff=-1,ind=3)
>>> print(f)
 -1.0*f(i+3)
>>> g=fd.GridPoint(coeff=2,ind=3)
>>> print(g)
  2.0*f(i+3)
>>> h=fd.Num([f,g])
>>> print(h)
  1.0*f(i+3)

Finally, the class FinDiff is the class used to obtain the finite difference approximation. FindDiff makes internal use of the rest of the classes defined in the module (parameter method can be one of backward, forward and central):

class FinDiff(object):
 def __init__(self, method = 'central', order = 1, point = GridPoint()):
...

FindDiff class has as attributes numer and denom that are instances of classes Num and Den, respectively, which must be printed separately:

>>> i=fd.FinDiff(method='forward',order=1,point=fd.GridPoint())
>>> print(i.numer)
  1.0*f(i+1) -1.0*f(i)
>>> print(i.denom)
(dh)^1

Results

Hoffmann and Chiang offer the following tables with the coefficients of the different grid point evaluations involved  in the finite difference representations of the derivatives, for every method, with the order ranging from 1 to 4.

Forward difference

Order	$f_i$	$f_{i+1}$	$f_{i+2}$	$f_{i+3}$	$f_{i+4}$
1	-1	1
2	1	-2	1
3	-1	3	-3	1
4	1	-4	6	-4	1

Lets see the results using FinDiff:

>>> fd1=fd.FinDiff(method='forward',order=1,point=fd.GridPoint())
>>> fd2=fd.FinDiff(method='forward',order=2,point=fd.GridPoint())
>>> fd3=fd.FinDiff(method='forward',order=3,point=fd.GridPoint())
>>> fd4=fd.FinDiff(method='forward',order=4,point=fd.GridPoint())
>>> fd5=fd.FinDiff(method='forward',order=5,point=fd.GridPoint())
>>> print(fd1.numer)
  1.0*f(i+1) -1.0*f(i)
>>> print(fd2.numer)
  1.0*f(i+2) -2.0*f(i+1)+  1.0*f(i)
>>> print(fd3.numer)
  1.0*f(i+3) -3.0*f(i+2)+  3.0*f(i+1) -1.0*f(i)
>>> print(fd4.numer)
  1.0*f(i+4) -4.0*f(i+3)+  6.0*f(i+2) -4.0*f(i+1)+  1.0*f(i)
>>> print(fd5.numer)
  1.0*f(i+5) -5.0*f(i+4)+ 10.0*f(i+3)-10.0*f(i+2)+  5.0*f(i+1) -1.0*f(i)

The results agree up to fourth order and additionally the result for order 5 is shown.

 Backwards difference

Order	$f_{i-4}$	$f_{i-3}$	$f_{i-2}$	$f_{i-1}$	$f_{i}$
1				-1	1
2			1	-2	1
3		-1	3	-3	1
4	1	-4	6	-4	1

>>> bd1=fd.FinDiff(method='backward',order=1,point=fd.GridPoint())
>>> bd2=fd.FinDiff(method='backward',order=2,point=fd.GridPoint())
>>> bd3=fd.FinDiff(method='backward',order=3,point=fd.GridPoint())
>>> bd4=fd.FinDiff(method='backward',order=4,point=fd.GridPoint())
>>> bd5=fd.FinDiff(method='backward',order=5,point=fd.GridPoint())
>>> print(bd1.numer)
  1.0*f(i) -1.0*f(i-1)
>>> print(bd2.numer)
  1.0*f(i) -2.0*f(i-1)+  1.0*f(i-2)
>>> print(bd3.numer)
  1.0*f(i) -3.0*f(i-1)+  3.0*f(i-2) -1.0*f(i-3)
>>> print(bd4.numer)
  1.0*f(i) -4.0*f(i-1)+  6.0*f(i-2) -4.0*f(i-3)+  1.0*f(i-4)
>>> print(bd5.numer)
  1.0*f(i) -5.0*f(i-1)+ 10.0*f(i-2)-10.0*f(i-3)+  5.0*f(i-4) -1.0*f(i-5)

Central difference

Order	$f_{i-2}$	$f_{i-1}$	$f_{i}$	$f_{i+1}$	$f_{i+2}$
1		-0.5	0	0.5
2		1	-2	1
3	-0.5	1	0	-1	0.5
4	1	-4	6	-4	1

>>> cd1=fd.FinDiff(method='central',order=1,point=fd.GridPoint())
>>> cd2=fd.FinDiff(method='central',order=2,point=fd.GridPoint())
>>> cd3=fd.FinDiff(method='central',order=3,point=fd.GridPoint())
>>> cd4=fd.FinDiff(method='central',order=4,point=fd.GridPoint())
>>> cd5=fd.FinDiff(method='central',order=5,point=fd.GridPoint())
>>> print(cd1.numer)
  0.5*f(i+1)+  0.0*f(i) -0.5*f(i-1)
>>> print(cd2.numer)
  1.0*f(i+1) -2.0*f(i)+  1.0*f(i-1)
>>> print(cd3.numer)
  0.5*f(i+2) -1.0*f(i+1)+  0.0*f(i)+  1.0*f(i-1) -0.5*f(i-2)
>>> print(cd4.numer)
  1.0*f(i+2) -4.0*f(i+1)+  6.0*f(i) -4.0*f(i-1)+  1.0*f(i-2)
>>> print(cd5.numer)
  0.5*f(i+3) -2.0*f(i+2)+  2.5*f(i+1)+  0.0*f(i) -2.5*f(i-1)+  2.0*f(i-2) -0.5*f(i-3)