This section presents the basic properties of solutions to the Laplace and Helmholtz equations.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part VI of the course APMA0340
Introduction to Linear Algebra with Mathematica

Laplace's equation (also called the potential equation or harmonic equation) is a second-order partial differential equation named after Pierre-Simon Laplace who, beginning in 1782, studied its properties while investigating the gravitational attraction of arbitrary bodies in space. However, the equation first appeared in 1752 in a paper by Euler on hydrodynamics. Laplace's equation is often written as:

\[ \Delta u ({\bf x}) =0 , \qquad {\bf x} \in \Omega \subset \mathbb{R}^n , \]
where \( \Delta = \nabla^2 = \frac{\partial^2}{\partial x_1^2} + \frac{\partial^2}{\partial x_2^2} + \cdots + \frac{\partial^2}{\partial x_n^2} \) is the Laplace operator or Laplacian in n dimensional space. Here Ω is some region in ℝn with smooth boundary. In two dimensions, the Laplace equation in rectangular coordinates becomes
\[ \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} =0 \qquad\mbox{or for short} \qquad u_{xx} + u_{yy} =0 . \]
Correspondingly, in three dimensions, Laplace's equation is
\[ u_{xx} + u_{yy} + u_{zz} = 0 , \]
Any real-valued function that satisfies Laplace's equation ∇² u = 0 in some domain Ω is referred to as being harmonic or potential function in Ω.

Another very important version of the above equation is the so called the Helmholtz equation

\[ \Delta u ({\bf x}) + \omega^2 u =0 \qquad \mbox{in domain } {\bf x} \in \Omega \subset \mathbb{R}^n \qquad\mbox{or} \quad {\bf x} \in \mathbb{R}^n \setminus \Omega , \]
named after the German physician and physicist Hermann von Helmholtz (1821--1894). If the region Ω extends to infinity, then it is assumed that a solution to the Laplace's or Helmholtz's equation satisfies additional constraint: it should approach zero at infinity for dimensions larger than 2. In two dimensions, it is sufficient to require soundness of the solution.

A nonhomogeneous version of the Laplace equation

\[ \Delta u ({\bf x}) = f({\bf x}) \qquad \mbox{in domain } \Omega \subset \mathbb{R}^n \qquad\mbox{or} \quad {\bf x} \in \mathbb{R}^n \setminus \Omega , \]
where f is a given smooth function, is called the Poisson's equation. The equation is named after the French mathematician, geometer, and physicist Siméon Denis Poisson (1781--1840). The Poisson's equation in three dimensional case is
\[ u_{xx} + u_{yy} + u_{zz} = f(x,y,z) , \]
where f is a given smooth function. A particular solution of the Poisson equation \( \Delta u = f \) is given by
\[ u(x,y) = \frac{1}{2\pi}\,\iint_{\Omega} f(\xi , \eta )\, \ln r\,{\text d}S , \qquad r^2 = (x-\xi )^2 + (y- \eta )^2 , \]
in two dimensions. For three dimensional space, we have a particular solution to the Poisson equation ∇²u = f:
\[ u(x,y,z) = \frac{-1}{4\pi}\,\iiint_{\Omega} f(\xi , \eta , \zeta )\, \frac{1}{r}\,{\text d}\xi\, {\text d}\eta \, {\text d}\zeta , \qquad r^2 = (x-\xi )^2 + (y- \eta )^2 + (z-\zeta )^2 . \]
In general n-dimensional case, the following function
\[ u({\bf x}) = \frac{-1}{\omega_n}\,\int_{\| {\bf x} - {\bf y} \| \le 1} f({\bf y})\, r^{2-n} ({\bf x}, {\bf y})\, {\text d}V_y , \qquad n \ge 3, \]
gives a particular solution to the Poisson's equation for the Laplacian, where ωn is the surface area of the unit sphere \( x_1^2 + x_2^2 + \cdots + x_n^2 =1 . \)

The Laplacian occurs in different situations that describe many physical phenomena, such as electric and gravitational potentials, the diffusion equation for heat and fluid flow, wave propagation, and quantum mechanics. The Laplacian represents the flux density of the gradient flow of a function. For instance, the net rate at which a chemical dissolved in a fluid moves toward or away from some point is proportional to the Laplacian of the chemical concentration at that point. We summarize the standard examples of its applications in the following table.

   Application       u       f
Incompressible irrotational flow in two or three dimensions u = velocity potential, velocity = q = ∇u Mass source strength/unit volume
Electrostatics u = electrostatic potential; electric field = E = ∇u Charge/unit volume
Steady state temperature in a solid u = temperature Heat source strength/unit volume
Two dimensional incompressible steady flow u = stream function = ψ = velocity = q = ψxi - ψyj -Vorticity
Static deflection of a thin membrane in two dimensions u = deflection Pressure
Newtonian gravitation u = gravitational potential; force of gravity = F = - ∇u Mass density

Laplace's or Helmholtz equations have infinite many solutions. To take out one of them, we need to impose some conditions that usually follow from physical counterpart. Since there is no time dependence in the Laplace's equation or Poisson's equation, there is no initial conditions to be satisfied by their solutions. However, there should be certain boundary conditions on the boundary curve or surface \( \partial\Omega \) of the region Ω in which the differential equation is to be solved. Typically, there are known three types of boundary conditions. The problem of finding a solution of Laplace's equation that takes on given boundary values is known as a Dirichlet problem. However, if the values of the normal derivative are prescribed on the boundary, the problem is said to be a Neumann boundary value problem. Physically, it is plausible to expect that three types of boundary conditions will be sufficient to determine the solution completely. Indeed, it is possible to establish the existence and uniqueness of the solution of Laplace's (and Poisson's) equation under the first and third type boundary conditions, provided that the boundary \( \partial\Omega \) of the domain Ω is smooth (have no corner or edge). However, we avoid explicit statements and their proofs because the material is beyond the scope of the tutorial. Instead, we concentrate our attention on some typical problems that could be solved by means of separation of variables.

Actually, there are known two kinds of boundary value problems: internal or inner when we are after solutions inside the given domain Ω, and external or outer when solutions are to be found outside the given domain, namely in \( {\mathbb R}^n \setminus \Omega . \)

Uniqueness Theorems:

Dirichlet problem. If two functions are harmonic (Δu = 0) in finite domain Ω and coincide on its smooth boundary ∂Ω, they are identical in Ω.

Neumann problem. If two one-valued functions u1 and u2 are harmonic (Δu = 0) in a finite domain Ω and their normal derivatives coincide on its boundary ∂Ω, then u1 and u2 differ by at most a constant in Ω.

Mixed boundary-value problem. If two functions u1 and u2 are harmonic (Δu = 0) in a finite domain Ω and satisfy the same mixed boundary conditions, then u1 and u2 coincide in Ω.    ▣

Note that the above uniqueness theorems are valid also for external boundary value problems. Often, the Dirichlet boundary value problems arise with discontinuous boundary conditions. If we assume that the function u satisfies Laplace's equation inside the domain Ω, contiguously adjoining the boundary values at the points of continuity of the latter, and bounded in neighborhoods of points of discontinuity, then the solution to the first type boundary value problem is unique.

Gauss' Integral Theorem: If Δu = 0 in a finite domain Ω ⊂ ℝ³ with smooth boundary ∂Ω, then
\[ \iint_{\partial\Omega} \frac{\partial u}{\partial {\bf n}} \, {\text d}A = 0 , \]
where dA is the element of surface and n is a unit outward normal to Ω.    ▣

Mean-Value Theorem: Let P(x,y,z) be a fixed point inside a finite domain Ω ⊂ ℝ³ with smooth boundary ∂Ω, and let Δu = 0 in Ω. Let S be a sphere of radius r centered at P with boundary ∂S lying entirely inside Ω. If Δu = 0 in Ω, then
\[ u(P) = \frac{1}{4\pi r^2}\iint_{\partial S} u \, {\text d}A \qquad\mbox{and in $n$ dimensional case} \qquad u (P) = \frac{1}{s_1 r^{n-1}}\int_{S_r} u \left( Q \right) {\text d}S_Q , \]
where Sr is the sphere of radius r centered at P and s1 is the area of the unit sphere. For two dimensional case (on the plane), we have for a harmonic function (\( \phi_{xx} + \phi_{yy} = 0 \) )
\[ \phi (x,y) = \frac{1}{2\pi r}\int_0^{2\pi} \phi \left( x+ r\,\cos\theta , y + r\,\sin\theta \right) {\text d}\theta . \]
Thus, the value of a harmonic function at any point P is the average of the values it takes on any sphere surrounding that point.    ▣

The Principle of Maximum Value: If the function u(P), defined and continuous in the closed region Ω+∂Ω, satisfies Laplace's equation ∇²u = 0 inside Ω, then maximum and minimum values of the function u(P) are attained on the smooth surface ∂Ω.

Dirichlet’s Principle: Any classical solution u of the Dirichlet problem for Laplace's equation
\[ \Delta u = 0, \qquad \left. u \right\vert_{\partial\Omega} = f , \]
where Ω is a bounded domain in ℝn and f is a given function on the smooth boundary ∂Ω of Ω, minimizes the Dirichlet integral
\[ I(u) = \int_{\Omega} \left\vert \nabla u(x) \right\vert^2 {\text d}x \]
among all smooth functions taking f on the boundary ∂Ω.

It will be shown in next sections that boundary value problems can be solved without much difficulty for simple two-dimensional regions, such as rectangle, half-plane, circle, and so on, the problem is much more difficult for general domains. There are, of course, methods for approximating and solving boundary value problems such as the method of Balayage (see Tsuji's book), the Perron method, finite-difference methods, finite element method, Monte Carlo method, Galerkin method, and boundary element method to name some, but nevertheless the problem is difficult.

An important technique to solve boundary value problems is integral transformation. We will demonstrate later its applications in some relatively simple domain. For example, the Fourier integrals can be used to determine functions of Laplace's or Helmholtz operators in infinite space ℝn. The Laplace operator Δ in the space ℝn has only negative eigenvalues, so its negation will be a positive one. We can define its square root (using Fourier transformation) as
\[ \left(\left( - \Delta \right)^{1/2} f \right) ({\bf x}) = C_n \int_{\mathbb{R}^n} \frac{f({\bf x}) - f({\bf y})}{\| {\bf x} - {\bf y} \|^{n+1}} \,{\text d}{\bf y} , \]
where x ∈ ℝn and Cn is a constant.

The real and imaginary part of any holomorphic function yield harmonic functions on ℝ² (these are said to be a pair of harmonic conjugate functions). the Laplace operator on the plane can be decomposed as \( \displaystyle \nabla^2 = \left( \frac{\partial}{\partial x} + {\bf j}\,\frac{\partial}{\partial y} \right) \left( \frac{\partial}{\partial x} - {\bf j}\,\frac{\partial}{\partial y} \right) \) where j is the unit vector in positive vertical direction on the complex plane ℂ, so j² = -1. Therefore, the solution of Laplace's equation in two dimensional domains is intimately related to the theory of analytic functions of a complex variable. This is immediately seen observing that, writing z = x + jy, the complex function g(z) = uxj uy is holomorphic in Ω because it satisfies the Cauchy–Riemann equations. Hence, g has locally a primitive f, and u is the real part of f up to a constant, as ux is the real part of f′ = g. In fact, this topic occupies a significant portion of texts on complex variables, and we are not after it.

Finally, we present two important methods to solve boundary value problems in arbitrary domain. We start with Green's function method, and then go to the reduction to integral equations---the technique used in finite element method.

Although Green's function method could be applied to arbitrary boundary value problem for an elliptic operator

\[ \begin{split} L[u] &= f(P), \qquad P \in \Omega , \\ \left( a(P)\,u(P) + b(P)\,\frac{\partial u}{\partial {\bf n}} \right)_{\partial\Omega} &= \varphi (P), \end{split} \]
we will mostly demonstrate it for Laplace's equation when \( L[u] = \nabla^2 u . \) As usual, n denotes a unit outward vector to the boundary ∂Ω of the region Ω.

The Green function method in two and three dimensions is based on the formulas:

\[ \Delta \left( \ln \frac{1}{r_{PQ}} \right) = -2\pi \,\delta \left( P, Q \right) \qquad\mbox{and} \qquad \Delta \left( \frac{1}{r_{PQ}} \right) = -4\pi \,\delta \left( P, Q \right) , \]
respectively. Here rPQ is the distance between points P(x,y,z) and Q(ξ,η,ζ) from ℝ³, and \( \delta (P,Q) = \delta(x,\xi )\,\delta(y,\eta )\,\delta(z,\zeta ) \) is the delta function of Dirac. For two dimensional case, points P(x,y) and Q(ξ,η) depend on two variables.
The Green function for the Poisson equation subject to the boundary conditions
\[ L \left[ u \right] = f \left( P\right) , \qquad \left( a(P)\,u(P) + b(P)\,\frac{\partial u}{\partial {\bf n}} \right)_{P\in\partial\Omega} = 0, \]
where L is an elliptic operator (of which we consider only either the Laplacian (L = Δ) or Helmholtz operator Δ + k²), is the function G(P.Q) of two variables/points that is the solution to the specific nonhomogeneous equation and homogeneous boundary conditions:
\[ L \left[ G \right] = \delta \left( P,Q \right) , \qquad \left( a(P)\,G(P,Q) + b(P)\,\frac{\partial G}{\partial {\bf n}} \right)_{P\in\partial\Omega} = 0. \]

Once Green's function is obtained, the solution of the corresponding boundary value problem \( \mbox{div}\left( k\,\nabla \,u \right) - q\,u = f , \quad \left( a(P)\,u(P) + b(P)\,\frac{\partial u}{\partial {\bf n}} \right)_{P\in\partial\Omega} = \varphi \) can be expressed explicitly:

\[ u(Q) = \int_{\Omega} G(P,Q)\, f(M)\,{\text d}V_P - \int_{\partial\Omega} k \left( G\,\frac{\partial u}{\partial {\bf n}} - u\, \frac{\partial G}{\partial {\bf n}} \right) {\text d}S_P . \]
In particular, for the first type boundary value problem (Dirichlet problem, when a = 1 and b = 0),
\[ \mbox{div}\left( k\,\nabla \,u \right) - q\,u = f , \qquad \left. u \right\vert_{\partial\Omega} = \varphi , \]
we have its solution expressed through corresponding Green's function
\[ u(Q) = \int_{\Omega} G(P,Q)\, f(M)\,{\text d}V_P + \int_{\partial\Omega} k\varphi (P)\, \frac{\partial G}{\partial {\bf n}}\,{\text d}S_P . \]
For the second type boundary value problem (Neumann problem, when a = 0 and b = 1),
\[ \mbox{div}\left( k\,\nabla \,u \right) - q\,u = f , \qquad \left. \frac{\partial u}{\partial {\bf n}} \right\vert_{\partial\Omega} = \varphi , \qquad \int_{\partial\Omega} \varphi (P)\,{\text d}S_P = 0 , \]
we have
\[ u(Q) = \int_{\Omega} G(P,Q)\, f(P)\,{\text d}V_P - \int_{\partial\Omega} k\varphi (P)\, G(P,Q)\,{\text d}S_P . \]
For the third type boundary value problem (when a ≠ 0 and b ≠ 0),
\[ \mbox{div}\left( k\,\nabla \,u \right) - q\,u = f , \qquad \left( a(P)\,u(P) + b(P)\,\frac{\partial u}{\partial {\bf n}} \right)_{P\in\partial\Omega} = \varphi (P) , \]
we have
\[ u(Q) = \int_{\Omega} G(P,Q)\, f(M)\,{\text d}V_P - \int_{\partial\Omega} \frac{k\varphi (P)}{b(P)}\, G(P,Q)\,{\text d}S_P . \]

Green's functions are all symmetric for a self-adjoint elliptic operator \( L[u] = \mbox{div} \left( k\,\nabla\,u \right) - q\, u : \)
\[ G(P,Q) = G(Q,P) . \]
Its proof is based on application of the second Green's formula:
\[ \int_{\Omega} \left\{ u\,L[v] - v\,L[u] \right\} {\text d}V = \int_{\partial\Omega} k \left( u\,\frac{\partial v}{\partial {\bf n}} - v\,\frac{\partial u}{\partial {\bf n}} \right) {\text d}S . \]

The Green function can be constructed explicitly

\[ P(P,Q) = \frac{-1}{4\pi \,r_{PQ}} +v (P,Q) \]
for three dimensional case, and
\[ P(P,Q) = \frac{-1}{2\pi}\,\ln \frac{1}{r_{PQ}} +v (P,Q) \]
for two dimensional case. In the above equations, v(P,Q) is a harmonic functions satisfying the appropriate boundary conditions to compensate the influence of the singular term:
\[ \left( a(P)\,v(P,Q) + b(P)\,\frac{\partial v}{\partial {\bf n}} \right)_{P\in\partial\Omega} = \left( a(P)\,\frac{1}{4\pi\,r_{PQ}} + b(P)\,\frac{\partial}{\partial {\bf n}} \left( \frac{1}{4\pi\,r_{PQ}} \right) \right)_{P\in\partial\Omega} \]
for three dimensional case, and
\[ \left( a(P)\,v(P,Q) + b(P)\,\frac{\partial v}{\partial {\bf n}} \right)_{P\in\partial\Omega} = \left( a(P)\,\frac{1}{2\pi}\ln \frac{1}{r_{PQ}} + b(P)\,\frac{\partial}{\partial {\bf n}} \left( \frac{1}{2\pi}\ln \frac{1}{r_{PQ}} \right) \right)_{P\in\partial\Omega} \]
for two dimensional case. The corresponding first kind and third kind boundary value problems for function v have unique solutions. The formulas and statements are valid for outer boundary value problems as well.

Note: Green's function for the Neumann problem cannot be expressed through the above formula \( G(P,Q) = \frac{-1}{4\pi\,r_{PQ}} + v \) because such function v does not exist because v must be a solution of second kind boundary value problem

\[ \Delta v =0 , \qquad \left. \frac{\partial v}{\partial {\bf n}} \right\vert_{\partial \Omega} = \varphi (M) = \frac{1}{4\pi} \, \frac{\partial}{\partial {\bf n}} \left( \frac{1}{r_{MP}} \right)_{\partial\Omega} . \]
Such function φ does not exist since necessary condition \( \int_{\partial\Omega} \varphi \,{\text d}S =0 \) fails. In this case, Green's functions is defined upon solving the Neumann problem:
\[ \Delta G = \delta (P,Q) , \qquad \left. \frac{\partial G}{\partial {\bf n}} \right\vert_{\partial \Omega} = \frac{1}{\mbox{area of }\partial\Omega} . \qquad \triangleleft \]

Since Green's function is harmonic everywhere except at one point, and because it vanishes at infinity, it has the physical interpretation of being the potential of the electrostatic field due to a point charge. This physical interpretation is of use in constructing Green's functions for some simple domains. Here we demonstrate the method of reflection for finding Green's function for half plane in the following example.

Example: Let us consider first two dimensional case, namely upper half-plane. Green's functions for such region may be interpreted as the deflection of a membrane (which is clamped all along the x-axis) as measured at a point P = (x,y) due to a unit concentrated force at Q = (ξ,η). It is clear from symmetry that the zero boundary condition on y = 0 can be achieved by adding an image force of unit negative strength at point R = (ξ,-η). Thus, we have

\[ G(P,Q) = \frac{1}{2\pi} \,\ln r_{PQ} - \frac{1}{2\pi} \,\ln r_{PR} , \]
\[ r_{PQ} = = \sqrt{(x-\xi )^2 + (y-\eta )^2} ; \qquad r_{PR} = = \sqrt{(x-\xi )^2 + (y+\eta )^2} . \]
harrow = Graphics[Arrow[{{-2, 0}, {1.5, 0}}]];
varrow = Graphics[Arrow[{{0, -0.2}, {0, 1.6}}]];
line = Line[{{1, -1.5}, {-1.8, 0.5}, {1, 1.5}}]
ll = Graphics[{Thick, line}];
vline = Line[{{1, -1.5}, {1, 1.5}}]
vl = Graphics[{Dashed, vline}];
p = {{1, -1.5}, {1, 1.5}}
point = Graphics[{PointSize[Large], Point[p]}];
xtext = Text[Style["x", FontSize -> 14, Black], {1.6, 0.1}]
ytext = Text[Style["y", FontSize -> 14, Black], {0.1, 1.7}]
ptext = Text[Style["P = (x,y)", FontSize -> 14, Red], {-2.2, 0.55}]
qtext = Text[ Style["Q = (\[Xi],\[Eta])", FontSize -> 14, Blue], {1.5, 1.5}]
rtext = Text[ Style["R = (\[Xi],-\[Eta])", FontSize -> 14, Blue], {1.5, -1.5}]
rpqtext = Text[Style["\!\(\*SubscriptBox[\(r\), \(PQ\)]\)", FontSize -> 14, Purple], {-0.5, 1.15}]
rprtext = Text[Style["\!\(\*SubscriptBox[\(r\), \(PR\)]\)", FontSize -> 14, Purple], {-0.5, -0.6}]
Show[vl, ll, varrow, harrow, point, Graphics[{xtext, ytext, ptext, qtext, rtext, rpqtext, rprtext}]]

Neumann's Green function in the upper half-plane is simply

\[ G(P,Q) = \frac{1}{2\pi} \,\ln r_{PQ} + \frac{1}{2\pi} \,\ln r_{PR} . \]
It may be interpreted as the velocity potential due to a two-dimensional unit source at Q and an image source of equal strength at R. Now \( \partial G/\partial{\bf n} = - \partial G/\partial y , \) and it is clear that on y = 0 the normal components of velocity due to the two sources are exactly opposite and cancel out. Note that since the region is infinite, integral over the boundary \( \int_{-\infty}^{\infty} G(P,Q)\,{\text d}x \) does not exist and we have normalized G by requiring its minimum value on the boundary y = 0 to occur at P = (ξ,0).

Now we go to the three dimensional case. Consider the Dirichlet problem in half space:
\[ \nabla^2 u =0 , \quad -\infty < x < \infty , \quad -\infty < y < \infty , \quad 0 < z < \infty , \qquad u(x,y,0) = f(x,y) . \]
The solution of this boundary value problem is given by Poisson's integral:
\[ u(x,y,z) = \frac{z}{2\pi} \iint_{{\mathbb R}^2} \frac{f(\xi , \eta )\,{\text d}\xi \,{\text d}\eta}{\left[ (x-\xi )^2 + (y-\eta )^2 + z^2 \right]^{3/2}} . \]

Example: Let us consider a unit sphere and two points in spherical coordinates P = (r,θ,φ) and Q = (ρϑ,ψ). Then we choose point R = (1/ρ,ϑ,ψ) outside the sphere symmetrical to Q. If we denote the angle POQ by γ, the distance PQ by rPQ and the distance PR by rPR, we have

\[ \cos \gamma = \cos\theta \,\cos\vartheta + \sin \theta \,\sin \vartheta \,\cos \left( \varphi - \psi \right) . \]
The corresponding distances become
\[ r^2_{PQ} = r^2 + \rho^2 - 2r\rho \,\cos \gamma , \]
\[ r^2_{PR} = r^2 + \frac{1}{\rho^2} - 2\,\frac{r}{\rho} \,\cos \gamma , \]
It is easy seen that if P is on the surface of the sphere (r = 1), then
\[ \left. \frac{r_{PQ}}{r_{PR}} \right\vert_{r=1} = \rho . \]
Using this limit value, we can immediately derive the following Green's functions:
  • Interior of unite circle: \( \displaystyle G(P.Q) = \frac{1}{2\pi}\,\ln \frac{r_{PQ}}{\rho\, r_{PR}} . \)
  • Exterior of unite circle: \( \displaystyle G(P.Q) = \frac{1}{2\pi}\,\ln \frac{\rho\, r_{PR}}{r_{PQ}} . \)
  • Interior of unite sphere: \( \displaystyle G(P.Q) = \frac{-1}{4\pi} \left[ \frac{1}{r_{PQ}} - \frac{1}{\rho\,r_{PR}} \right] . \)
  • Exterior of unite sphere: \( \displaystyle G(P.Q) = \frac{-1}{4\pi} \left[ \frac{1}{r_{PR}} - \frac{\rho}{r_{PQ}} \right] . \)
  • Poisson formula for the circle: the solution to the Dirichlet problem
    \[ \Delta u \equiv u_{rr} + \frac{1}{r}\, u_r + \frac{1}{r^2}\,u_{\theta\theta} =0 \]
    subject to the Dirichlet boundary condition \( \displaystyle u(R,\theta ) = f(\theta ) \) in the interior domain rR is
    \[ u(r, \theta ) = \frac{R}{2\pi}\, \int_0^{2\pi} \,{\text d}\vartheta \,\frac{f(\vartheta ) \left( R^2 - r^2 \right)}{R^2 + r^2 -2rR\,\cos\left( \vartheta - \theta \right)} . \]
    This is called Poisson's formula for circle of radius R.
  • Poisson formula for the sphere: the solution to the Dirichlet problem
    \[ \Delta u \equiv u_{rr} + \frac{2}{r}\, u_r + \frac{1}{r^2}\,u_{\theta\theta} + \frac{1}{r^2 \sin^2 \theta} \,u_{\phi\phi} + \frac{\cot\theta}{r^2}\,u_{\theta} =0 \]
    subject to the Dirichlet boundary condition \( \displaystyle u(R,\theta , \phi ) = f(\theta , \phi ) \) in the interior domain rR is
    \[ u(r, \theta , \phi ) = \frac{R}{4\pi}\,\int_0^{\pi} {\text d}\psi \int_0^{2\pi} \sin\vartheta \,{\text d}\vartheta \,\frac{f(\vartheta , \psi ) \left( r^2 - R^2 \right)}{\left[ R^2 + r^2 -2rR\,\cos\left( \vartheta - \theta \right) \right]^{3/2}} . \]
    This is called Poisson's formula
It is possible to show by direct substitution that if u = f(r,θ,φ) is harmonic inside the sphere of radius R, then v = (R/r)f(R/r,θ,φ) is harmonic outside, and v → 0 as r → ∞.
r^2*D[R/r*f[R/r], {r, 2}] + 2*r*D[R/r*f[R/r], {r, 1}] Simplify[%]
(R^2 (2 r Derivative[1][f][R/r] + R (f^\[Prime]\[Prime])[R/r]))/r^3
because Laplace's equation with axial symmetry (uφ = 0) is
\[ u_{rr} + \frac{2}{r}\, u_r + \frac{1}{r^2}\,u_{\theta\theta} + \frac{\cot \theta}{r^2}\,u_{\theta} = 0. \]

From second Green's formula for Laplace's equation in some domain Ω ⊂ ℝ³ with smooth boundary ∂Ω, it follows

\[ u(x,y,z) = \frac{=1}{4\pi} \iint_{\partial\Omega} \left[ u\,\frac{\partial}{\partial {\bf n}} \left( \frac{1}{r} \right) - \frac{1}{r}\, \frac{\partial u}{\partial {\bf n}} \right] {\text d}S , \qquad r^2 = (x-\xi )^2 + (y-\eta )^2 + (z-\zeta )^2 . \]
Similar formula is valid for plain domain Ω ℝ²:
\[ u(x,y) = \frac{1}{2\pi} \oint_{\partial\Omega} \left[ u\,\frac{\partial}{\partial {\bf n}} \left( \ln r \right) - \ln r\, \frac{\partial u}{\partial {\bf n}} \right] {\text d} s , \qquad r^2 = (x-\xi )^2 + (y-\eta )^2 . \]
This is a very interesting formula. It expresses the value of a harmonic function u at any point P inside a domain Ω with smooth boundary in terms of the values of u and \( \partial u/\partial{\bf n} \) on ∂Ω. However, we know that knowledge of the boundary values of u along as well as its normal derivative is sufficient to determine u in the interior. Therefore, we need a way of eliminating either u or its normal derivative \( \partial u/\partial{\bf n} \) from the above formula. This is the idea for reducing a boundary value problem to an integral equation.

We work with Neumann and Dirichlet boundary value problems in two dimensions:

\[ \begin{split} u_{xx} + u_{yy} &= 0 \qquad\mbox{in } \Omega , \\ \frac{\partial u}{\partial {\bf n}} &= g \qquad\mbox{on } \partial\Omega ; \end{split} \]
\[ \begin{split} u_{xx} + u_{yy} &= 0 \qquad\mbox{in } \Omega , \\ u &= f \qquad\mbox{on } \partial\Omega . \end{split} \]
In the above problems, we assume that u, f, and g are suitably differentiable, and the boundary ∂Ω is smooth.

The mathematical tool required to recast boundary value problems as integral equations is the boundary or surface potential. Although there are different types of boundary potentials, we mostly will deal either with the single-layer boundary potential

\[ S(x,y) = \int_{\partial\Omega} \sigma (s)\,\ln \frac{1}{r}\, {\text d}s , \qquad r^2 = (x-\xi )^2 + (y-\eta )^2 ; \]
and the double-layer boundary potential
\[ D(x,y) = \int_{\partial\Omega} \mu (s)\,\frac{\partial}{\partial {\bf n}} \left( \ln \frac{1}{r}\right) {\text d}s , \qquad r^2 = (x-\xi )^2 + (y-\eta )^2 ; \]
The importance of these boundary potentials lies in the fact that first, they both satisfy Laplace's equation ∇² = 0 for any density σ(s) and μ;(s) in such a way that S(x,y) also satisfies the Neumann boundary conditions, and D(x,y) satisfies the Dirichlet boundary conditions. Namely, the following theorems are valid.
Single-layer boundary potential: Let ∂Ω be a smooth curve containing a region Ω and σ(s) be a single-layer density distribution along ∂Ω. Under these conditions the following properties hold for the single-layer potential S(x,y).
  1. S(x,y) satisfies Laplace's equation ∇²S = 0 inside and outside the boundary ∂Ω.
  2. The single-layer potential S(x,y) is continuous across the boundary ∂Ω.
  3. The outward normal derivative of S(x,y) has a jump discontinuity of 2πσ(ξ,η) across ∂Ω at the point (ξ,η) ∈ ∂Ω.
  4. The property
    \[ \lim_{(x,y)\to (\xi ,\eta )} \,\frac{\partial\,S(x,y)}{\partial {\bf n}} = \frac{\partial\,S(\xi , \eta )}{\partial {\bf n}} - \pi\,\sigma (\xi , \eta ) \]
    holds for points (x,y) approaching boundary points (ξ,η) from inside ∂Ω. Correspondingly, when points (x,y) approaching boundary points (ξ,η) from outside ∂Ω, we have
    \[ \lim_{(x,y)\to (\xi ,\eta )} \,\frac{\partial\,S(x,y)}{\partial {\bf n}} = \frac{\partial\,S(\xi , \eta )}{\partial {\bf n}} + \pi\,\sigma (\xi , \eta ) . \]

Similar statemet is valid for two dimensional single-layer boundary potential

\[ S(P) = \iint_{\partial\Omega} \sigma (Q)\, \frac{1}{r_{PQ}}\, {\text d}S_Q , \qquad r^2 = (x-\xi )^2 + (y-\eta )^2 + (z- \zeta )^2 ) = \| P-Q\|^2 , \]
except the last property that reads as
\[ \begin{split} \left( \frac{\partial\,S(x,y,z)}{\partial {\bf n}} \right)_i = \left( \frac{\partial\,S(x,y,z)}{\partial {\bf n}} \right)_0 - 2\pi\,\sigma (P) , \\ \left( \frac{\partial\,S(x,y,z)}{\partial {\bf n}} \right)_i = \left( \frac{\partial\,S(x,y,z)}{\partial {\bf n}} \right)_0 + 2\pi\,\sigma (P) . \end{split} \]
Here index i indicates approaching point at the boundary from inside, while inswz e indicatesd approacing boundary point from outside. Subscript 0 servs for the point at the boundary.
Double-layer boundary potential: Let ∂Ω be a smooth curve containing a region Ω and μ(s) be a double-layer density distribution along ∂Ω. Under these conditions the following properties hold for the double-layer potential D(x,y).
  1. D(x,y) satisfies Laplace's equation ∇²D = 0 inside and outside the boundary ∂Ω.
  2. The double-layer potential D(x,y) is continuous across the boundary ∂Ω.
  3. The outward normal derivative of D(x,y) is continuous across the boundary ∂Ω.
  4. The property
    \[ \lim_{(x,y)\to (\xi ,\eta )} \,D(x,y) = D(\xi , \eta ) + \pi\,\mu (\xi , \eta ) \]
    holds for points (x,y) approaching boundary points (ξ,η) from inside ∂Ω. Correspondingly, when points (x,y) approaching boundary points (ξ,η) from outside ∂Ω, we have
    \[ \lim_{(x,y)\to (\xi ,\eta )} \,D(x,y) = D(\xi , \eta ) - \pi\,\mu (\xi , \eta ) , \]
Similar statemet is valid for two dimensional double-layer boundary potential
\[ D(x,y) = \iint_{\partial\Omega} \mu (s)\,\frac{\partial}{\partial {\bf n}} \left( \frac{1}{r}\right) {\text d}s , \qquad r^2 = (x-\xi )^2 + (y-\eta )^2 \]
except the last property that reads as
\[ \begin{split} D_i (x,y,z) = D(P) - 2\pi\,\mu (P) , \qquad P(x,y,z) \in \partial \Omega , \\ D_e (x,y,z) = D(P) + 2\pi\,\mu (P) , \qquad P(x,y,z) \in \partial \Omega , \end{split} \]
where Di(x,y,z) is the limiting value of the potential of double layer approaching the point P from inside, and De(x,y,z) is the limiting value of the potential of double layer approaching the point P from outside.


The potentials of a sinle and double layer at points of the surface ∂Ω are improper integrals. These integrals converge for a special class of C¹ surfaces called Lyapunov surfaces.
A surface is called a Lyapunov surface, if the following conditions are fulfilled.
  1. At every point of surface ∂Ω there exists a unique normal and tangent surface.
  2. There exists a number d > 0 such that straight lines, parallel to the normal at some point P of the surface ∂Ω intersect a part SP of the surface ∂Ω lying inside a sphere of raius ,i.d,/i. with center P, not more than once. These parts of the surface SP are called the Lyapunov Neighbourhoods.
  3. The angle γ(P,Q) between normal vectors at pomits P and Q satisfies the following inequality
    \[ \gamma (P,Q) < A\,r^{\delta} , \]
    where r is the distance between the points P and Q, A is some positive constant, and 0 < δ ≤ 1.

Example: Assuming that the density is constant μ(&zi;,η) ≡ 1, then the double layer potential will be

\[ D(x,y) = \begin{cases} 2\pi , & \ (x,y) \mbox{ inside } \partial\Omega , \\ \pi , & \ (x,y) \mbox{ on } \partial\Omega , \\ 0, & \ (x,y) \mbox{ outside } \partial\Omega . \end{cases} \]
So it experiebces a jump of discontinuity of 2π across ∂Ω.    ■
We are now in a position to convert the Neumann and Dirichlet problems to integral equations of Fredholm type. So consider the Dirichlet boundary value problem
\[ \begin{split} \nabla^2 u &= \qquad \mbox{on } \Omega , \\ u&=f \qquad \mbox{on } \partial\Omega , \end{split} \]
where Ω is a general domain in two dimensions with Lyapunov boundary ∂Ω. We seek its solution as a double layer potential
\[ D(x,y) = \int_{\partial\Omega} \mu (s)\,\frac{\partial}{\partial {\bf n}} \left( \ln \frac{1}{r}\right) {\text d}s , \qquad r^2 = (x-\xi )^2 + (y-\eta )^2 . \]
Since D(x,y) automatically satisfies Laplace's equation in Ω for any dipole density μ(s), we seek μ(s) so that D(x,y) also satisfies the boundary condition D = f. To do this, recall that D(x,y) has a jump discontinuity across the boundary ∂Ω.

Since we require the solution of the Dirichlet problem to be continuous out to the boundary, we can replace the boundary condition D = f by

\[ \lim_{(x,y)\to (\xi ,\eta )} \,D(x,y) = D(\xi , \eta ) + \pi\,\mu (\xi , \eta ) = f(\xi , \eta ). \]
Then we get
\[ D(\xi , \eta ) + \pi\,\mu (\xi , \eta ) = f(\xi , \eta ) , \]
and finally, replacing D(ξ,η) by its integral form, we arrive at
\[ \int_{\partial\Omega} \frac{\cos\varphi}{r}\,\mu(s)\,{\text d}s + \pi\,\mu (\xi , \eta ) = f(\xi , \eta ) , \]
where φ is the angle between the normal at the fixed point P(ξ,η) and the line PQ connecting P with varing point Q of integration. The above equation is a Fredholm integral equation of the second kind.

Similarly, we can solve a Neumann problem using a single layer potential S(x,y). Then using jump discontiuity property of its normal derivative, we arrive at the integral equation

\[ \int_{\partial\Omega} \sigma (s)\,\ln \frac{1}{r}\, {\text d}s -\pi\,\sigma (\xi , \eta ) = g(\xi , \eta ) , \]
from which we cansolve for σ(s).    ■


  1. Gustafson, K. and Abe, T., The third boundary condition—was it robin’s? The Mathematical Intelligencer, March 1998, Volume 20, Issue 1, pp 63--71.
  2. Kellogg, O.D., Foundations of potential theory, 1929, Springer, Berlin.
  3. Markina, Irina, Potential theory: the origin and applications, Department of Mathematics, University of Bergen, Johannes Brunsgate 12, Bergen 5008, Norway.
  4. Seeley, R.T., Complex powers of an elliptic operator, American Mathematical Society, 1967, 20 pages.
  5. Seeley, R.T., The resolvent of an elliptic boundary problem, American Journal of Mathematics, 1969, Vol. 91, No 4, pp. 889--920.
  6. Tsuji, M., Potential theory in mordern function theory, Chelsea, reprint, 1975.
  7. A. D. Wentzell (Ventcel’), On boundary conditions for multidimen-sional diffusion processes, Teoriya Veroiatnostei i ee Primeneniya, translation from Russian in Theory of Probability & Applications, 1959, Vol. 4, Issue 2, pp. 164–177.
  8. Yu, C.L., Reflection principle for solving of higher order elliptic equations with analytic coefficients, SIAM Journal on Applied Mathematics, 1971, Vol. 20, No. 3, pp.358--363.


Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions