PDE Relaxation

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2007

Introduction

Partial differential equations (PDE) are usually much more difficult to solve than ordinary differential equations and there has accumulated a vast literature on many different symbolic and numerical techniques for attacking various classes of PDEs.  We obviously cannot review and compare these techniques in a course that covers other topics as well.  For our purposes it is sufficient to discuss one simple technique that is applicable to many of the PDEs encountered in physics.

Rather than attempt a general discussion, we consider the specific problem of Laplace's equation

\(TraditionalForm\`∂\^2 ϕ\/∂x\^2 + ∂\^2 ϕ\/∂y\^2 == 0\)

with ϕ[x, y] specified on a set of internal surfaces and the boundary to a two-dimensional region.  We will attempt to solve the problem using a uniform rectangular grid

(Subscript[x, i]) = (Subscript[x, 0]) + i δ
(Subscript[y, j]) = (Subscript[y, 0]) + j δ
(Subscript[ϕ, i,j]) = ϕ[(Subscript[x, i]), 
        (Subscript[y, j])]

where δ is the grid spacing.  To discretize the equation we need to express second derivatives in terms of finite differences.  The first derivative is represented by the discrete rule:

rule1 = 
        (Subscript[f, i])_ -> (((Subscript[f, i + 1]) - (Subscript[f, i - 1]))/(2  δ));

Applying this rule twice,

(Subscript[f, i]) /. rule1 /. rule1 // Simplify

we find that the second derivative is represented by a simple formula based upon a step size Δ -> 2  δ.  Notice that the Laplacian only relates grid points whose indices differ by even numbers.  Hence, the discrete Laplacian operator in two dimensions is represented as

DiscreteLaplacian2D[ϕ_, i_, 
        j_] := ((ϕ[[i + 1, j]] + ϕ[[i - 1, j]] + ϕ[[i, j + 1]] + ϕ[[i, j - 1]])/4) - ϕ[[i, j]]

where the grid has been rescaled.  Therefore, the discrete Laplacian at a particular grid point is represented by the difference between the value of the function at that point and the average of its values at the four nearest grid points.  In this sense one can describe a solution to Laplace's equation as the function which interpolates between boundary values as smoothly as possible.

A simple algorithm, called the relaxation method, now presents itself.  We start by establishing an array ϕ[[i,j]] that satisfies the boundary conditions and contains reasonable guesses for the function at grid points not on the boundary surfaces.  We then systematically replace each point outside the boundary surfaces by the average of its four nearest neighbors.  This process is repeated until the solution converges to acceptable accuracy.  As the local averaging procedure is iterated, errors in the initial guesses are reduced as the solution gradually relaxes toward the correct solution to the discretized version of Laplace's equation.  Obviously, the smaller the grid spacing and the greater the number of grid points the more accurately will the discretized solution represent the continuous solution, at a cost in computation time proportional to (Δ^(-2)).  Errors in the initial solution can be described as tension in the proposed solution.  As the iteration proceeds, these errors are smoothed out and the function relaxes toward the true solution.  Corrections diffuse through the grid slowly, spreading to nearest neighbors in each iteration.

Unfortunately, like most simple algorithms, the direct relaxation method converges too slowly to be of much practical utility.  Fortunately, a relatively simple modification improves the performance of the relaxation method considerably.  The trick is to anticipate the changes in the local average that will result from changes to neighboring values during the relaxation process and thus to overcorrect each point.  The resulting algorithm is known as simultaneous over-relaxation (SOR).  One can show that for the optimal choice of the over-relaxation parameter, ω, the computation time becomes linear in (Δ^(-1)); hence, the SOR is faster by about a factor of 100 for a typical 100×100 grid.  The basic SOR algorithm is now

ϕ[[i, j]]->ϕ[[i, 
          j]] + ω DiscreteLaplacian[ϕ, i, 
            j]

where ω is adjusted as the iteration progresses.  The primary difficulty in implementing the SOR algorithm is optimization of the over-relaxation parameter — the gain in efficiency is realized only when ω is rather close to its optimal value.

The over-relaxation parameter is based upon a quantity

\(\(\(	\)\(ρ\_Jacobi = 
      1 - π\^2\/\(2  n\^2\)\)\)\)

known in PDE theory as the spectral radius of Jacobi iteration (a simpler algorithm), where n×n is the dimensionality of the grid (assumed square).  Rather than process the grid points sequentially, it is more efficient to divide each iteration into two halves, adjusting the even points first followed by the odd points.  For the first half-iteration, we use ω = 1.  For the second half of the first iteration, we use ω -> 1/(1 - (ρ^2)/2).  For all subsequent iterations, ω -> 
        1/(1 - (ρ^2) ω/4).  One can then show that as the iteration proceeds, the over-relaxation parameter approaches the optimal value

\(ω -> 
      2\/\(\(1 + \@\(1 - ρ\_Jacobi\%2\)\)\(	\)\)\)
(Subscript[ρ, SOR]) = (((((Subscript[ρ, Jacobi])/(1 + (Sqrt[1 - (Subscript[ρ, Jacobi]^2)])))))^2)

This strategy is known as Chebyshev acceleration of the SOR algorithm.  For further details about the SOR method and references to original literature see Press et al., Numerical Recipes.

Initialization

ClearAll["Global`*"];
Off[General::spell,  General::spell1];
$TextStyle = {FontFamily -> "Times", 
          FontSize -> 12};
Needs["Graphics`"]

Implementation

We assume that the problem is formulated on a square array and that the boundary surfaces are parallel to the coordinate axes.  We must supply an array of initial values for the potential function and an array of True|False values which specify the points that lie on the boundary — the potential for points on the boundary does not change.  The relaxation algorithm is iterated until either the maximum number of iterations is exceeded or until the root-mean-square Laplacian

√((1/n) ∑(((φ[[i, j]] - (φ[[i + 1, j]] + φ[[i - 1, j]] + φ[[i, j + 1]] + φ[[i, j - 1]])/4. ))^2)) < conv

falls below a specified convergence criterion.

\(\(\(Clear[relax];\)\(
\)
    \)\)
\(relax[ϕ_, boundary_, nmax_, conv_] := 
      Module[{φ, n, m, ρ, ω, ξ, niter, 
          test}, 
{n, m} = Dimensions[ϕ]; ρ = 
          1.  - π\^2\/\(2. \ n\^2\); ω = 
          1. ; φ = ϕ; niter = 0; 
        test = 100. ; 
While[
          test > conv\  && \ niter < 2  nmax, test = 0. ; 
          Do[If[Not[boundary[[i, j]]] && 
                Mod[i + j, 2] == 
                  Mod[niter, 2], ξ = φ[[i, 
                    j]] - \((φ[[i + 1, 
                          j]] + φ[[i - 1, 
                          j]] + φ[[i, j + 
                            1]] + φ[[i, j - 1]])\)/4. ; 
              test = test + ξ\^2; φ[[i, 
                  j]]\  = φ[[
                    i, j]] - ω\ ξ], {i, 2, 
              n - 1}, {j, 2, m - 1}]; niter = niter + 1; 
          If[niter == 1, ω = 
              1. /\((1.  - ρ\^2/2. )\), ω = 
              1. /\((1.  - \(ρ\^2\) ω/4. )\)]; 
          test = \@\(test\/n\)]; 
        Print["number of iterations = ", niter/2]; 
        Print["residual = ", \ test]; Return[φ]]\)

Example: parallel-plate capacitor in two dimensions

The electric field is constant inside an infinite parallel-plate capacitor and vanishes outside.  The field for a capacitor of finite dimensions is more complicated, bowing outward near the edges and connecting the outer edges with wide, somewhat circular, arcs.  For simplicity we consider a finite parallel-plate capacitor in two dimensions, assuming that it remains very long in the third dimension.

Initialize potential

First, we establish 100×100 grids for the initial potential and the boundary specification.

Clear[boundary, ϕ, ϕ0];
boundary = Table[False, {i, 100}, {j, 100}];
ϕ = Table[0.0, {i, 100}, {j, 100}];

Next we declare the two plates to have constant potentials

Do[If[
        i ≥ 31 && 
          i ≤ 70 && (j == 41  || j == 60), 
        boundary[[i, j]] = True, 
        boundary[[i, j]] = False], {i, 
        1, 100}, {j, 1, 100}]

and the initial potential to vary linearly between them.  Some care in the choice of coordinates is needed to center the plates within the grid.

Do[If[
        i ≥ 31 && i ≤ 70 && j ≥ 41 && 
          j ≤ 60, ϕ[[i, 
            j]] = (2  j - 101)/
            19. , ϕ[[i, j]] = 
          0.0], {i, 1, 100}, {j, 1, 100}]

A better initial guess would include the region outside the plates, but we do not bother guessing there.

Generate relaxed solution

Now we simply invoke the relaxation function and ask it to produce a fairly accurate solution for the entire region.  Be patient — this can take a lot of time!

ϕ0 = ϕ;

ϕ = relax[ϕ0, boundary, 1000, (10^(-4))];

Display potential and field

It is instructive to examine orthogonal slices across the capacitor

ListPlot[
        Table[ϕ[[50, j]], {j, 1, 
            100}], PlotRange -> All];

ListPlot[
        Table[ϕ[[i, 45]], {i, 1, 
            100}]];

and to make a three-dimensional plot.

ListPlot3D[ϕ, PlotRange -> All];

It is also useful to construct an interpolation function which can be used to produce a contour plot.

f = ListInterpolation[ϕ];

ContourPlot[f[x, y], {x, 2, 99}, {y, 2, 99}, PlotPoints -> 100, 
        Contours -> Range[-0.9, 0.9, 0.2], 
        ContourShading -> False, Frame -> False, 
        Epilog -> {RGBColor[0, 0, 1], Line[{{31, 41}, {70, 41}}], 
            RGBColor[1, 0, 0], Line[{{31, 60}, {70, 60}}]}];

Finally, we use PlotGradientField to display the electric field.

PlotGradientField[-f[x, y], {x, 1, 100}, {y, 5, 100}, 
        PlotPoints -> 15, ScaleFactor -> 5, 
        Epilog -> {RGBColor[1, 0, 0], Line[{{31, 41}, {70, 41}}], 
            Line[{{31, 60}, {70, 60}}]}];

Capacitance

According to Gauss' law, specialized to two-dimensions, the total charge per unit transverse length on an electrode is

\(Q\  = \ \(ε\_0\) \(∮ \(E\& ⇀ \)·d\(S\& ⇀ \)\)\)

where (Subscript[ε, 0]) is the electric permittivity, \(TraditionalForm\`E\&⇀ = \(-∇\&⇀ ϕ\)\) is the electric field, \(TraditionalForm\`dS\&⇀ = \(n\&^\) dS\), dS is a differential length for any closed curve that encompasses the electrode (and none other), and \(TraditionalForm\`n\&^\) is the outward normal unit vector for that curve.  Assuming that there are only two electrodes, the capacitance (per unit length) for this system is then

C = Q/V 
\(V\  = \ \(∫\_2\%1 \(E\& ⇀ \)·d\(s\& ⇀ \)\  = \ Δ ϕ\_1 - Δ ϕ\_2\)\)

where V is the potential difference between the two electrodes.

For the present problem, we choose to enclose the upper electrode by a rectangle, as sketched below, and we know from the chosen boundary conditions that the potential difference between the two electrodes is 2 V.  Note that we prefer to use the exact voltage rather than compute the potential difference using a line integral because the discontinuities across the plates cause numerical instabilities in that integration.

PlotGradientField[-f[x, y], {x, 1, 100}, {y, 5, 100}, 
        PlotPoints -> 15, ScaleFactor -> 5, 
        Epilog -> {RGBColor[1, 0, 0], Line[{{31, 41}, {70, 41}}], 
            Line[{{31, 60}, {70, 60}}], GrayLevel[0], Dashing[{0.02, 0.02}], 
            Line[{{20, 50}, {80, 50}, {80, 70}, {20, 70}, {20, 
                  50}}]}];

Ex[x_, y_] = -D[f[x, y], x];
Ey[x_, y_] = -D[f[x, y], y];

Cap =  ((Subscript[ε, 0])/2) (Integrate[
              Ey[x, 70] - Ey[x, 50], {x, 20, 80}] + 
            Integrate[Ex[80, y] - Ex[20, y], {y, 50, 70}])

It is instructive to compare this result with the textbook analysis of an ideal capacitor without fringe fields.  We imagine that the electric field is uniform between the plates and vanishes outside the capacitor.  The voltage is then V = E d, where d is the separation between the plates, and the charge on each plate is Q = (Subscript[ε, 0]) E w, where w is the width of the plate.  Thus, the ideal capacitance is C = (Subscript[ε, 0]) w/d, which takes the value 2  (Subscript[ε, 0]) for the present parameters.  Therefore, the actual capacitance is about 1.8 times the ideal capacitance without edge effects and fringe fields.  With the ratio of width to separation being only 2, this is a relatively open geometry that permits significant fringe fields.  Instead of being distributed uniformly on inside surfaces only, the charge density is largest near the edges and is also present on the outside surfaces, such that the charged storage for given voltage is substantially larger than the simple textbook estimate.  This estimate only becomes reasonable when w/d becomes very large.  Numerical methods are obviously required for many problems!  These issues are explored at greater depth in the exercises.

Comments about speed and efficiency

This method is best suited for computation of the potential and field relatively close to the electrodes where edge effects caused by artificial confinement to a rectangular grid at zero potential are small.  The edge effects can become important near the boundary of the computational grid.  Furthermore, if the artficial box compresses the equipotential surfaces, we may overestimate the fringe fields and, in the case of the parallel-plate capacitor, the capacitance of the system.  If we do require accurate fringe fields, we must expand the grid and endure longer computations.  If we are more interested in far field solutions, we should probably evaluate multipole amplitudes and employ a multipole expansion in that region instead.

Generating an accurate solution to even a simple two-dimensional pde obviously requires a considerable amount of cpu time; imagine how much longer it would take using the direct relaxation method instead of the more efficient simultaneous over-relaxation method.  For problems exhibiting symmetries in their boundary configuration and values, one can often reduce the amount of computation time by performing the evaluations within a suitable subregion only and using symmetry to map that portion of the solution onto the entire region.  Thus, for the parallel-plate capacitor we need only compute one quadrant to obtain the entire solution.  However, this strategy generally requires a specialized program for every problem.  One then needs to determine whether the human time needed to specialize the algorithm is more valuable than the computer time needed to execute the more general but less efficient program.  Unless the computation times become prohibitively long, I tend to value human resources over computer resources and recommend using well-tested general programs instead of custom programs for each problem.  After all, computers are getting faster all the time but humans are not!  Nor does the computer require as many coffee breaks.

Interpreted languages, like Mathematica, are not necessarily the best choice for numerical algorithms that require many iterations — compiled languages, like Fortran or C, will generally perform such calculations much more rapidly.  Although the Mathematica Compile function can often improve the performance of numerical computations, the present implementation of the SOR algorithm is not particularly amenable to Compile.  Perhaps a better approach would be to write the relaxation algorithm in C and to link that program to Mathematica using MathLink.  We could then use Mathematica for display but a faster compiled language for the cpu-intensive computations.

Electrodes outlined as polygons

Suppose that we wish to compute the electrostatic potential outside a system of electrodes, given as closed polygons represented by a list of points in which the first and last coincide.  To initialize the potential array, we need a function that determines a point is enclosed by an electrode.  The query InPolygonQ was developed in a separate notebook and will serve this purpose nicely.

InPolygonQ[pt_, poly_] := 
      Module[{t, cpoly, n, k}, 
 (* 
          Prepare closed figure. *) 
cpoly = 
          If[Last[poly] ≠ First[poly], Append[poly, First[poly]], 
            poly]; 
n = 
          Length[cpoly] - 1; 
If[
 (* 
            Check whether point is on the boundary. *) 
MemberQ[
Table[
              If[0 ≤ (t /. 
                      Flatten[
                        Solve[pt == 
                            t*cpoly[[
                                  k]] + (1 - t)*
                                cpoly[[
                                  k + 1]], 
                          t]]) ≤ 1, True, False, False], {k, 1, 
                n}], True], 
          True, 
 (* 
            Use winding number to test whether the point is inside. *) 
Chop[(1/(2  π)) 
                Sum[Arg [(cpoly[[k + 1]] - pt) . {1. , I}/((cpoly[[k]] - pt)) . {1. , I}], {k, 1, n}]] ≠ 0]
]

example: right triangle

First we construct the boundary.

Clear[boundary, ϕ, ϕ0];
boundary = Table[False, {i, 1, 100}, {j, 1, 100}];
ϕ = Table[0.0, {i, 1, 100}, {j, 1, 100}];

poly = {{41, 41}, {60, 41}, {41, 60}, {41, 
            41}};
Show[Graphics[Line[poly]], 
        AspectRatio -> Automatic];

Then we evaluate and display the potential and field using standard steps.

Do[
        boundary[[i, j]] = 
          InPolygonQ[{i, j}, poly], {i, 1, 100}, {j, 1, 100}];
Do[ϕ[[i, j]] = 
          If[boundary[[i, j]], 1. , 
            0. ], {i, 1, 100}, {j, 1, 100}];
ϕ0 = ϕ;

ϕ = relax[ϕ0, boundary, 1000, (10^(-4))];

f = ListInterpolation[ϕ];

ContourPlot[f[x, y], {x, 20, 80}, {y, 20, 80}, 
        PlotPoints -> 100, AspectRatio -> Automatic, 
        ContourShading -> False, 
        Contours -> Range[0.35, 0.95, 0.1], Frame -> False, 
        Epilog -> {GrayLevel[0.75], Polygon[poly]}];

PlotGradientField[Evaluate[-f[x, y]], {x, 30, 70}, {y, 32, 70}, 
        PlotPoints -> 15, ScaleFactor -> 5, 
        Epilog -> {GrayLevel[0.75], Polygon[poly]}];

DisplayTogether[
ContourPlot[
          f[x, y], {x, 22, 80}, {y, 22, 80}, PlotPoints -> 100, 
          AspectRatio -> Automatic, ContourShading -> False, 
          Contours -> Range[0.35, 0.95, 0.1], Frame -> False, 
          Epilog -> {GrayLevel[0.75], 
              Polygon[poly]}], 
PlotGradientField[
          Evaluate[-f[x, y]], {x, 22, 80}, {y, 22, 80}, 
          PlotPoints -> 15, ScaleFactor -> 5]];

Notice that the innermost equipotentials tend to follow the shape of the electrode,while the outer ones tend to be circular until one approaches the artificial square edge of the computational mesh,which distorts the ideal circular shape in the far field.

Electrode defined by parametric curve

Consider an electrode that can be described by a closed parametric curve of the form 𝒞[t] = {x[t], y[t]} where 𝒞[
          t + 2  π n] = 𝒞[t] for integer n.  We can approximate the potential surrounding this electrode by discreting the boundary and applying our present tools.  For example, the polygon sketched below approximates an electrode in the shape of a lemniscate.  It is defined by 50 points, but looks smooth enough for our purposes.  (The more points, the longer the calculation.)

poly = 
        Table[10.  {Sin[t], Sin[2  t]} + {50.5, 50.5}, {t, 0, 
            2  π, ((2  π)/50)}];
ListPlot[poly, PlotJoined -> True];

The construction of the boundary is now somewhat time-consuming

Clear[boundary, ϕ, ϕ0];
boundary = Table[False, {i, 1, 120}, {j, 1, 120}];
ϕ = Table[0.0, {i, 1, 120}, {j, 1, 120}];

Do[
        boundary[[i, j]] = 
          InPolygonQ[{i, j}, poly], {i, 1, 120}, {j, 1, 120}];

DisplayTogether[ListPlot[Position[boundary, True]], 
        ListPlot[poly, PlotJoined -> True]];

Notice that the result can be inaccurate if the mesh does not align well with the polygon vertex points.  If that occurs, we must use a finer mesh that is better able to describe curved surfaces.  Here the mesh does not reproduce the left and right sides of the electrode particulary well and is somewhat ragged elsewhere.  However, we will not spend the time now that would be required for a better approximation to this electrode.

We now compute the potential and field using familiar steps.

Do[ϕ[[i, j]] = 
          If[boundary[[i, j]], 1. , 
            0. ], {i, 1, 120}, {j, 1, 120}];
ϕ0 = ϕ;
ϕ = 
        relax[ϕ0, boundary, 1000, 
          (10^(-4))];
f = ListInterpolation[ϕ];

ContourPlot[f[x, y], {x, 20, 80}, {y, 20, 80}, 
        PlotPoints -> 200, AspectRatio -> Automatic, 
        ContourShading -> False, 
        Contours -> Range[0.49, 0.99, 0.1], Frame -> False, 
        Epilog -> {GrayLevel[0.75], Polygon[poly]}];

DisplayTogether[
ContourPlot[
          f[x, y], {x, 22, 80}, {y, 20, 80}, PlotPoints -> 200, 
          AspectRatio -> Automatic, ContourShading -> False, 
          Contours -> Range[0.49, 0.99, 0.1], Frame -> False, 
          Epilog -> {GrayLevel[0.75], 
              Polygon[poly]}], 
PlotGradientField[
          Evaluate[-f[x, y]], {x, 22, 80}, {y, 20, 80}, 
          PlotPoints -> 15, ScaleFactor -> 8]];

The large spacing between equipotentials within the neck region shows that the electric field is especially small there because repulsion of like charges produces a relatively small surface charge density near the neck.

Interior shielding by Faraday pail (thick U)

The thick U-shaped electrode shown below is sometimes called the Faraday pail because it tends to expel electric fields from the interior of the pail even though it is a fairly open geometry.  This behavior is illustrated by placing a circular electrode over the mouth of the pail.  First, we construct and display the electrodes.

poly1 = {{49, 60}, {49, 40}, {71, 40}, {71, 60}, {68, 60}, {68, 
            43}, {52, 43}, {52, 60}, {49, 60}};
poly2 = 
        Table[4.  {Cos[t], Sin[t]} + {60, 70}, {t, 0, 
            2  π, ((2  π)/100)}];
Show[Graphics[{Line[poly1], Line[poly2]}], 
        AspectRatio -> Automatic, PlotRange -> All];

Then, we evaluate the potential and field using familiar steps.  The circle is charged to 1 V while the pail is at ground (as is the edge of the computational region).

Clear[boundary, ϕ, ϕ0];
boundary = Table[False, {i, 1, 120}, {j, 1, 120}];
ϕ = Table[0.0, {i, 1, 120}, {j, 1, 120}];

Do[
        boundary[[i, j]] = 
          InPolygonQ[{i, j}, poly1], {i, 1, 120}, {j, 1, 
          65}];
Do[
        boundary[[i, j]] = 
          InPolygonQ[{i, j}, poly2], {i, 1, 120}, {j, 65, 
          120}];
Do[ϕ[[i, j]] = 
          If[boundary[[i, j]], 0. , 
            0. ], {i, 1, 120}, {j, 1, 65}];
Do[ϕ[[i, j]] = 
          If[boundary[[i, j]], 1. , 
            0. ], {i, 1, 120}, {j, 65, 120}];
ϕ0 = ϕ;

ϕ = 
        relax[ϕ0, boundary, 1000, 
          (10^(-4))];
f = ListInterpolation[ϕ];

ContourPlot[f[x, y], {x, 1, 120}, {y, 1, 120}, 
        PlotPoints -> 200, PlotRange -> All, 
        AspectRatio -> Automatic, ContourShading -> False, 
        Contours -> Range[0.05, 1. , 0.1], Frame -> False, 
        Epilog -> {RGBColor[1, 0, 0], Polygon[poly1],  
            RGBColor[0, 0, 1], Polygon[poly2]}];

Notice that potential is pretty uniform and the electric field is quite small in the bottom of the pail, despite the charged object above it.

PlotGradientField[-f[x, y], {x, 42, 80}, {y, 30, 83}, 
        PlotPoints -> 15, ScaleFactor -> 4, 
        Epilog -> {RGBColor[1, 0, 0], Polygon[poly1],  
            RGBColor[0, 0, 1], Polygon[poly2]}];

\(TraditionalForm\`∂\^2 ϕ\/∂x\^2 + ∂\^2 ϕ\/∂y\^2 == \(-ρ\)\)

\(TraditionalForm\`∂\^2 ϕ\/∂r\^2 + \(1\/r\) ∂ϕ\/∂r + ∂\^2 ϕ\/∂z\^2 == 0\)