A Matlab-Based Finite Difierence Solver for the Poisson ...

1 A Matlab-Based Finite Difference Solver for the Poisson Problemwith Mixed Dirichlet-Neumann Boundary ConditionsAshton S. Reimera), Alexei F. Cheviakovb)Department of Mathematics and Statistics, University of Saskatchewan, Saskatoon, S7N 5E6 CanadaApril 17, 2012 AbstractA Matlab-Based Finite -difference numerical Solver for the Poisson equation for a rectangle anda disk in two dimensions, and a spherical domain in three dimensions, is presented. The Solver isoptimized for handling an arbitrary combination of Dirichlet and Neumann boundary conditions,and allows for full user control of mesh refinement. The Solver routines utilize effective andparallelized sparse vector and matrix operations. Computations exhibit high speeds, numericalstability with respect to mesh size and mesh refinement, and acceptable error values even ondesktop SUMMARY/NEW VERSION PROGRAM SUMMARYM anuscript Title:A Matlab-Based Finite -difference Solver for the Poisson Problem in Two and Three Di-mensions, Optimized for Strongly Heterogeneous Boundary ConditionsAuthors:Ashton S.

2 Reimer and Alexei F. CheviakovProgram Title:FDMRP - Finite -Difference with Mesh Refinement Poisson equation solverJournal Reference:Catalogue identifier:Licensing provisions:GNU General Public License language:Matlab 2010aComputer:PC, MacintoshOperating system:Windows, OSX, LinuxRAM:8GB (8,589,934,592 bytes)Number of processors used:Single Processor, Dual CoreKeywords: Poisson problem, Finite -difference Solver , Matlab, Strongly heterogeneous boundary conditions,Narrow Escape Differential EquationsNature of problem:To solve the Poisson problem in a standard domain with patchy surface -type (stronglyheterogeneous) Neumann/Dirichlet boundary method: Finite difference with mesh :Spherical domain in 3D; rectangular domain or a disk in features:Choice between mldivide/iterative Solver for the solution of large system of linear algebraicequations that arise.

3 Full user control of Neumann/Dirichlet boundary conditions and mesh time:Depending on the number of points taken and the geometry of the domain, the routine maytake from less than a second to several hours to IntroductionThe Laplacian operator is a fundamental component in a large number of multi-dimensional linearmodels of mathematical physics; to name a few, the Laplacian arises in description of variousa)Electronic mail: mail: and quantum wave phenomena, as well as models involving diffusion and viscosity time-independent boundary value problems (BVP) for Laplace and Poisson equations con-stitute a special class of important problems. The Laplace operator is separable in many classicaland esoteric coordinate systems [7,8]. The corresponding linear problems can be solved analyticallyfor specific domains which correspond to parallelepipeds in the respective coordinate systems, pro-vided that boundary conditions are sufficiently simple.

4 For a general BVP, however, one normallyhas to rely on numerical methods are common numerical methods for solution of linear second-order time-independent partial differential equations. Finite -difference methods involve discretization of thespatial domain, the differential equation, and boundary conditions, and a subsequent solution of alarge system of linear equations for the approximate solution values in the nodes of the numericalmesh. Simplest discretizations of second-order differential operators commonly have first or secondorder accuracy. The resulting large system of linear equations involves a sparse matrix and aresolved by iterative methods (Jacobi, Gauss-Seigel, etc.) or Gaussian elimination/LU decomposition,which have been significantly optimized for sparse current work is motivated by BVPs for the Poisson equation where boundary correspond toso called patchy surfaces , , are strongly heterogeneous, involving combination of Neumann andDirichlet boundary conditions on different parts of the boundary.

5 Models involving patchy surfaceBVPs are found in various fields. For example, models of vacuum tubes with patchy surfacesin electrostatics were considered in [5, 6]. Cell signalling models in biology are concerned withmotion of ions that seek an open channel located in the cell membrane; the latter is modelled as acombination of homogeneous reflecting (Neumann) and absorbing (Dirichlet) boundary conditions( , [1,10]). In particular, the continuum limit, the mean first passage time (MFPT)v(x) requiredfor a Brownian particle to leave a domain with a reflecting boundary and absorbing traps on theboundary satisfies the mixed Dirichlet-Neumann problem (cf. [3]) v= 1D, x ,v= 0, x a= Nj=1 j, j= 1, .., N; nv= 0, x r,( )whereDis the diffusion coefficient and may be constant or a function of space. Dirichlet boundaryparts jcorrespond to small absorbing traps, and the boundary part rcorresponds to reflective(Neumann) boundary conditions.

6 In order to determine the applicability limits of the approximatesolutions derived in [2] and [9], a numerical solution to the MFPT problem was required for com-parison. The numerical Solver described below was initially developed to solve MFPT problems( ).The programs presented in the current work are rather general; they provide numerical solutionsof Poisson boundary value problems v=f(x), x ,v=vD(x), x D, nv=vN(x), x N,( )involving the combination of Dirichlet boundary conditions on the part Dof the boundary, andNeumann boundary conditions on the part Nof the boundary. The forcing termf(x) is arbitrary; is either a rectangle or a disk in two dimensions, or a sphere in three programs were implemented inMatlab, and are respectively optimized to avoid loops anduse vector and matrix operations even at the stage of problem formulation, which allowed for asignificant speedup of computation.

7 The programs are compatible withOctave, the open-sourceMatlabequivalent, with the exception that the iterative solving options are unavailable many general purpose Finite -element software packages exist ( ,ELMER,COMSOL, etc.), such2 black-box packages require a significant investment of time and/or money, not offering sufficientinsight into the numerical methods used and their implementation. TheMatlab-based numericalsolvers described in the current contribution offer a transparent, simple-to-use way to solve Poissonproblems in simple geometries with a Finite -difference , a second-order symmetric discretization of the Laplacian operator was used. Theprograms allow for mesh refinement through variable step sizes in any direction (in which case thescheme becomes generally first-order accurate). Mesh refinement allows one to concentrate themesh around given point(s) within the computation domain or on its boundary, where higherresolution is needed.

8 For example, for MFPT problems with small boundary or volume traps, meshrefinement was used to resolve the trap regions accurately, and have coarser mesh far from trapswhere the solution changes slowly. The properties of the Finite -difference approximation, includingimplications of singularities in curvilinear coordinates, are discussed in Section a two- or three-dimensional domain, the discretization of the Poisson BVP ( ) yields a systemof sparse linear algebraic equations containingN=LMequations for two-dimensional domains,andN=LMNequations for three-dimensional domains, whereL, M, Nare the numbers of steps inthe corresponding directions. In order to solve the large sparse system, the programs take advantageof highly effectiveMatlabroutines:mldivideandbicg. For example, on a system with 8 gigabytesof memory, for the spherical domain, it was possible to use a total ofN 107data points,using the valuesL= 400, M=N= 200 (in the directions of the polar angle, radius, and azimuthalangle, respectively).

9 Depending on the number of points used and hardware configurations, theruntime can vary from fractions of a second to examples for problems in a rectangle, disk, and sphere with homogeneous and nonhomoge-neous boundary conditions are presented in Sections 3 5. The examples include error plots for trialsolutions which were known exactly. For all computed configurations where exact solutions are notknown, stability of the numerical results with respect to mesh size and mesh refinement functionswas experimentally presented programs can be easily generalized to other sets of curvilinear coordinates usingsimilar considerations and a corresponding different set of scaling (Lam e) The Finite -Difference ApproximationFinite-difference methods approximate the derivative of a function at a given point by a finitedifference. In one dimension,a x b, consider a functionu(x) and a numerical mesha=x1<.

10 < xn=b. Letu(xi) ui. The mesh step sizes are given byhi=xi+1 xi. A derivative ofu(x)atx=xican be approximated by, for example, a forward difference( u x)i=ui+1 uihi( )or a central difference( u x)i=ui+1 ui 1hi+hi 1,( )and the second derivative by a central difference( 2u x2)i=(2hi+hi 1)(ui+1 uihi ui ui 1hi 1).( )Formulas ( ) and ( ) have second-order accuracy if allhi=h= const, and are first-orderaccurate now present the discretized versions of the Dirichlet and Neumann boundary conditions andthe Poisson equations operator in Cartesian, polar and spherical Approximation of Boundary ConditionsFor a Dirichlet boundary condition, ,u(b) =A, one one simply sets the approximate numericalsolution at a given boundary point to equal the boundary value:u1= order to apply Neumann boundary conditions of the form u x(a) =Ain one dimension, thesimplest choice is to use a first-order approximation given by a backward difference, and hence onehasu1=u2 formula naturally generalizes to multiple dimensions and curvilinear geometries, where aNeumann boundary condition corresponds to the value of the normal derivative prescribed at agiven boundary it is well known, solutions to Poisson equation with a mixture of Dirichlet and Neumann bound-ary conditions can involve mild, square-root-type singularities, such that the solution derivative canbecome unbounded at a boundary point where the boundary condition changes from Dirichlet toNeumann (see, , [4]).