MathSource Reviews  7(4)

Solving the Laplace Equation

Matthew M. Thomas

In the previous issue of this journal, we saw how the heat equation in one dimension found a noteworthy application in economics, therein providing a bridge from the physical sciences to their dismal counterpart, and earning Nobel Laureate status for bridge-builders Scholes, Merton, and (posthumously) Black. In this issue, we focus our attention on the physical sciences side of that bridge. We also broaden our perspective beyond a single dimension, and we constrain ourselves to steady state. Our focus is on the 2-D Laplace equation, which is of course the 2-D heat equation at steady state (although its physical meaning goes well beyond heat transfer). We begin by pondering Laplace and his contemporaries, contemplating their accomplishments in the turbulent wake of their time and place in history. We then narrow our scope to the Laplace equation itself, though our narrowed scope blinds us not to the broad physical meaning previously implied. From there, our attention turns to Gennady Stupakov and his contribution to MathSource: A Laplace equation solver for electrostatic problems, given Dirichlet boundary conditions. The Stupakov MathSource items comprise a Mathematica package and a version 3.0 notebook; we will examine both. Finally, our summary will view the strengths and weaknesses of the Stupakov contribution, from what might be described as different sets of eyes.

Pierre-Simon Laplace

Arguably, the worst place for one of professional rank in the 20th century was Phnom Penh, Cambodia in the mid- to late-1970s. In establishing whatever it was that they sought to establish, Pol Pot and his Khmer Rouge chose a totalitarian agrarian existence for nascent Democratic Kampuchea. Alas, those of professional rank were deemed unworthy of the collective. Ever the innovators, the Khmer Rouge used the professionals not for solving the problems they were trained to solve, but for lining the Killing Fields with their skeletal remains. Bone supersedes brain -- 'twas the Khmer approach. To say that Paris was to the early 1790s what Phnom Penh was to the late 1970s is an overstatement, in terms of unfettered butchery. Nevertheless, overactive Jacobin guillotines under the Committee of Public Safety were hardly positive reinforcement for prominent French savants of the day. And what a prominent lot they were: first, Lagrange, Lavoisier, Laplace, and Legendre; later, Fourier, Biot, and Poisson. On the whole, they not only avoided the blade (save for Lavoisier), but also made contributions to their trades that are withstanding the test of time. We cannot know what would have sprung forth from the erstwhile French protectorate in southeast Asia two decades ago; we do know and continue to benefit from what the motherland produced two centuries ago. 

Born in 1749, Pierre-Simon Laplace got his start as a mathematics professor at the Military School of Paris, thanks to the recommendation of Jean le Rond d'Alembert (he of wave equation and mechanics fame). His first feat was to prove the stability of the (idealized) solar system in 1773; the proof won him entry into the Academy of Sciences. In 1780, he and Antoine-Laurent Lavoisier (that "Father of Modern Chemistry" whose tax collecting won him a date with the guillotine) collaborated on an ice calorimeter invention. For the most part, his work on the equation that now bears his name took place in the 1780s: More on that and the potential function in due course. In 1785, he both won full membership in the Academy and examined a noteworthy young candidate -- Napoleon Bonaparte -- for admission to the Military School. During the Revolution, both Laplace and Lagrange literally saved their necks through military service, by calculating artillery trajectories and supervising gunpowder manufacturing. Afterwards, Laplace took the political route, serving first Napoleon and later the reinstituted monarchy with an equanimity that would swell Talleyrand with pride. After additional noteworthy efforts in probability and celestial mechanics, Laplace died in 1827. He lives on, through his accomplishments, and through biographical essays [Bell 1937] and books [Gillispie 1997] on his life.

Whereas Joseph-Louis Lagrange came to be regarded as a pure mathematician, Laplace was seen as a physicist for whom mathematics was but a means to an end. Since Laplace undertook the greater challenges in celestial mechanics, the upstart Siméon-Denis Poisson was said to hold Laplace in slightly higher esteem than he held Lagrange. As a mathematical physicist himself, fellow upstart Jean-Baptiste-Joseph Fourier also appreciated Laplace's practical perspective, yet gave due credit to Lagrange for the depth and accuracy of Lagrange's derivations. (A Laplace derivation, by contrast, often contained the wild card "Il est aisé à voir" ["It is easy to see;" Bell 1937, pp. 177-178] -- surely a precursor to the "proof shall be left as an exercise for the student" phrase that blights modern textbooks). Two centuries later, ranking the Revolutionary-era French savants seems frivolous: Each -- be it Lagrange for his method of variation of parameters, Lavoisier for his quantitative chemistry, Laplace for what is to come, Adrien-Marie Legendre for his ODE and polynomials, Fourier for the heat equation (presaged and championed by Laplace [Gillispie 1997, p. 251; Grattan-Guinness and Ravetz 1972, pp. 444-452]), Jean-Baptiste Biot for relating magnetic induction to electric current, Poisson for extending the Laplace equation, and all again for works not cited here -- is highly regarded. Herein, we narrow our focus to Laplace.

The Laplace Equation; Potential

Equation (1) is the Laplace equation in two dimensions.

(∇^2) u  =  
      (∂^2) (u/∂)(x^2) +  
          (∂^2) (u/∂)(y^2) = 
        0                                                                                                                                                          (1)

Of course, this second-order elliptic partial differential equation bears a strong similarity to the 2-D heat equation. Change the right-hand side of equation (1) to the partial derivative of u with respect to time t, and you have the dimensionless heat equation. But to say that equation (1) followed from the heat equation is to state an anachronism. Credit for the heat equation and its solution lies with Fourier's efforts in the early 1800s. Laplace brought forth equation (1) in the 1780s, not in a heat transfer context, but in an analysis of Saturn's rings [Laplace, 1895, p. 278 -- originally published in 1789]. Equation (1) is a wonderful tool for describing steady-state thermal gradients. But to constrain this equation to a heat transfer context is to ignore the larger scope of its application. Let us thus consider the larger scope of its application.

Laplace in the 1780s was motivated by Newton's work in gravitation. His results, applied initially to celestial mechanics (and applicable to the aforementioned steady-state heat transfer), were a major contributor to the rise of potential theory as an independent branch of mathematics. Potential theory has applications in electrostatics, magnetostatics, electrodynamics, complex analysis, conformal mappings, and differential geometry. It is in an electrostatics context that Stupakov presents the Laplace Equation Solver MathSource item. In equation (1), dependent variable u[x,y] is labeled the potential function or the harmonic function; while the latter is perfectly legitimate, we will use the former. Advanced mathematics textbooks [e.g. Kreyszig 1979, pp. 551-557; Spiegel 1981, pp. 566-567] offer a basic treatment on the concept of the potential. A colorful treatment [Gellert et al 1989, pp. 693-698] of potential theory assigns equation (1) to Laplace circa 1782, the term "potential" to Gauss circa 1840, and the linkage between the partial derivatives of u and the components of gravitational or electrostatic forces to Lagrange circa 1773.

In the wake of this treatment, we turn to E. T. Bell: "To call a spade a spade, Laplace stole outrageously, right and left, wherever he could lay his hands on anything of his contemporaries and predecessors which he could use. From Lagrange, for example, he lifted the fundamental concept of the potential." Bell also scolds Laplace for taking "whatever he needed in the way of analysis" from Legendre, and for citing no one but Newton in Laplace's works on celestial mechanics [Bell 1937, p. 174]. Perhaps word was out among pre-Revolutionary aristocrats to count the silverware after a visit from Pierre-Simon Laplace.

No discussion of equation (1) -- which uncontestedly belongs to Laplace -- is complete without also discussing the boundary conditions that are in effect. Within the domain that equation (1) applies, assume an open region of three-dimensional space bounded by a smooth surface that is either open or closed. Let n denote a spatial coordinate in the direction normal to the surface. There then are three types of boundary conditions that can apply. Dirichlet conditions specify potential u on the boundary (often, u = 0). Neumann conditions specify the partial derivative of u with respect to n on the boundary. Cauchy or mixed conditions specify a linear combination of u and its partial derivative on the boundary. For elliptic PDEs such as equation (1), Dirichlet conditions will lead to a unique, stable solution that depends continuously on the boundary conditions when the boundary is closed. A useful table of PDE solution types for boundary condition type vs. open or closed boundary vs. type of PDE (elliptic, hyperbolic, or parabolic) is available [Zwillinger 1996, pp. 417-418]. These boundary conditions in conjunction with Green's formula have also been discussed [Gellert et al 1989, pp. 696-697]. In the following treatment, we limit the discussion to a 2-D electrostatic system with closed boundaries and Dirichlet conditions in effect.

Laplace Equation Solver: package

MathSource directory www.mathsource.com/Content/Applications/Physics/0208-460 contains two items: package LESolver.m (6 kB) and version 3.0 notebook Documentation.nb (1.5 MB). Both are the handiwork of Gennady Stupakov, who had been at the Stanford Linear Accelerator Center at least through 1996. The version of LESolver.m in MathSource as of late September 1998 was version 1.00, dated November 1996. Two email addresses for Stupakov are listed in the package and the notebook: stupakov@wco.com, and stupakov@slac.stanford.edu. The Laplace Equation Solver abstract states that its Mathematica code solves the 2-D Laplace equation for a closed boundary described by Dirichlet conditions.

If you are running version 3.0 of Mathematica, you can use LESolver.m to characterize electrostatic potential both inside and outside the boundary. If you are running an earlier version, you can perform an interior boundary characterization with the package as is, but you will need a special MathLink program for exterior characterization. The reason: An enhanced version 3.0 SingularValues[] function, which is better than its 2.x counterpart, which was far better than its disastrous 1.x counterpart (which upon decomposition belched forth singular values or positive eigenvalue square roots not necessarily sorted in descending order -- this to the extreme consternation of principal component analysts, partial least squarers, and other practitioners of continuum regression).

The LESolver.m package itself does not contain much by way of documentation, but the code is clean and the appropriate usage statements are intact. The package comprises four primary functions: Prepare, DrawPlot, Solver, and Potential. Prepare[] uses Map[] and Table[] to perform triage, preparing a mesh and lists for use by the other functions. Given the Prepare[] output, DrawPlot[] illustrates that mesh, as well as the boundary (and its normal vectors) supplied by the user. Solver[] generates an integral equation solution on the interior or the exterior of the boundary; it also makes use of the Prepare[] output. Potential[] uses the Prepare[] output, the Solver[] output, and user-defined x- and y-coordinate parametric equations to compute electrostatic potential at a given x,y coordinate. (The parametric equations define the boundary.) As the notebook shows, ListPlot3D[] and ListContourPlot[] are used to illustrate Potential[] output. Il est aisé à voir that LESolver.m is making use of a numerical method [e.g. Jeffrey 1996, pp. 792-801] to characterize electrostatic potential; alas, there is no documentation to divulge details.

Notice: As of late September 1998, there were four drastic flaws in the version of LESolver.m found in MathSource. The four flaws were four misplaced carriage returns (CRs) within the Potential[] part of the package. The extraneous CRs split variable names ttab1 in line 173, xxmin in line 196, and yymin in line 198; another CR splits function Table[] in line 206. These CRs cause Mathematica to treat the split term portions as multiplicand and multiplier: ttab1 becomes t * tab1, xxmin becomes x * xmin, yymin becomes y * ymin, and Table becomes Tab * le to the confused Mathematica kernel. Upon your loading the flawed LESolver.m version, Mathematica will generate four Syntax::newl messages to warn you of the flaws. But if you delete the extraneous CRs prior to loading LESolver.m, the package will perform as intended.

Laplace Equation Solver: notebook

While the package contains relatively few comments, the Documentation.nb notebook does much to compensate for that deficiency. This notebook explains proper package use by example, first through a case whose boundary is trivial, then through a case whose boundary is non-trivial. In the trivial case, the boundary is a unit circle, centered on the coordinate-system origin; the electrostatic potential characterization is interior to the closed boundary. In the non-trivial case, the boundary is a slab with rounded ends, again centered on the coordinate-system origin; the electrostatic potential characterization is exterior to the closed boundary. In each case, the user defines the x- and y-coordinate parametric equations corresponding to the closed boundary, and defines boundpoten[], a parametric equation describing the Dirichlet boundary condition. The user then employs Prepare[], DrawPlot[], Solver[], and Potential[] in conjunction with Mathematica's list-plotting capabilities to illustrate electrostatic potential for the given boundary and Dirichlet condition.

In the trivial case, the LESolver.m solution can be compared to a well-established analytical solution; while the LESolver.m solution will push no state-of-the-art, a high degree of similarity between it and its analytical cousin would reassure users of LESolver.m. For this problem, the boundary is an origin-centered unit circle, and we will characterize the electrostatic potential difference within the boundary interior. We first define the x- and y-coordinate parametric equations xx[s_] := Cos[s] and yy[s_] := Sin[s], and Dirichlet boundary condition boundpoten[s_] := Log[(0.5-xx[s])^2 + yy[s]^2]. These three equations are sent to Prepare[], with n = 40 specifying 40 mesh points on the boundary; the output then goes to DrawPlot[], which produces Figure 1. The output also goes to Solver[], which solves the 2-D Laplace equation for the boundary interior (based upon a Domain -> Interior option). Solver[] and Prepare[] outputs then become inputs to Potential[], which computes electrostatic potential at point x,y. A list ("potenarray") is then generated to hold, at each x,y point inside the boundary, the difference between Potential[] output and the boundpoten[] value. Figures 2 and 3 offer 3-D and contour perspectives, respectively, of those list values, therein indicating the electrostatic potential difference at any point within the boundary.

For this particular problem, there is an analytical solution to the 2-D Laplace equation with a Dirichlet boundary condition [Zwillinger 1996, p. 420; Sokolnikoff and Redheffer 1966, pp. 474-476]. The analytical solution, appropriately enough, is an integral named for Laplace protégé Siméon-Denis Poisson; equation (2) provides that integral:

u[r, θ] =(1/(2  π))(∫_0^(2  π))((1 - (r^2))/(1 - 2 r Cos[θ - ϕ] +  (r^2)))f[ϕ]\[DifferentialD]ϕ                                                       (2)

In equation (2), f[ϕ] = u[1,θ] for 0 <= θ < 2 π, which is the Dirichlet boundary condition defined in boundpoten[]. The Documentation.nb notebook compares Potential[] output to equation (2) output for specific x,y points within the boundary. The resulting electrostatic potentials are indeed acceptably similar (i.e., within 1% of each other). 

The approach to the non-trivial case is identical to the trivial case approach. For the non-trivial case, however, the x- and y-coordinate parametric equations are more complicated. (They have to be; they are now representing not a circle but a straight slab with rounded ends. The Which[] function sees use in each.) The Dirichlet boundary condition in this case is boundpoten[t_] := -xx[t]. The calls to Prepare[] and DrawPlot[] are as before, but with double the number of mesh points (n = 80) on the boundary. Figure 4 is the DrawPlot[] output; compare it to Figure 1. For this case, Solver[] is called using the Domain -> Exterior option. The balance of the non-trivial exterior case procedure is identical to its trivial interior counterpart, and Figures 5 and 6 are the equivalents of Figures 2 and 3 (respectively). Stupakov does well to present this case in Documentation.nb: It illustrates to the user how a non-trivial boundary must be represented, both by the x- and y-coordinate parametric equations and within the "potenarray" list. The Exterior option worked well in Mathematica version 3.0; your reviewer did not assess its performance in version 2.x.

Notice: As of late September 1998, there was also a drastic flaw in the version of Documentation.nb found in MathSource. The flaw was an undefined function "g" within a call to Prepare[] in the non-trivial example. The notebook offers the graphics that would result were the flaw not there; however, it is indeed a hapless user who tries to generate the slab boundary using the Prepare[] function input as is. Fortunately, this flaw is also easy to fix. Just above the flawed function input is a definition of boundpoten[] ... as a function of the x-coordinate parametric equation for this case. Replace "g" with "boundpoten" in the flawed input, and the Prepare[] and DrawPlot[] functions will duly create the slab boundary you see in both Documentation.nb and Figure 4. One suspects that the problematic "g" function is a remnant of some draft version of this Mathematica notebook -- a stowaway on the trip to MathSource. You now know where the stowaway hides; cast him overboard and seat the rightful passenger.

Summary

If we were to summon the French savants from two centuries ago, and let them work over Gennady Stupakov's MathSource package and notebook, what (given what we think we know about the savants) might be the result? Conjecture: The meticulous Lagrange ferrets out the extraneous carriage returns in LESolver.m and the undefined "g" function in Documentation.nb. The levying Lavoisier contemplates the duty were these MathSource items not free. The insouciant Laplace thinks nothing of the numerical solution method documentation absence from LESolver.m; it is, after all, easy to see. The pensive Legendre ponders an extension of LESolver.m to spherical Laplace equation solutions, so that the eponymous developer of the Legendre ODE might make use of his own equation. The dynamic Fourier feels ill at ease and bored without a non-zero partial derivative (with respect to time) on the right-hand size of Laplace's equation. The slighted Biot wants an explicit magnetostatics example worked into Documentation.nb. The contented Poisson, his own integral at work in verifying the methods of LESolver.m in Documentation.nb, finds the latter to be just fine as is ... Monsieur Biot notwithstanding. And your reviewer, standing well clear of the Frenchmen, who have after all been dead for lo these many years, applauds Gennady Stupakov for an imperfect but nevertheless useful package on 2-D Laplace equation solutions under Dirichlet boundary conditions.

References

Bell, E. T. "From Peasant to Snob," Men of Mathematics, Chapter 11. Simon and Schuster, New York (1937).

Gellert, W., Gottwald, S., Hellwich, M., Kästner, H. and H. Küstner, editors. VNR Concise Encyclopedia of Mathematics. Second Edition. Van Nostrand Reinhold, New York (1989).

Gillispie, Charles Coulston. Pierre-Simon Laplace -- 1749-1827: A Life in Exact Science. Princeton University Press, Princeton, NJ (1997).

Grattan-Guinness, I. in collaboration with J. R. Ravetz. Joseph Fourier -- 1768-1830. MIT Press, Cambridge, MA (1972).

Jeffrey, Alan. Mathematics for Engineers and Scientists. Fifth Edition. Chapman & Hall, London (1996).

Kreyszig, Erwin. Advanced Engineering Mathematics. Fourth Edition. John Wiley and Sons, New York (1979).

Laplace, Pierre-Simon. Oeuvres Complètes de Laplace. Tome Onzième. Gauthier-Villars et Fils, Paris (1895).

Sokolnikoff, I. S. and R. M. Redheffer. Mathematics of Physics and Modern Engineering. Second Edition. McGraw-Hill, New York (1966).

Spiegel, Murray R. Applied Differential Equations. Third Edition. Prentice-Hall, Englewood Cliffs, NJ (1981).

Zwillinger, Daniel, editor. CRC Standard Mathematical Tables and Formulae. Thirtieth Edition. CRC Press, Boca Raton, FL (1996).

About the Author

Matthew M. Thomas completed his doctorate in chemical engineering at Washington University in St. Louis, in December 1995. He and compatriots Thad Podgajny and Steve Buckner are the co-inventors responsible for U.S. Patent 5,712,613, issued on 27 January 1998; the patent features Mathematica in a product that reduces electromagnetic wave scatter.

Figures

Figure 1 - Origin-centered unit circle boundary, with 40 mesh points on the boundary and a normal vector illustrated at each.

Figure 2 - Electrostatic potential difference superimposed upon the closed boundary from Figure 1, as generated by ListPlot3D[] and ParametricPlot3D[], respectively.

Figure 3 - Contour plot representation of the electrostatic potential difference and closed boundary from Figures 1 and 2; ListContourPlot[] supplants ListPlot3D[] and ParametricPlot[] supplants ParametricPlot3D[] in the generation of this representation.

Figure 4 - Boundary for a straight slab with rounded ends, with 80 boundary mesh points and a normal vector at each.

Figure 5 - Electrostatic potential difference in 3-D, superimposed upon the closed boundary from Figure 4. Note how, in this case, that difference is zero within the boundary interior.

Figure 6 - Contour plot representation of the electrostatic potential difference and closed boundary from Figures 4 and 5. Again: The potential difference is zero within the boundary interior.