# Poisson solver using Mathematica

I am looking for some help with a Poisson solver I am writing in Mathematica. The code is quite long with Arrays plugged in, so the full details can be found at http://pastebin.com/uSrSDcW6

I am calculating voltages given charge densities using the central difference method derived from Poisson's equation. After calculating the voltage, I test the data set for convergence. I am setting convergence thresholds on the order of 10^-1000+. I have the loop set up to kick out after 10000 iterations in case something goes awry, as a fail safe. I have a loop counter in place for sanity. The program seems to run fine as long as the convergence threshold is set to 10^-10.

My question is this: No matter what I update the threshold to, for example, 10^-100 or 10^-150, the computation stops after 633 iterations and kicks out of the loop. I would appreciate any help with this, I am completely stuck. I've added comments to the program that should be explanatory for anyone on this forum.

*Update 10/9/12** I've narrowed my problem down to Machine Precision. I need to expand the Machine Precision to my Max Machine Precision of ~10^-309. How do I set this globally within my program when Im using 10x10x34 arrays initialized at ~15 digit precision?

• @Nasser do you have any code allowing one to solve the Poisson pde in 3D, in particular in spherical coordinates? I was looking through the demonstrations that you posted online but can only find solutions in one or two dimensions. – Vincent Tjeng Apr 20 '13 at 4:22

## 1 Answer

I can't find the actual code in your linked data file, but it may be worth posting my own solution for a 2D Poisson problem here. It is copied from my web page. I'm using a maximum of 100000 iterations by default. From your description, it sounds as if you could try to re-write your loops using constructs such as Fold, Nest or - as I do below - FixedPoint. The latter is a natural choice for relaxation methods.

I don't have much time right now, but I hope this will provide some ideas to you since my code is relatively short and contains all the logic plus all the data.

The code (main program is poissonSolver) allows you to specify charge densities as well as conductors at fixed voltages. Moreover, there is also the possibility to specify an arbitrary spatially varying electric susceptibility.

My central focus in this solution was to find a good way to enter all the spatially varying functions in a graphical way. This means that I didn't focus on numerical refinements. But the speed of the solver is perfectly fine for practical applications.

The input of the program is actually a set of Graphics objects, which are then assembled in a routine createLandscape. The advantage is that you can create those graphics representing your conductors, dielectrics and charges in any graphics program of your choice, and paste them right into the notebook.

Here, I'm drawing everything within the notebook to remain self-contained.

All the ingredients are supplied separately as 2D Graphics objects:

• conductors (Black = no conductor, all other GrayLevel values give the value of the fixed potential). The potential on each conductor should be constant. Because I use GrayLevel to encode the potential, all potentials are positive (as can always be arranged with a suitable additive constant).
• chargeMinus (Black = zero charge, White = large negative charge)
• chargePlus (Black = zero charge, White = large positive charge)
• susceptibility (susceptibility; Black = 0)

The code follows:

defaultGridSize = 120;

step = Compile[{{phi, _Real, 2}, {orig, _Real, 2}, {mask, _Real,
2}, {dchiX, _Real, 2}, {dchiY, _Real, 2}},
Module[{f = (RotateRight[phi] + RotateLeft[phi] +
RotateRight[phi, {0, 1}] + RotateLeft[phi, {0, 1}])*mask +
orig, ex, ey},
ex = (RotateRight[phi, {0, 1}] - RotateLeft[phi, {0, 1}]);
ey = (RotateRight[phi] - RotateLeft[phi]);
f - (dchiX*ex + dchiY*ey)*mask]];

iterate =
Compile[{{gridArray, _Real, 2}, {originalArray, _Real,
2}, {maskArray, _Real, 2}, {dchiX, _Real, 2}, {dchiY, _Real,
2}, {tol, _Real}},
FixedPoint[step[#, gridArray, maskArray, dchiX, dchiY] &,
originalArray, 100000,
SameTest -> (Max@Abs@Flatten[#1 - #2] < tol &)]];

digitize[gr_, n_] :=
N@ImageData@
ColorConvert[
Image[Show[gr, Background -> Black,
BaseStyle -> {Antialiasing -> False}], ImageSize -> n],
"GrayScale"];

createLandscape[conductors_, chargePlus_, chargeMinus_,
suceptibility_, nGrid_] :=
Module[{gridConductors, gridRho, gridChi, maskList},
gridConductors = digitize[conductors, nGrid];
maskList = N[1. - Unitize[gridConductors]];
gridRho = (digitize[chargePlus, nGrid] -
digitize[chargeMinus, nGrid])*maskList;
gridChi = digitize[suceptibility, nGrid]*maskList;
{gridConductors, gridRho, gridChi, maskList}];

poissonSolver[conductors_, chargePlus_, chargeMinus_, suceptibility_,
nGrid_: defaultGridSize, tolerance_: 10^(-6)] :=
Block[{averagePotential, gridConductors, gridRho, gridChi, gridEps,
gridList, maskList, mask4List, dChiYList, dChiXList,
initialGrid}, {gridConductors, gridRho, gridChi, maskList} =
createLandscape[conductors, chargePlus, chargeMinus,
suceptibility, nGrid];
averagePotential = Mean[Select[Flatten@gridConductors, Positive]];
initialGrid = averagePotential*maskList + gridConductors;
gridEps = 1. + gridChi;
gridList = gridConductors + gridRho/(4.*gridEps);
mask4List = maskList/4.;
dChiYList = (RotateLeft[gridChi] -
RotateRight[gridChi])/(2. gridEps);
dChiXList = (RotateLeft[gridChi, {0, 1}] -
RotateRight[gridChi, {0, 1}])/(2. gridEps);
Reverse@
iterate[gridList, initialGrid, mask4List, dChiXList, dChiYList,
tolerance]];

conductors =
Graphics[{{GrayLevel[.5],
Rectangle[1.1 {-1, -1}, 1.1 {1, 1}]}, {Black,
Disk[{0, 0}, 1]}, {GrayLevel[.1],
Disk[{-.25, 0}, .15]}, {GrayLevel, Disk[{.1, -.25}, .1]}},
PlotRangePadding -> 0, ImagePadding -> None];

chargePlus =
Graphics[{{GrayLevel[.1], Disk[{-.5, .3}, .05]}},
PlotRange -> (PlotRange /. FullOptions[conductors])];

chargeMinus =
Graphics[{{GrayLevel[.1], Disk[{0., .4}, .05]}},
PlotRange -> (PlotRange /. FullOptions[conductors])];

susceptibility =
Graphics[{GrayLevel[.4], Rectangle[{0.3, -1}, {1, 1}]},
PlotRange -> (PlotRange /. FullOptions[conductors])];

Timing[
potential =
poissonSolver[conductors, chargePlus, chargeMinus, susceptibility];]

(* ==> {2.95401, Null} *)

ListPlot3D[potential, PlotRange -> All,
PlotStyle -> {Orange, Specularity[White, 10]}] In this plot, the circular boundary (at constant potential) is visible as a slight crease. The charges create the two spikes (one positive and one negative), and the conductors create flat plateaus. You can also identify a faint, straight-line kink in the potential where I've added a dielectric interface (it runs parallel to the 6-th grid line from the right).

• Thanks for the great comments, and what you have with the graphical representation is where I would like to get once I get my program computing. My issue is an exceedence of the machine precision at 16 digits. Mathematica isn't very helpful when it comes to increasing the machine precision. It bascially says "N[MachinePrecision, 50]". What do I need to do so that precision is applied to my entire Mathematica program? I need to test convergence on the order of <10^-100 and my machine max precision is ~10^300, so I got the computing power to do it. See pastebin.com/uSrSDcW6 – Joel D Oct 9 '12 at 23:46