Solving an elliptic partial differential equation using InterCall

Abstract

Numerically solving elliptic partial differential equations often requires solving a very large, but sparse, linear system. This notebook contains an example in which InterCall is used to solve a PDE by accessing the ITPACK public domain library to handle the sparse linear system. Arbitrary two dimensional regions (including holes and cracks) as well as rectangular three dimensional regions can also be handled with this technique. Some exercises are given.

ITPACK is a public domain fortran library that contains various iterative methods for solving sparse linear systems. It is particularly useful for matrix problems that arise from elliptic partial differential equations with Dirichlet, Neuman or mixed boundary conditions on arbitrary two dimensional regions (including holes and cracks) as well as rectangular three dimensional regions.

The problem

An elliptic partial differential equation

Solve uxx + 2 uyy = 0 on S = [0,1] X [0,1] with  u = 1 + sin(π x) y  on the boundary of S.
Use the standard 5-point central finite difference formula, ie:
                              -2
                          -1  6 -1
                              -2
This will lead to a system involving a large positive definite sparse matrix with just 5 nonzero elements per row. The sparse system can be solved using an iterative method such as successive over-relaxation.

Various iterative methods are contained in the ITPACK library, and these can be accessed from Mathematica using InterCall.

The method

One method for solving an elliptic partial differential equations is to discretize the equation and then solve the resulting large, but sparse, matrix system. Implementing a sparse linear system solver in Mathematica would however be slow and so we shall import a fortran routine called SOR from the ITPACK library using InterCall to solve the system.

Startup InterCall and import ITPACK

Load in the InterCall package

<<InterCall.m

Declare defaults for the ITPACK routine that does successive over-relaxation

Read the ITPACK manual to understand the arguments and defaults below.

AddDefault[
   SOR[ $N->Max[$JA],$IA->In,$JA->In,$A->In,$RHS->In,
        $U->0:>Out,$IWKSP:>Null,$NW->$N,$WKSP:>Null,
        $IPARAM->{100,0,0,6,1,1,1,0,-1,0,0,0},
        $RPARAM->{5.*10^-6,0.,0.,.75,1.,0.,.25,
                  100.*10^-16,0.,0.,0.,0.},
        $IER:>Null,$NZ->Last[$IA]-1],
   S[I,I[$N+1],I[$NZ],R[$NZ],R[$N],R[$N],I[3*$N],I,
     R[$NW],I[12],R[12],I,I]
         ];

Check the default setting

The default shows that we must supply the positions and values of the non-zero elements, and also the rhs. The solution u is returned.

GetDefault[sor]

Import the routine into Mathematica

Import[{SOR}];

Use a 19x19 grid for the solution

This will lead to a sparse matrix with 1729 nonzero elements and 361 unknowns.

n = 19; nx = ny = n; h = 1./(n+1);

Setup the sparse matrix and right hand side

The sparse matrix will be represented by the position and value of the nonzero elements.

{NORTH,  WEST, CENTRE, EAST,   SOUTH} = {1,2,3,4,5};
pos = Table[{{i,j+1},{i-1,j},{i,j},{i+1,j},{i,j-1}},{i,1,nx},{j,1,ny}];
val = Table[{  -2,     -1,     6,    -1,     -2   },{i,1,nx},{j,1,ny}];
rhs = Table[                   0,                   {i,1,nx},{j,1,ny}];

Setup the boundary conditions

The following adjusts the values of the matrix and sets the rhs along the boundary of the region.

f[x_,y_] := N[ 1 + Sin[Pi x] y ];
Do[ (* along the NORTH wall *) x = i h; y = 1.;
   rhs[[i,ny      ]] += -val[[i,ny,NORTH]]*f[x,y];
   pos[[i,ny,NORTH]]=0;  val[[i,ny,NORTH]]=0., {i,1,nx}];
Do[ (* along the SOUTH wall *) x = i h; y = 0.;
   rhs[[i, 1      ]] += -val[[i, 1,SOUTH]]*f[x,y];
   pos[[i, 1,SOUTH]]=0;  val[[i, 1,SOUTH]]=0., {i,1,nx}];
Do[ (* along the  EAST wall *) x = 1.; y = j h;
   rhs[[nx,j      ]] += -val[[nx,j, EAST]]*f[x,y];
   pos[[nx,j, EAST]]=0;  val[[nx,j, EAST]]=0., {j,1,ny}];
Do[ (* along the  WEST wall *) x = 0.; y = j h;
   rhs[[ 1,j      ]] += -val[[ 1,j, WEST]]*f[x,y];
   pos[[ 1,j, WEST]]=0;  val[[ 1,j, WEST]]=0., {j,1,ny}];

Flatten out various values to be passed to the SOR routine

Consult the ITPACK manual to find out the meaning of the arguments to SOR.

ja1= Flatten[pos /. {i_,j_} -> (i-1)*ny+j, 1];
ia = FoldList[Plus, 1, Length /@ ja1];
ja = ja1//Flatten;
a  = val//Flatten;
b  = rhs//Flatten;

Solve the sparse system using ITPACK's SOR routine

This takes about 10 seconds of real time.

u = Partition[sor[ia, ja, a, b], ny];

Append the known boundary values to the solution

u = MapIndexed[
        Module[{i=#2[[1]]-1},Join[{f[i h,0.]},#1,{f[i h,1.]}]]&,
        Join[{Table[f[0.,j h],{j,1,ny}]}, u,
             {Table[f[1.,j h],{j,1,ny}]}] ];

Plot the solution

ListPlot3D[ u//Transpose, AxesLabel->{"x","y","u"},
                          MeshRange->{{0,1},{0,1}} ];

Check the relative accuracy of the solution

The argument RPARAM contains pieces of information about the performance of SOR. In particular its eleventh component gives the number of digits of relative accuracy in solution.

Peek[$RPARAM][[11]]//Round

Exercises

Modify the above method to solve the same problem on a region with a square hole removed. Take u=0 along the square bounded by x=0.25, x=0.5, and y=0.5, y=0.75 . Here is the solution you should get:

How would you remove a circular hole?

Try a problem with mixed Dirichlet and Neuman boundary conditions.

Try a 3-dimensional problem. Think of a way to display the answer.

References

ITPACK

ITPACK is a public domain fortran library written by David R. Kincaid, John R. Respess, and David M. Young of the University of Texas at Austin and Roger G. Grimes of the Boeing Computer Services Company. It is avaliable via netlib.

InterCall

InterCall was developed by Terry Robb. For more information on InterCall send email to intercall_forward@wri.com. An order form is available via MathSource – email the line send 0202-587-0011 to mathsource@wri.com for details. Information can also be posted or faxed to Analytica below.

Analytica                     Fax: +61 (9) 386 5666
P.O. Box 522
Nedlands 6009
Perth WA
AUSTRALIA.

Appendix

There is no appendix.