The following sections are devoted to Laplace and Helmholtz equations as typical representatives of the elliptic partial differential equations. We mostly deal with the plane case when the number of independent variables is restricted to be two. However, some three dimensional cases are also included in the tutorial. Elliptic differential equations occur in many applications; its theory contains many important results of which only few are included.

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

Elliptic Equations

Let A(x,y), B(x,y), and C(x,y) be known differentiable functions in some domain Ω ⊂ ℝ² with smooth boundary ∂Ω. The differential operator of second order

\[ L[u] \equiv \left( A\,u_x \right)_x + \left( B\,u_y \right)_x + \left( B\,u_x \right)_y + \left( C\,u_y \right)_y \]
is called elliptic if these functions satisfy the inequalities:
\[ A > 0 \qquad\mbox{and} \qquad B^2 < A\,C . \]
Here and throughout subscripts specify partial derivatives, so \( u_x = \partial u/\partial x \) and so on. The most general linear elliptic second order partial differential equation in n independent variables is
\[ \sum_{i,j=1}^n a_{ij}({\bf x}) \,\frac{\partial^2 u}{\partial x_i\,\partial x_j} + \sum_{k=1}^n b_k ({\bf x}) \,\frac{\partial u}{\partial x_k} + c({\bf x})\,u = f({\bf x}) , \]
where the symmetric matrix \( \left[ a_{ij}({\bf x}) \right] \) is positive definite, and boldface is used for vectors.

Typical examples of elliptic equations give famous Laplace's equation

\begin{equation} \Delta u =0 \qquad \mbox{or} \qquad \nabla^2 u \equiv \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} =0 \label{EqElliptic.1} \end{equation}
in two dimensions and
\begin{equation} \Delta u =0 \qquad \mbox{or} \qquad \nabla^2 u \equiv \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} =0 \label{EqElliptic.2} \end{equation}
in three dimensions. Other examples of elliptic equations include the Helmholtz equation
\begin{equation} \Delta u + \omega^2 u =0 , \label{EqElliptic.3} \end{equation}
where ω² is a given function; and equations generated by the powers of the Laplacian such as the biharmonic equation \( \Delta^2 u = 0 . \) Solutions of Laplace's equation are referred to as harmonic functions. Note that the operator ∇² is commonly written as Δ by mathematicians and it is called the Laplacian or the Laplace operator. Their nonhomogeneous counterparts are usually refer to as Poisson's equation; so \( \nabla^2 u = f \) and \( \nabla^2 u + \omega^2 u = f \) are their corresponding examples. Some general properties of elliptic equations are formulated for the two dimensional case below.
Dirichlet principle: If u(x,y) is a two-dimensional harmonic function defined over Ω ⊂ ℝ² with smooth boundary ∂Ω, then at all points of Ω, \( u_{xx} + u_{yy} =0 . \) Let v(x,y) be any non-harmonic smooth function that has the same prescribed boundary conditions as u on ∂Ω so that at all points of ∂Ω, u = v. If v possesses first order derivatives in Ω, then
\[ \iint_{\Omega} \left( u_x^2 + u_y^2 \right) {\text d}x\,{\text d} y \le \iint_{\Omega} \left( v_x^2 + v_y^2 \right) {\text d}x\,{\text d} y \]

This principle can be generalized for arbitrary elliptic operator,

Generalized Dirichlet principle: If u(x,y) is a two-dimensional function that is defined over Ω ⊂ ℝ² with smooth boundary ∂Ω and satisfy the equation
\[ L[u] \equiv \left( A\,u_x \right)_x + \left( B\,u_y \right)_x + \left( B\,u_x \right)_u + \left( C\,u_y \right)_y =0 \]
at all points of Ω. Let v(x,y) be a smooth function that has the same prescribed boundary conditions as u on ∂Ω so that at all points of ∂Ω, u = v. If v possesses first order derivatives in Ω, then
\[ \iint_{\Omega} \left( A\,u_x^2 + 2B\,u_x u_y + C\,u_y^2 \right) {\text d}x\,{\text d} y \le \iint_{\Omega} \left( A\,v_x^2 + 2B\,v_x v_y + C\,v_y^2 \right) {\text d}x\,{\text d} y . \]

Statement: Consider the linear elliptic equation \begin{equation} \sum_{i,j =1}^n a_{ij} ({\bf x})\,u_{x_i x_j} + \sum_{i=1}^n b_i ({\bf x})\,u_{x_i} + c({\bf x})\,u = f({\bf x}), \label{EqElliptic.4} \end{equation} where x = (x1, x2, … , xn) ∈ Ω ⊂ ℝn, the coefficients 𝑎ij, bi, c, and function f = f(x) are defined in an open domain Ω. If \[ c = c({\bf x}) < 0 \qquad\mbox{and} \qquad f = f({\bf x}) \ge 0 \quad (\mbox{or } \le 0) \] hold in Ω, then the solution u = u(x) of Eq.\eqref{EqElliptic.4} regular in Ω cannot attain a positive maximum (or negative minimum) in Ω. If \[ c \le 0 \qquad \mbox{and} \qquad f > 0 )\mbox{or } < 0) \] holds in Ω, then u = u(x) cannot attain a non-negative maximum (or non-positive minimum) in Ω.

Example: Consider the Helmholtz equation (42, page 76) \[ u_{xx} + u_{yy} + 2u = 0 \qquad (0 < x < \pi , \quad 0 < y < \pi ). \tag{Helmholtz.1} \] Then the function u = u(x, y) = sin(x) sin(y) is a solution of Eq.(Helmholtz.1) in the square S = [0, π} × [0, π] and vanish on the perimiter of S without vanishing inside S.    ■

Example: Consider the Helmholtz equation (42, page 184) \[ u_{xx} + u_{yy} - 4u = 0 \qquad (-\infty < x < \infty , \quad 0 < y < \infty ). \tag{Helmholtz.2} \] in half plane \[ \Omega : \quad y \ge 0 \quad \mbox{and} \quad x \in \mathbb{R} = (-\infty , \infty ). \] Eq.(Helmholtz.2) has a regular solution in Ω that vanish on the boundary ∂Ω :    y = 0 and x ∈ ℝ without identically vanishing in Ω.

In gact, such a solution is of the form \[ u = u(x,y) = \sinh \left( x \sqrt{2} \right) \sinh \left( y \sqrt{2} \right) . \]    ■

Example: Consider the Laplace equation \eqref{EqElliptic.1} \( u_{xx} + u_{yy} = 0 \) in the unbounded domain (16, page 20) \[ \Omega = \left\{ (x,y) \in \mathbb{R}^2 \,: \ -\infty < x < \infty , \quad 0 < y < \pi \right\} . \] This Laplace equation has regular solution \[ u = u(x,y) = e^x \,\sin y \] in Ω that vanishes on the boundary ∂Ω : y = 0 and y = π,     −∞ < x < ∞, without identically vanishing in Ω.    ■

Green's second identity: If u(x,y) and v(x,y) are both twice continuously differentiable over Ω ⊂ ℝ² with smooth boundary ∂Ω, then
\[ \iint_{\Omega} \left( u\,\nabla^2 v - v\,\nabla^2 u \right) {\text d}x{\text d}y = \oint_{\partial\Omega} \left( u\,\frac{\partial v}{\partial n} - v\, \frac{\partial u}{\partial n} \right) {\text d}S , \]
where ∂/∂n is a derivative along the outward normal. This identity is named after the British mathematician George Green (1793--1841); it can be extended for a three-dimensional case.

This formula can be generalized for arbitrary self-adjoint elliptic operator \( L[u] = \mbox{div}\left( k\,\nabla \,u \right) - q\, u \) and for biharmonic operator.

Generalized Green's second identity:
\[ \int_{\Omega} \left( u\,L[v] - v\,L[u] \right) {\text d}V = \oint_{\partial\Omega} k \left( u\,\frac{\partial v}{\partial n} - v\, \frac{\partial u}{\partial n} \right) {\text d}S , \]
for \( L[u] = \mbox{div}\left( k\,\nabla \,u \right) - q\, u \) and
\[ \int_{\Omega} \left( u\,\nabla^4 v - v\,\nabla^4 u \right) {\text d}V = \oint_{\partial\Omega} \left( u\,\frac{\partial}{\partial n}\,\nabla^2 v - \nabla^2 v\, \frac{\partial u}{\partial n} + \nabla^2 u\,\frac{\partial v}{\partial n} -v\,\frac{\partial}{\partial n}\, \nabla^2 u \right) {\text d}S . \]

Elliptic partial differential equations are typically accompanied by boundary conditions. To be more specific, let Ω be domain (finite or infinite) in n-dimensional space ℝn with smooth boundary ∂Ω. There are known several boundary conditions, out of them we mostly concentrate on three of them.

Dirichlet's boundary conditions or boundary conditions of the first kind:

\[ \left. u({\bf x}) \right\vert_{{\bf x}\in \partial \Omega} = f({\bf x}) . \]
Neumann boundary conditions or boundary conditions of the second kind:
\[ \left. \frac{\partial u({\bf x})}{\partial {\bf n}} \right\vert_{{\bf x}\in \partial \Omega} = f({\bf x}) , \]
where n is a unit outward normal vector to the boundary ∂Ω.

Third kind boundary conditions (that are also frequently called mixed):

\[ \left. \left[ \alpha ({\bf x})\,u({\bf x}) + \beta ({\bf x})\,\frac{\partial u({\bf x})}{\partial {\bf n}} \right]\right\vert_{{\bf x}\in \partial \Omega} = f({\bf x}) , \]
where α and β are some smooth functions on the boundary ∂Ω. These boundary conditions are mistakenly associated with the French mathematician Victor Gustave Robin (1855–1897). There is no evidence that Robin ever used the boundary conditions of third kind in his research. This term "Robin boundary conditions" was first introduced by the Russian mathematician Sergei Natanovich Bernstein (1880--1968), who most likely had difficulties with reading articles in French (remember, Google was not available to him; see a historical survey by Gustafson and Abe). The boundary value problem satisfying a mixed boundary conditions was first solved in 1910 by the Polish mathematician Stanisław Zaremba (1863--1942) for the Laplace equation.

When we consider heat transfer problems, the Dirichlet boundary conditions correspond to prescribing temperature on the boundary. The Neumann boundary conditions correspond to prescribing heat flux on the boundary. The convective boundary conditions leads to the third kind boundary conditions. An elliptic partial differential equation with one of corresponding boundary conditions is called the boundary value problem. Out of these, there are two important classes of boundary value problems: interior or inner boundary value problems and exterior or outer boundary value problems. They are ussually referred to a finite domain Ω ⊂ ℝn and a solution is seeking either inside Ω or outside it.

There are two main reasons for considering these three kinds of boundary conditions. First, they are mostly used in practical applications; and second, the corresponding boundary value problems have solutions. These boundary value problems also have unique solutions except Neumann's conditions. On the plane (the number of independent variables is 2), the Neumann boundary value problem for the Laplace equation is determined up to arbitrary additive constant. Moreover, the boundary function f(x,y) should also satisfy some auxiliary condition for existence of the solution. In multi-dimensional case, the situation is slightly different, but still requires some additional conditions. All these statements are valid only when the boundary ∂Ω of the finite domain Ω is smooth. When the boundary has singular points where its curve/surface is not differentiable (we call such points wedge or corner points), then the given boundary value problem may have multiple solution, and additional conditions must be imposed to guarantee its uniqueness. On the plane, it is required that the solution should be bounded when a point approaches the corner. Such condition is called the wedge condition. In arbitrary n-dimensional case (when n ≥ 3), these conditions will be different from the plane case. In general, the lack of boundary regularity or the transmission conditions result in singular behavior of the solution near such point. The boundary value problems with singular points on the boundary were studied extensively. A thorough survey can be found in Kondrat'ev--Oleinik papers. When tangent lines from both sides of the corner point coincide (as for instance in the heart shape domain ❤), the boundary value problem has no solution.

Consider a domain \( \Omega \subset \mathbb{R}^n \) in n dimensional space, and a function \( u\,:\, \overline{\Omega} \to \mathbb{R} \) defined on Ω and its boundary. We assume that the boundary \( \partial\Omega \) of the domain and the function are sufficiently smooth. In particular, the boundary should not have corner points---this case will be analyzed separately.

The Dirichlet integral (functional) of u over Ω is defined by

\[ {\cal F} (u) = \frac{1}{2} \int_{\Omega} \left\vert \nabla \,u({\bf x}) \right\vert^2 {\text d} {\bf x} . \]
We try to minimize the Dirichlet integral in the space of functions that satisfy either the Dirichlet boundary condition
\[ u = f \qquad \mbox{on } \partial \Omega \]
or the Neumann boundary condition
\[ \frac{\partial u}{\partial {\bf n}} = f \qquad \mbox{on } \partial \Omega . \]
Here f is a given function defined on the boundary \( \partial\Omega \) of Ω, and n is the exterior normal unit vector to the boundary.

Let \( h\,:\,\overline{\Omega} \to \mathbb{R} \) be a smooth function that satisfies either homogeneous Dirichlet boundary condition \( \left. u \right\vert_{\partial \Omega} =0 \) or homogeneous Neumann condition \( \left. \partial u/\partial{\bf n} \right\vert_{\partial \Omega} =0 . \) Since the further derivations only slightly differ from each other in the case of Dirichlet boundary condition and the Neumann boundary condition, we consider for simplicity only the latter one. Calculation of the variational derivative of the Dirichlet functional yields

\[ {\text d} {\cal F} (u)\, h = \lim_{\varepsilon \to 0} \,\frac{\text d}{{\text d}\varepsilon}\,\frac{1}{2} \int_{\Omega} \left\vert \nabla u + \varepsilon \nabla h \right\vert^2 {\text d} {\bf x} = \int_{\Omega} \nabla u \cdot \nabla h \, {\text d} {\bf x} . \]
Thus, any minimizer of the Dirichlet functional must satisfy
\[ \int_{\Omega} \nabla u \cdot \nabla h \, {\text d} {\bf x} =0 \]
for all smooth functions h that vanish on the boundary.

Using the identity
\[ \nabla \cdot \left( h \nabla u \right) = h\,\Delta u + \nabla u \cdot \nabla h \]
and the divergence theorem, we get
\[ \int_{\Omega} \nabla u \cdot \nabla h \, {\text d} {\bf x} = - \int_{\Omega} \left( \Delta u \right) h \, {\text d} {\bf x} + \int_{\partial\Omega} h\, \frac{\partial u}{\partial {\bf n}} \, {\text d} S . \]
Since h=0 on the boundary, the integral over the boundary is zero, and we have
\[ {\text d} {\cal F}(u) \, h = - \int_{\Omega} \left( \Delta u \right) h \, {\text d} {\bf x} . \]
Thus, the variational derivative of the Dirichlet functional, defined by
\[ {\text d} {\cal F}(u) \, h = \int_{\Omega} \frac{\delta {\cal F}}{\delta u} \, h \, {\text d} {\bf x} , \]
is given by
\[ \frac{\delta {\cal F}}{\delta u} = - \Delta u . \]
Therefore, a smooth minimizer u of the Dirichlet functional satisfies Laplace’s equation
\[ \Delta u = 0 . \]
Similarly, a minimizer of the functional
\[ {\cal F} (u) = \int_{\Omega} \left\{ \frac{1}{2} \left\vert \nabla u \right\vert^2 - f\,u \right\} {\text d} {\bf x} , \]
where \( f\,:\,\overline{\Omega} \to \mathbb{R} \) is a given function, satisfies Poisson’s equation
\[ \Delta u = f \qquad \mbox{in } \Omega . \qquad \blacksquare \]
      Consider a domain Ω consisting of two adjacent domains Ω1 and Ω2, separated by set S (see Figure 1), so that
\[ \Omega = \Omega_1 \cup \Omega_2 \cup S . \]
Assume that the closure of S is a surface of class C¹. (42, pages 235--237).
domain = Plot[{x^{1/20}*(2*Pi - x)*Sin[x/2]/15 + x^2*Sin[x/2]/20, -x^{1/20}*(2*Pi - x)*Sin[x/2]/15 - x^2*Sin[x/2]/20}, {x, 0, 2*Pi}, PlotRange -> {{-0.15, 2*Pi + 0.5}, {-1.1, 1.1}}, Axes -> False, PlotStyle -> {{Thickness[0.01], Black}, {Thickness[0.01], Black}}, Filling -> {1 -> {2}}];
line = Plot[0, {x, 0, 2*Pi}, PlotStyle -> {Thickness[0.01], Blue}, Axes -> False];
ar1 = Graphics[{Black, Arrowheads[0.08], Arrow[{{-0.1, 0}, {2*Pi + 0.5, 0}}]}];
ar2 = Graphics[{Black, Arrowheads[0.08], Arrow[{{0, -1.1}, {0, 1.1}}]}];
t1 = Graphics[ Text[Style[ ToString[Subscript[\[CapitalOmega], 1], TraditionalForm], FontSize -> 18], {5.4, 0.3}]];
t2 = Graphics[ Text[Style[ ToString[Subscript[\[CapitalOmega], 2], TraditionalForm], FontSize -> 18], {5.4, -0.3}]];
s = Graphics[ Text[Style[ToString["S", TraditionalForm], FontSize -> 18], {3.4, 0.2}]];
x = Graphics[ Text[Style[ToString["x", TraditionalForm], FontSize -> 18], {6.5, 0.2}]];
y = Graphics[ Text[Style[ToString["y", TraditionalForm], FontSize -> 18], {0.3, 1.02}]];
Show[domain, line, ar1, ar2, t1, t2, s, x, y]
       Domain Ω = Ω1 ∪ Ω2S.            Mathematica code

Theorem: A function u = u(x, y) ∈ C²(ω) harmonic in each domain Ω1 and Ω2 satisfies the Laplace equation 7nabla;²u = 0 in Ω.

Example: Let S = (0, 1) be the interval on the abscissa and Ω1 = y > 0; correspondingly, Ω2 = y < 0. Let u(x, y) satisfies Laplace's equation uxx + uyy = 0 and the initial conditions

\[ \left. \begin{array} u(x,0) = 0 , \\ u_y (x,0) = f(x), \end{array} \right\} \]
where f(x) is a continuously differentiable function not analytic in the interval (0, 1). The boundary ∂Ω1 contains the segment S. According to the reglection principle, the function u(x, y) can be extended in anti-symemtrix way
\[ u(x,y) = - u(x,-y) \]
into the lower half-plane Ω2. As a result, we get a harmonic function in Ω = ℝ². But is the function is harmonic in Ω, it should be harmonic and analytic within the interval (0, 1). This leads to conclusion that f(x) is an analytic function---contradiction our assumption. As a particular example, you may want to choose f(x) = |x|.    ■


  1. Elliptic problems on polyhedral domains, TAMU.
  2. Grisvard, P.: Elliptic problems in non-smooth domains. London: Pitman 1985
  3. Karl Gustafson and Takehisa Abe. The third boundary condition—was it robin’s? The Mathematical Intelligencer, March 1998, Volume 20, Issue 1, pp 63--71.
  4. Kondratiev, V.A., Boundary value problems for elliptic equations in domains with conical or angular points, Transactions of the Moscow Mathematical Society, 1967, Vol. 16, 227--313.
  5. Kondrat'ev, V.A., Oleinik, O.A.: Boundary-value problems for partial differential equations in non-smooth domains. Russian Mathematical Survey, 1983, Vol. 38, No. 2, pp. 1–86.
  6. Moon, P., and Spencer, D.~E., The meaning of the vector Laplacian, Journal of The Franklin Institute, 1953, 256, Issue 6, 551--558.
  7. Schechter, M., General boundary value problems for elliptic partial differential equations, Bulletine of the American Mathethatical Society, 1959, Volume 65, Number 2, pp. 70--72.
  8. Serkh, K. and Rokhlin, V., On the solution of the Helmholtz equation on regions with corners, Proceedings of the National Academy of Sciences of the United States of America, 2016, Vol. 113, pp. 9171--9176.
  9. Spence, E.A., Boundary Value Problems for Linear Elliptic PDEs, 2009, Dissertation, Cambridge.
  10. Zhevandrov, P. and Merzon, A., On the Neumann problem for the Helmholtz equation in a plane angle, Journal of Applied Mathematics and Mechanics, 2000, Vol. 23, pp. 1401--1446.


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