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.