Laplace Equation Solver

by Gennady Stupakov, Stanford Linear Accelerator Center, Stanford
University. Please direct questions and comments to stupakov@wco.com. This note must not be removed.;[s]3:0,0;117,1;134,0;167,-1;2:2,13,9,Times,0,12,0,0,
0;1,13,9,Times,1,12,0,0,0;

Introduction

This notebook discusses how to use the package LESolver.m developed by the
author for solution of two dimensional Laplace equation.

The package LESolver.m contains a Mathematica code for solution of 2D
Laplace equation in a domain enclosed by an arbitrary smooth curve with a
Dirichlet boundary condition. Both interior (in the region encompassed by
the curve) and exterior (outside of the curve) problems can be solved with
the package. However, solving the exterior problem uses a matrix function
SingularValues[] available only in version 3.0 of Mathematica. If you don't
have the version 3.0, you can contact the author and get a MathLink program
that finds singular values in version 2.x and enables the code to solve the
exterior problems with earlier versions. The interior problem can be solved
with versions 2.x without invoking MathLink.

The method used for solution is based on reduction of the Dirichlet problem
to a 1D Fredholm integral equation of the second type. The integral
equation is discretized on a mesh generated on the boundary of the domain
and solved as a system of linear equations. The method assumes that the
boundary is represented by a smooth Jordan curve, however, with a proper
caution, it actually works if the boundary has angles (e.g., rectangle). In
the last case, the solution is apparently not as accurate as for a smooth
boundary with the same number of mesh points.

A few remarks about the hardware requirements are in place. For a serious
calculation you need a Pentium PC, Power Mac or UNIX workstation. Keep in
mind that with the parameter n given you will be solving a linear system of
n*n  equations. I was able to do n as large as 100 on my Power Mac 9500/120
with 32 Mb of RAM. Another possible limitation of the speed comes out from
multiple calculations of the value of the solution in many space points for
graphical presentation of the result. If you want to compute the potential
in 10000 points to plot a nice 3D surface, it might be a time consuming
procedure.

A nice feature of the approach used is that after finding the solution with
the Solver, you can save it, clean the memory, and then continue
calculation of the potential in a new session.

Interior Problem

Finding solution

First, read in the package.

Needs["LESolver`"]

Change this command accordingly if you put the package in a subdirectory of
your package directory.

To illustrate how to use the package, we first chose an example in which
the Laplace equation is solved inside a circle of unit radius. Giving a
physical meaning to the solution, we will be talking about electrostatic
potential inside the circle, determined by its value on the boundary.

First, we need to specify the boundary of the domain as a parametric curve
given as a function of its arc length s

xx[s_]:=Cos[s];
yy[s_]:=Sin[s];

It is important that in all cases the parameter s must represent the arc
length measured along the curve. When parameter s varies from its minimal
value smin to its maximal value smax, (smax>smin) the point passes the
curve in anticlockwise direction. Since the curve is supposed to be closed,
the values of xx and yy at smin must be equal to their respective  values
at smax.

smin=0;
smax=2 Pi//N;

We also have to define the number of points for the mesh on the boundary.
Let us chose

n=40;

Finally, we need to specify the value of the potential on the circular
boundary, again as a function of arc length s, boundpoten[s]. As an
example, we choose the following potential on the boundary

boundpoten[s_]:=Log[(0.5-xx[s])^2+yy[s]^2]

This potential is generated by a line charge with a negative unit linear
charge density located at the position x=0.5, y=0.

Now we need to prepare the input data for the solver. This is done by a
special command from the package, Prepare[], which should be given as
arguments the names of the functions xx and yy, the name of the function
boundpoten for the potential on the boundary, minimal and maximal values of
the parameter s combined in a list, and, finally, the number of points n:

data=Prepare[xx,yy,boundpoten,{0,smax},n];

This function returns a list of parameters which is saved in a variable data.

To visualize the domain in which the solution is searched, one can use a
command DrawPlot[] from the package which takes the argument data,

DrawPlot[data]

The plot shows the boundary with the mesh points indicated by red points
and blue normals at each mesh point pointing inside the area encompassed by
the curve. If you find that the normals point to the outward of the curve,
it means that your curve is passed clockwise and you have to reparametrize
it.

Now we can call the solver and save the solution in the variable solution

solution=Solver[data, Domain->Interior];

The function Solver can take an option Domain which is either Interior, or
Exterior. By setting Domain->Interior  you indicate to the solver that you
are finding the potential in the interior of the region bounded by the
curve; the option Exterior indicates that the potential is searched outside
of the curve. An example of using the Solver to the external solution is
given below. Solution of the interior problem can be performed with
Mathematica versions 2.x,however, solution of the exterior problem requires
using of the function SingularValues[] available only in the version 3.0 of
Mathematica. If you don't have the version 3.0 and you need to solve an
exterior problem, you can contact the author and obtain a MathLink program
which does the singular value decomposition.

After solution is found, calling the function Potential[x,y,solution,data]
returns the value of the potential at the point x, y

Potential[.3,.4, solution,data]

For this particular case, an analytical solution to our problem is known as
a Poisson's integral. In cylindrical coordinates r, fi the analytical
solution is given by the following formula

analytsol[r_,fi_]:=1/(2 Pi//N)*
	NIntegrate[boundpoten[s] (1.0-r^2)/(1+r^2-2 r Cos[fi-s]), {s,0,2 Pi}]

Compare the numerical solution with the analytical one for r=0.5, fi=0
(x=0, y=0.5)

analytsol[0.5,0]

Potential[0.5,0, solution,data]

The numerical solution is accurate to within better than 1%.

Plotting the result

We now want to plot the potential of a charged wire with a unit charge per
unit length located inside the circle. This potential will be given by the
difference of the Potential[] and Log[(0.5-xx[s])^2+yy[s]^2]. Let us define
the function poten[x_,y_] which is equal  to this difference inside the
circle and vanishes outside of it

poten[x_,y_]:=If[x^2+y^2>=1,0,
	Potential[x, y, solution, data]-Log[(0.5-x)^2+y^2]]

and store the values of the potential in the array potenarray

potenarray=Table[poten[x,y],{x,-1.1,1.1,.035},{y,-1.1,1.1,.035}];

Now we plot the result as a 3D surface and show the boundary with
ParametricPlot3D

ParametricPlot3D[{Cos[t],Sin[t],0},{t,0,2 Pi},Boxed->False,
	DisplayFunction->Identity];
ListPlot3D[potenarray,
	ClipFill->None,
	Boxed->False,
	Mesh->False,
	MeshRange->{{-1,1},{-1,1}},
	AxesLabel->{x,y,Potential},DisplayFunction->Identity];
Show[%,%%,DisplayFunction->$DisplayFunction];

We can also use ListContourPlot for plotting the contour levels of the
potential. The red circle show the boundary.

ParametricPlot[{Cos[t],Sin[t]},{t,0,2 Pi},
PlotStyle->{RGBColor[1,0,0]},
DisplayFunction->Identity];

surf=ListContourPlot[potenarray,Contours->15,ContourShading->False,
MeshRange->{{-1,1},{-1,1}},
DisplayFunction->Identity];

Show[%,%%,DisplayFunction->$DisplayFunction]

Exterior problem

Here is an example of exterior problem - a conducting slab with rounded
ends immersed in an external homogeneous vertical electric field. The slab
is parametrically determined by the functions xx[t] and yy[t], where t is
the arc length.

rad=0.2;
pi=N[Pi];
Clear[n,xx,yy,boundpoten]

yy[t_]:=Which[
t<2.,(t-1)/Sqrt[2.]-rad/Sqrt[2.],
t< (2 + pi rad), rad Sin[(t-2-pi/4 rad)/rad]+1/Sqrt[2.],
t< (4 + pi rad), -(t-3-pi rad)/Sqrt[2.]+rad/Sqrt[2.],
t<=(4 + 2 pi rad), rad Sin[(t-4-1 pi/4 rad)/rad]-1/Sqrt[2]
]//N;

xx[t_]:=Which[
t<2.,(t-1)/Sqrt[2.]+rad/Sqrt[2.],
t< (2 + pi rad), rad Cos[(t-2-pi/4 rad)/rad]+1/Sqrt[2],
t< (4 + pi rad), -(t-3-pi rad)/Sqrt[2.]-rad/Sqrt[2.],
t<=(4 + 2 pi rad), rad Cos[(t-4-1 pi/4 rad)/rad]-1/Sqrt[2]
]//N;

tmax=(4 + 2 Pi rad)//N;

The potential on the boundary is equal to the minus x- coordinate

boundpoten[t_]:=-xx[t];

Chose n=80 and plot the slab with the mesh,

n=80;
Clear[data]
data=Prepare[xx,yy,g,{0,tmax},n];
DrawPlot[data]

Find the solution with the option Domain->Exterior

solution=Solver[data, Domain->Exterior];

Here is what stored in the solution,

solution

Plotting Results

To plot the result, we define a function poten[x_,y_] that is equal to the
found potential plus the potential of
the exterior field (equal to x) outside of the slab, and vanishes inside.
The values of this function are saved in the array
potenarray

Clear[poten,potenarray]
poten[x_,y_]:=If[((y>x-rad Sqrt[2.]) &&
	(y<x+rad Sqrt[2.]) &&
	(y<-x+Sqrt[2.]) &&
	(y>-x-Sqrt[2.])) ||
	(((x-1/Sqrt[2.])^2+(y-1/Sqrt[2.])^2)<rad^2) ||
	(((x+1/Sqrt[2.])^2+(y+1/Sqrt[2.])^2)<rad^2)
	,0,Potential[x, y, solution, data]+x]

potenarray=Table[poten[x,y],{x,-1.5,1.5,.05},{y,-1.5,1.5,0.05}];

Now we plot the result as a 3D surface and show the boundary with
ParametricPlot3

ParametricPlot3D[{xx[t],yy[t],0.05},{t,0,tmax},
	Boxed->False,
	DisplayFunction->Identity];
ListPlot3D[potenarray,
	ClipFill->None,
	Boxed->False,
	Mesh->False,
	MeshRange->{{-1.5,1.5},{-1.5,1.5}},
	AxesLabel->{x,y,Potential},DisplayFunction->Identity];
Show[%,%%,DisplayFunction->$DisplayFunction]

Anothe view of the result

ListPlot3D[potenarray,ClipFill->None,Mesh->False,
LightSources->{{{1,0,1},RGBColor[1,1,1]}},
Boxed->False]

We can also use ListContourPlot for plotting the contour levels of the
potential. The red circle shows
the boundary.

ParametricPlot[{xx[t],yy[t]},{t,0,tmax},
	PlotStyle->{RGBColor[1,0,0]},
	DisplayFunction->Identity];

ListContourPlot[potenarray,Contours->20,ContourShading->False,
	MeshRange->{{-1.5,1.5},{-1.5,1.5}},
	DisplayFunction->Identity];

Show[%,%%,DisplayFunction->$DisplayFunction]