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
<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]
SOR[ (* TYPE=S *)
$N -> Max[$JA], (* DATA=I *)
$IA -> In, (* DATA=I[1 + $N] *)
$JA -> In, (* DATA=I[$NZ] *)
$A -> In, (* DATA=R[$NZ] *)
$RHS -> In, (* DATA=R[$N] *)
$U -> 0 :> Out, (* DATA=R[$N] *)
$IWKSP :> Null, (* DATA=I[3*$N] *)
$NW -> $N, (* DATA=I *)
$WKSP :> Null, (* DATA=R[$NW] *)
$IPARAM -> {100, 0, 0, 6, 1, 1, 1, 0, -1, 0, 0, 0},
(* DATA=I[12] *)
$RPARAM -> {5./10^6, 0., 0., 0.75, 1., 0., 0.25,
100./10^16, 0., 0., 0., 0.}, (* DATA=R[12] *)
$IER :> Null, (* DATA=I *)
$NZ -> Last[$IA] - 1 (* DATA=I *)
] (* CODE="LIBRARY" *)
sor[$IA_, $JA_, $A_, $RHS_] -> $U
Import the routine into Mathematica
Import[{SOR}];
InterCall::opened: Opened connection to host solar.
InterCall::import: Importing: {SOR}
InterCall::linked: Using remote driver version 2.0 on host solar.
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
5
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.