Chapter 3

Linear Second-Order Equations

In this chapter, we derive the heat equation in general dimension, the equation of wave propa-

gation in an elastic membrane, and the Poisson's equation. We present the properties of each

equation and then classify the linear second-order partial diﬀerential equations into elliptic,

parabolic, and hyperbolic equations. The main reason for t his classiﬁcation is that ea ch class

of equation has distinct solutions that behave diﬀerently. For instance, elliptic equations

have s o luti o ns that are smooth and have no singularities, whereas hyperbolic equations have

solutions that have s ingularities and can exhibit wave-like b eh avior. Understanding these

diﬀerences is crucial for solving and interpreting the solutions of partial diﬀerential equations

in various ﬁelds of science and engineering.

3.1 Heat equation

The heat equation is a fundamental equation in physics, engineering, and applied mathe-

matics that is used to model a wide range of phenomena, inclu ding heat transfer, thermal

diﬀusion, and mass diﬀusion. In this section, we will ﬁrst de rive the heat equation for solid

conductive continua, and then i ntroduce the advection-diﬀusion equation for ﬂuids.

3.1.1 Derivation of heat eq u ation

In the ﬁrst chapter, we demonstrated the derivation of the equation for heat ﬂow in a

conductive rod. In this chapter, we pres ent the derivation in three-dimensional space (R

)

to gain a deeper understanding of it s physics.

Consider a solid conductive continuum Ω⊂R

, which can be a bounded or unbounded set,

with or without a boundary. If Ω is a bounded domain, we suppose that the only possibility

that Ω exchange heat with the ambient space is through its boundary bnd(Ω). Let u(x; t)

denote th e temperature of an arbitrary point x2Ω at time t. For any bounded region D ⊂Ω,

the thermal energy inside D at time t is deﬁned by the integral:

(t) =

ZZZ

q(x; t) dV ;

where q is the density function of the heat energy. Therefore, the increme nt in energy can

be expressed as:

(t + δt) ¡Q

(t) =

ZZZ

fq(x; t + δt) ¡ q(x; t)gdV :

In the limiting case a s δt!0, we obtain the fo llowing relation:

(t) =

ZZZ

ρc@

u(x; t) dV ;

where ρ is the density and c is the speciﬁc heat capacity of the continu um.

On the other hand, we can relate the rate o f change of thermal energy inside D to the heat

exchange through its boundary bnd(D). The rate of heat exch a nge through the boundary

of D is proportional to the temperature gradient, according to Fourier's law , which states

that th e heat ﬂux J(x; t) at a point x in D is given by J(x; t)=¡ α(x)ru(x; t), where α(x)

is the thermal conductivity at x. Note that the negative sign arises because the direction

of the tempe rature gradient is towards increasing temperature, whereas the direction of the

heat ﬂux is from hotter to colder regions. T he ne t ﬂux leaving D through bnd(D) is equal

to the surface integral



bnd(D)

J(z; t) ·ν(z) dS = ¡



bnd(D)

α(z) ru(z; t) ·ν(z) dS ;

where ν(z) is the outward-pointing unit normal vector at z 2 bnd(D). By the conservation

of energy principle, we have

(t) = ¡



bnd(D)

J(z; t) ·ν(z) dS ;

where the negative sign accounts for the fact that the integral measures the ﬂux leaving D

at time t. Substituting

by the triple integral and rearranging, we obtain

ZZZ

ρc@

u(x; t) dV =



bnd(D)

α(z) ru(z; t) ·ν(z) dS:

Using Gauss's theorem, we can convert the surface integral in the right-hand side of the

previous relation to a volume integral. Applying the theorem yields:

ZZZ

div (α (x)ru(x; t))dV =



bnd(D)

α(x)ru(x; t) ·ν(x)dS;

and y rearranging the terms, we arrive at the integral form of the heat equation:

ZZZ

Ω

[ρc@

u(x; t) ¡div (α(x) ru(x; t))] dV = 0:

Since this equation holds for arbitrary Ω, we obtain the follow ing diﬀerential form of the

heat equation

(x; t) =

ρc

div (α(x) ru(x; t)):

2 Linear Second-Order Equations

If α is constant, then we have

(x; t) = k ∆u;

where k =

ρc

> 0.

Remark 3.1. The he a t equation that we derived invo lves a ﬁrst-order partial derivative of u

with respect to time. Therefore, to fully describe the heat dynamics, we need to specify the

initial condition u(x; 0) for x 2 Ω. This is called the initial condition o f the heat equation.

In fact, the initial condition provides the initial thermal energy distribution for the sy stem,

which drives its subsequent behav ior. Additi o nally, in some physical systems, there may be

an internal source or sink of thermal energy that is not accounted for in our de rivation. In

the presence of such a source or sink, the heat equation takes the fo rm:

ρc

div (αru) + h;

where the function h represents the rate at which the internal source or sink produces or

absorbs energy.

3.1.2 Convection-diﬀusion equation

The convection-diﬀusion equation is a partial diﬀerential equation that is widely used to

model heat and mass transfer in ﬂuids, such as in chemical engineering and environmental

science. The equation takes into account both convection, which is the transport of heat

or mass by the ﬂuid ﬂow, and diﬀusion, which is the transpo rt of heat or mass due to a

concentration gradient.

In the modeling of the heat problem, we assumed th a t heat is transferred only by conduc-

tivity. This assumption make s sense if the medium is solid. If the medium is a ﬂuid or gas,

the heat can be transferred in the form of convection i n addition to diﬀusion. In this scenario,

the heat ﬂux J consists of two terms, the conductivity term and the convection term:

J = J

cond

+ J

conv

If the ﬂuid moves with velo city V , then J

conv

= uV , and thus



bnd(D)

J ·ν dS = ¡



bnd(D)

αru ·ν dS +



bnd(D)

uV ·ν dS:

Therefore, the heat equation reads

ρcu

+div (uV ) =div (αru):

The above equation is called the convection-diﬀusion equation. The ﬁrst term on the left-

hand side of the equation represents the change in temp erature with respect to time. The

second term represents the convective transport of heat by the ﬂuid ﬂow, while the third

term represents the diﬀusive transpo rt due to a co nc entration gradient. In the presence of

a source or sink, the heat equation takes the form:

ρcu

+div (uV ) =div (αru) + h:

3.1 Heat equation 3

The convection-diﬀusion is a fundamental equation in many areas of science and engineering

and has many applications, such as modeling air and water pollution, heat transfer in pipes,

and drug delivery in biological systems.

3.1.3 Boundary conditions for bounded domains

Continuing with o ur discussion of the heat equation on bounded domain, it is worth noting

that the solution to a partial d iﬀerential equation, and especially to the heat equation,

depends heavily on the boundary conditions imposed on the bounded domain Ω.

Recall that the Newton cooling law relates the te mperature inside and outside a region

Ω. Speci ﬁcally, if T denotes the tempe rature outside Ω, the law states that the heat ﬂux

through bnd(Ω), is proportional to the diﬀerence T ¡u, where u is the temperature function

inside Ω . This relationship is expressed mathematically by the Robin boundary condition

= κ(T ¡u);

where ν denotes the outward normal to the boundary, and κ is the conductivity factor of

the boundary. When κ = 0, the bounda ry is perfectly insulated, and the condition reduces

to the hom o g en eous Neumann boundary condition

= 0. More generally, the Neumann

boundary condition takes the form

= g on bnd(Ω), where g is a given function. In the

limiting case where κ is very large, the Ro bin condition can be rewritten as

T ¡u =

;

and taking the limit κ!1 yields the condition u =T on bnd(Ω), which is called the Dirichlet

boundary condition.

Dirichlet boundary condition. The Dirichlet boundary condition is used when the

temperature at the boundary of the domain Ω is known, such that u(z; t)= g(z; t) for

z 2bnd(Ω). This condition is named after the German mathematician Peter Gustav

Lejeune Dirichlet. For a heat equation, this means that the temperature at bnd(Ω) is

kept a t the known value g(z; t). If g is identically zero, then the condition is call ed

the homogeneous Dirichlet condition.

u = g

PDE inside Ω

In one-dimensional problems, such as a conductive rod extending from x = 0 to

x = L, the condition becomes u(0; t) = g

and u(L; t) = g

, where g

and g

can be

constants or functions of tim e.

4 Linear Second-Order Equations

Neumann boundary condition. This condition is named after the German math-

ematician Carl Neumann. Mathematically, the Neumann bou ndary c ondition is

expressed as

:= ru ·ν = g;

where

is the direct ional derivative of u along the unit normal direction ν on

bnd(Ω). In the context of a heat problem, this condition speciﬁes the heat ﬂux through

the boundary bnd(Ω) in terms of function g. In part icular, if g is identically zero,

the boundary is perfectly insulated.

= g

PDE inside Ω

For a one-dimensional problem where the co nductive rod extends from x = 0 to

x = L, the condition reads:

(0; t ) = g

; u

(L; t) = g

Robin's or mixed boundary condition. As we observed above, the Robin boundary

condition is a straightforward derivation of the Newton cooling law in the context of

a heat problem. Mathematically, it can be expressed as

αu (z; t) + β

(z; t) = g(z; t):

where α and β are some coeﬃcients, and g is a given function. In particular, if g

is identically zero, the condition is called the homogeneous Robin boundary condi-

tion. The Robin boundary condition is commonly used in problems that involve heat

transfer between a solid and a ﬂuid, where the heat transfer coeﬃcient (i.e., α / β)

depends on the properties of the ﬂuid and the surface of the solid. It is also used

to model heat tran sfer in biological tissues, where the heat transfer coeﬃcient can

depend on the blood ﬂow rate a nd other physiological factors. The Robin boundary

condition is named after the French mathematician Victor Gustave Robin.

αu + β

= g

PDE inside Ω

For a one-dimensional problem, the condition reads



u(0; t) + β

(0; t ) = g

u(L; t) + β

(L; t) = g

3.1 Heat equation 5

3.1.4 Fundamental solution of heat equatio n

Consider a con ductive rod of diﬀusivity factor k > 0, extended over the entire real line (¡1;

1). Suppose that an inﬁnitely small portion of the rod at x = 0 has a temperature of u = 1,

and the temperature i s zero everywhere else. As time progresses, the heat diﬀuses across the

rod, but the total thermal energy remains constant as it was at t = 0. It can be shown that

the heat proﬁle evolves according to the function:

Φ(x; t) =

4πkt

4kt

;

for t > 0. In fact, one can verify that Φ

= k Φ

through di rect calculation. The graph of

Φ(x;t) demon strates the evolution of heat proﬁle for t > 0. The following ﬁgure shows Φ(x;t)

for t = 0.1; 0.5 a nd 1 where k = 1:

-4 -2 0 2 4

0.2

0.4

0.6

0.8

That the total thermal energy remains constant over time follows from the fact that the

integral of Φ for a ny ﬁxed value o f t is constant:

¡1

Φ(x; t)dx = 1:

In fact, using the change of variable v =

4kt

, we have:

¡1

Φ(x; t)dx =

¡1

¡v

dv = 1:

This conﬁrms the conservation of thermal energy since the initial heat sou rce is a unit

pointwise source at x = 0. The function Φ is called the Gaussian or heat kernel of a heat

equation.

Exercise 3.1. By direct calculation, verify that Φ(x; t) satisﬁes the heat equation u

= ku

Now, consider the heat problem given by the partial diﬀerential equation

= ku

; x 2(¡1; 1); t > 0;

with initial condition u(x; 0) = f(x), where f(x) is a continuous and integrable fun ction in

(¡1; 1). The follow ing theorem gives the solution to this problem:

6 Linear Second-Order Equations

Theorem 3.1. The solution to the given problem is given by the convolution

u(x; t) = f ∗Φ :=

¡1

Φ(x ¡ y; t) f(y) dy:

Proof. We ﬁrst show that u deﬁned as above satisﬁes the partial diﬀerential equation. We

have

(x; t) =

¡1

(x ¡ y; t) f (y) dy;

and u

x x

(x; t) =

¡1

(x ¡ y; t) f(y) dy:

The above representations of the partial derivatives are justiﬁed by the assumption of inte-

grability of f, th a t is ,

¡1

jf(x)jdx < 1;

and the Dominated Convergence Theorem which allows us to take the limits from outsi de

the integrals to inside the integrals. Accordingly, we have

¡ku

¡1

[Φ

(x ¡ y; t) ¡kΦ

x x

(x ¡ y; t)] f(y) dy:

By the above exercise, we have Φ

¡k Φ

x x

= 0, and this completes the claim. It remains to

show that the convolution satisﬁes the initial condition for t !0. We have

lim

t!0

u(x; t) = lim

t!0

4πkt

¡1

(x¡y)

4kt

f(y)dy = lim

t!0

4πkt

¡1

4kt

f(x ¡ y)dy:

Using the substitution v =

4kt

, we can write

lim

t!0

u(x; t) = lim

t!0

¡1

¡v

f(x ¡ 4kt

v)dv:

If we pass the limit t!0 inside the integral (which is justiﬁed by the Dominant Convergence

Theorem), we obtain

lim

t!0

u(x; t) =

¡1

¡v

lim

t!0

f(x ¡ 4kt

v)dv:

Since f is continuous, we ob tain

lim

t!0

u(x; t) =

¡1

¡v

f(x)dv =

f(x)

¡1

¡v

dv = f (x);

which completes the proof. 

Example 3.1. Consider the following problem:

= u

x x

u(x; 0) =



1 ¡1 < x < 1

0 otherwise

3.1 Heat equation 7

The initial heat p roﬁle, denoted by f(x), is 1 for x in the interval (¡1; 1), and zero elsewhere.

The ﬁgure below show s the solution at time t =0.1. Note that the solution is smooth, despite

the initial condition having jumps a t x = ±1. This is a characteristic property of the heat

equation, w he re it immediately smooths out any discontinuities after t = 0.

-4 -2 0 2 4

0.2

0.4

0.6

0.8

In this equation, we have taken k = 1. Remember that k is proportional to the diﬀusivity

of the rod. Therefore, for larger values of k, the solution u(x; t) spreads out faster, as show n

in the ﬁgure below for k = 4.

-4 -2 0 2 4

0.2

0.4

0.6

0.8

Exercise 3.2. Another interesting property of the convolution solution of the heat equation is that

the initial temperature of distant points can have an i mmediate impact on the heat evolution of a ﬁxed

po int. This is because the solution invo lves an integral over the entire domain, and the contribution

from faraway points is not necess arily negligible. This is in contrast to so me other types of diﬀerential

equations, where the evolution of a point is only inﬂuenced by its immediate neighbors. This property

of the h eat equation m akes it particularly useful for modeling systems where heat can quickly spread

throughout a medium, such as in heat transfer pr ob lems in engineering or physics.

For the simple heat equation u

= u

, suppose the initial condition is given by f (x) = δ(x) + δ(x ¡3),

that is, two point heat sources at x = 0 and x =3. Find the convolution solution and then use the following

cod e in Matlab to draw the solution at t = 0.5:

x=-10:0.01:10;

k=1;t=0.5;

Fi=@(z) exp(-z.^2/(4*k*t))/sqrt(4*pi*k*t);

plot(x,Fi(x)+Fi(x-3))

8 Linear Second-Order Equations

Exercise 3 .3. From the ab ove exercise, we learned that the heat equation is a diﬀusion equation,

which means that heat is gradually spread throughout the rod. This spreading process is aﬀected by the

diﬀusivity factor k and the initial heat distribution, but it also depends on the position of the ﬁxed p oint

relative to all other points on the rod. In other words, the heat at distant points has a n immediate impact

on the heat evolution of a ﬁxed point. In addition to its smoothing-out property, another important

feature of the heat equation is its time-invariance. This means that the solutions to the system:



= ku

t > 0

u(x; 0) = f (x)

;

and



= ku

t > t

u(x; t

) = Φ(x; t

) ∗f(x)

;

are the same. In other words, we can deﬁne a ﬂow as

(f (x)) = Φ(x; t) ∗ f(x);

that satisﬁes the property

t+s

(f(x)) = φ

(φ

(f (x))):

This means that the heat evolution of a ﬁxed point depends only on the initial heat distribution over

the entire rod a nd the diﬀusivity factor k, and not on the speciﬁc time at which we consider the system.

The ﬂ ow φ

thus provides a way to study the time-evolution of the heat equation in a systematic way, by

tracking how initial heat distributions evolve over time. We will prove this fact in subsequent chapters

using the Fourier transform. However, we can verify t he fact for the Dirac delta function, f(x)=δ(x),

here. Consider the heat equation u

= u

with initial condition u(x; 0)=δ(x). Show that the solution

of this problem at time t + s is equal to the solution of the system



= u

u(x; s) = Φ(x; s)

;

at time t > 0. In other words, the following relation holds:

Φ(x; t) ∗Φ(x; s) = Φ(x; t + s);

for t; s > 0.

Exercise 3.4. The above exercise helps us to gain a better understanding of the solution to a heat

equation as a convolution integral. Speciﬁcally, a heat problem can be viewed as a system with the

response φ

(f (x)) for the initial condition f . When f (x) = δ(x), the imp ulse response of the system is

Φ(x; t). Recall that any continuous function f(x) can be expressed as the convolution

f(x) =

¡1

f(y) δ(x ¡ y) dy:

We can rewrite this integral as a Riemann sum

f(x) ≈

f(y

) δ(x ¡ y

) δy

The response of the heat system to this function is

(f(x)) ≈

f(y

) φ

[δ(x ¡y

)] δy

f(y

) Φ(x ¡ y

; t) δy

;

Returning to the integral formulation, we obtain

(f(x)) =

¡1

f(y) Φ(x ¡y; t) dy:

In the above proof, we have used the following properties

a) The solution of the system



= u

u(x; 0) = αδ(x)

;

for a constant α is u = αΦ(x; t).

3.1 Heat equation 9

b) The solution of the system



= u

u(x; 0) = δ(x ¡y)

;

is u = Φ(x ¡ y; t).

Justify both properties through direct calculation.

Exercise 3.5. We can compare the heat kernel Φ(x; t) to the Gaussian or normal probability density

function of a random variable X with expected value µ = 0 and variance σ

, given by

N(x; σ) =

σ 2π

2σ

Upon comparison, we see that the standard deviation of the Gaussian distribution is related to time as

σ = 2t

a) Use a Z-tabl e and ﬁnd a > 0 in terms of σ such that

P (¡a < X < a) = 0.95:

b) Since 95 percent energy of Φ is concentrated in (¡a; a), we can use this value to approximate the

convolution integral Φ(x; t) ∗ f(x) as

Φ(x; t) ∗ f (x) ≈

¡a

Φ(x ¡ y; t) f(y) dy:

c) Suppose jf (x)j≤M, show the inequality



Φ(x; t ) ∗ f (x) ¡

¡a

Φ(x ¡ y; t) f (y) dy



≤0.05M :

Exercise 3.6. Verify that the function Φ

given by the function

; :::; x

; t) = (4πkt)

¡n/2

¡(x

+···+x

)/4kt

;

solves the heat equation

= k∆u;

on R

, and furthermore the convolution

u(x; t) = Φ

(x; t) ∗ f (x);

solve the heat problem



= k∆u

u(x; 0) = f (x)

;

where f is a continuous and integrable function in R

. For the proof, we need the result

; :::x

; t)dV = 1;

which we accept without prof.

Exercise 3.7. Assume th at f is a piecewise continuous function and integrable in (¡1; 1). Suppose

that lim

x!0

f(x) = a and lim

x!0

f(x) = ¡a. Show the following relation

lim

t!0

4πkt

¡1

f(x) e

4kt

dx = 0:

Exercise 3.8. Consider the equation u

= u

for x 2 (¡1; 1), and t > 0. If u(x; t) and u

(x; t) are

bo unded and furthermore lim

x!±1

u(t; x) = 0, then show the following integral is decreasing

P (t) =

¡1

ju(x; t)j

dx:

Conclude that the following problem



= u

u(0; x) = f(x)

;

10 Linear Second-Order Equations

has a unique solution under the above as sumptions.

Exercise 3.9. Consider the equation u

=div ( α ru) on R

× R

. Suppose u(x; t) and ru(x; t) are

bo unded, and furthermore lim

jxj!±1

u(x; t) = 0. Show that the following integral is decreasin g

P (t) =

¡1

ju(x; t)j

dV :

Conclude that the following problem



=div (αru)

u(x; 0) = f (x)

;

has a unique solution under the above as sumptions.

3.1.5 Solution for problems in bounded domains

The convolution solution is only applicable for heat problems deﬁned on doma ins without

boundary. To see t his, consider the simple heat problem on the interval [0; 1]:

= ku

x x

x 2(0; 1); t > 0

u(0; t) = u(1; t) = 0

u(x; 0) = f (x)

Note that the convolution integral does not satisfy the prescribed boundary conditions at

x = 0; 1. Therefore, the heat ke rnel Φ (x; t) cannot be used to obtain the solution of heat

problems deﬁned on domains with boundaries. In subsequent chapters, we will study partial

diﬀerential equations deﬁned on bounded domains, and we will need to develop new tools

that allow us to express the solution o f PDEs as inﬁnite s eries o r improper integrals .

However, we can still obtain some impo rtant properties o f the heat equati on on bounded

domains by solving boundary value problems. T he following exercises illustrate some of these

properties.

Exercise 3.10. Consider the equation u

= k u

in (0; L) where k > 0 is a constant. Show th at the

integral

P (t) =

ju(x; t)j

dx;

is decreasing for the following boundary conditions

i. u(0; t) = u(L; t) = 0

ii. u(0; t) ¡u

(0; t) = 0, u(L; t) + u

(L; t) = 0.

Exercise 3.11. Let Ω be a bounded open set with smooth boundary bnd(Ω). Consider the following

Dirichlet problem

(

=div (αru)

bnd(Ω)

For the power function

P (t) =

Ω

ju(x; t)j

dV ;

show P (t) is decreasing.

Exercise 3.12. Consider the following heat problem

= ku

¡ru

(0; t) = u

(L; t) = 0

u(x; 0) = f (x)

3.1 Heat equation 11

where r is a constant. Let E(t) be the following energy function

E(t) =

u(x; t) dx:

Prove the following relation

E(t) =



f(x) dx



¡rt

Exercise 3.13. Repeat the above problem for the following equation

=div (αru) ¡ru



bnd(Ω)

= 0

u(x; 0) = f (x)

Exercise 3.14. Let Ω be an open bounded set with smooth boundary bnd(Ω). Consider the following

heat problem on Ω

= k∆u



bnd(Ω)

= 0

u(x; 0) = f (x)

;

a) Show the following relation

Ω

u(x; t) dV =

Ω

f(x) dV :

b) Now for the problem

= k∆u ¡ru



bnd(Ω)

= 0

u(x; 0) = f (x)

;

where r is a constant, show the following relation

Ω

u(x; t) dV =



Ω

f(x) dV



¡rt

3.2 Wave equation

3.2.1 Derivation of the eq u ation

In the ﬁrst chapter, we obtained the equation of motion for a vibrating string by modeling it

as an inﬁnite number of small mass-spring systems. We now extend this approach to derive

the wave equation for an elastic membrane in two-dimensional Euclidean space. The method

we use can be generalized to any dimensi o n R

, where n is greater than or equal to one.

Let Ω be an elastic membrane in the two-dimensional space, and let D be an inﬁnitesimal

element of area δA =δx δy. For the displacement function u(x; y; t), the tension on the

boundary of D is given by T =τ ru, where τ is called the stress and depend s on the phy sical

characteristics of the elastic material. The vector force exerted on D can be de ﬁned by the

integral

F =

τ ru ·νdl;

12 Linear Second-Order Equations

where ν is the unit exterior normal vector t o the boundary of D. Applying Newton's second

law to this small portion of the membrane, we obtain

ρδAu

= F ;

where ρ is the d en sity of the membrane. By applying the divergence theorem, we can write

lim

δA!0

δA

div (τ ru) dS =

div (τ ru):

If the medium is homogeneous with co nstant τ and ρ, we can write the derived equation

as u

= c

∆u, where c=

is the wave speed of the elastic medium. This equation is

known as the wave equation f o r elastic membranes. This method can b e extended to higher-

dimensional spaces, where the wave equation takes a similar form.

The derived formulation is incomplete without specifying the initial conditions for the

equation. Note that the equation c o ntains a second-order partial derivative of u with respect

to time. Therefore, the equation should be accompanied by two initial conditions: one to

specify the initial displacement u(x; y; 0) and the seco nd to specify the initial velocity u

(x;

y; 0). These initial conditions correspond, respectively, to the initial stretch potential energy

of the membrane and the initial kinetic energy. These energies cause a dynamic for the

membrane that forces it to oscillate like a wave function. Evidently, if u(x; y; 0) = 0 and

(x; y; 0) = 0, then u = 0 for all t. Therefore, we can write the wave equation in general

dimension as follows

= c

∆u

u(x; y; 0) = f(x; y)

(x; y; 0) = g(x; y)

;

for x 2Ω ⊂R

, where f and g are initial displacement and velocity respectively.

If Ω is a bounded domain, such as a string, then the natural boundary condition for the

equation is the homogeneo us Dirichlet condition. Physically, this conditio n means that the

endpoints of th e string are fastened and ﬁxed. However, other possible boundary conditions

include the Neumann an d Robin conditions, which we have seen previously for the heat

equation. The wave equation with homogeneous Robin's boundary condition on a bounded

domain Ω is

= c

∆u on Ω

αu + β

= 0 on bnd(Ω)

u(x; y; 0) = f (x; y)

(x; y; 0) = g(x; y)

3.2.2 Damped wave problem

In the formulation of the wave equation, we neglected the drag or resistance force that the

membrane may experience. To include such a possible forces, we assume that the drag forces

is proportional to the velocity of the displacement, that is,

drag

;

3.2 Wave equation 13

and write the damped wave equation as

+ ξu

= c

∆u;

where ξ > 0 is a constant. If in addition, the membrane is under an external force h, the

equation reads

+ ξu

= c

∆u + h:

3.2.3 Solution of 1D wave equation

The one-dimensional wave equation u

= c

x x

has a simple solution. By denoting the

time derivative with the operator @

and the space derivative with @

, we can represent

the equation as (@

¡c@

)(@

+ c@

) [u] = 0 or equivalently (@

+ c@

)(@

¡c@

) [u] = 0. This

representation al lows us to decompose the wave equation into two uncoupled ﬁrst-order

PDEs:



+ cu

= 0

¡cu

= 0

;

which can be solved for arbitrary functions F (x ¡ ct) and G(x + c t). Thus, the general

solution to the equation is u(x; t) = F (x ¡ct) + G(x + ct).

Now, consider the wave equation with initial conditions u(x;0)= f (x), and u

(x;0)= g(x).

It is possible to determine F ;G such that the general solution satisﬁes these initial conditions.

To ﬁnd F and G, we can u se the initial conditions to set up a system of equations. From u(x;

0)= f(x), we have F (x)+ G(x)= f(x). Taking the derivative with respect to t a nd evaluating

at t = 0, we get u

(x; 0) = F

(x) ¡G

(x) =

g(x), and ﬁnally we obtain the following system

(

F (x) + G(x) = f(x)

(x) ¡G

(x) =

g(x)

A simple manipulation gives the following solution

u(x; t) =

[f(x ¡ct) + f(x + ct)] +

x¡ct

x+ct

g(s) ds:

The above formula is common ly known a s the D'Alembert formula for solving 1D wave equa-

tions, named after the French mathematician Jean-Baptiste Rond d'Alembert. According to

the formula, the value of u at an arbitrary point (x; t) is determined by the values of f and

g in the segment [x ¡ct; x + ct]. This segment is referred to as the wave equation's inﬂuen ce

domain, as illustrated in the ﬁgure below. It is worth noting that if g ≡0, then the value o f

u(x; t) represents the average of the right wave f(x ¡ct) and the left wave f(x + ct), bo th

propagating with the same speed c.

In the wave equation, the value of u at an arbitrary point (x; t) is determined not only by

the initial conditions but also by the values of u a nd its derivatives in the segment [x ¡ct;

x + ct]. This means that the inﬂuence of far distant po ints on a ﬁxed point is not immediate,

but rather it takes time for the values to propagate and reach the ﬁxed point, and this time

is determined by the velocity c. So the behavior of the wave equation is very diﬀerent from

that of the heat equation, where the value of u a t a point is determined by the co nvolution

over entire axis x.

14 Linear Second-Order Equations

x + ctx ¡ct x

(x; t)

η = x + ct

ξ = x ¡ct

Example 3.2. Consider the following problem

= 9u

x x

u(x; 0) = f(x)

(x; 0) = 0

;

where f(x) is the function

f(x) =



1 ¡1 ≤x ≤1

0 otherwise

The ﬁgure below illustrates the solution at time t = 0 and t = 1 for c = 3.

-6 -4 -2 0 2 4 6

0.2

0.4

0.6

0.8

As it is observed, the initial condition is split into two branches, one moving to the right

and one moving to the left with the same speed c = 3. T his is because the solution to the

wave equation is a superposition of two waves propagating in opposite directions. Each

wave carries half of the energy of the initial function, which is why we see two identical

waves moving away from the origin in the ﬁgure. While the heat equation can sm ooth out

discontinuities in the initial condition, wave equations are unable to do so, which can lead

to issues with properly deﬁning the partial diﬀerential equation at these points.

Exercise 3.15. Consider the equation

¡c

= 0:

Take ξ = x ¡ct and η = x + ct. Show that the equation under the transformation reduces to

ξη

u = 0:

Solve the obtained equation and conclude that the solution of the equation is

u(x; t) = F (x ¡ct) + G(x + ct);

for arbitrary functions F ; G.

3.2 Wave equation 15

Exercise 3.16. Consider the non-homogeneous equation

(

= c

+ h(t; x)

u(x; 0) = u

(x; 0) = 0

a) Use the characteristic change of variable ξ = x ¡ct, η = x + ct to ch an ge the equation u

= c

to the form u

ξη

= 0.

b) Use this method to show that the solution to the above non-homogeneous problem is

u(x; t) =

Z Z

∆

h(t

; x

) dt

;

where ∆ is the inside of the triangle constructed on the characteristic ξ = x ¡ct, and η = x + ct

3.2.4 Solution on ﬁnite string

We will now derive the solution to a 1D wave problem on a string o f length L w hos e endpoints

at x = 0 and x = L are ﬁxed. The problem is given by the wave equation:

with the initial and bo undary conditions:

u(t; 0)=u(t; L)=0

u(x; 0)= f(x)

(x; 0)= g(x)

To solve this equation, we can use the following formula:

u(t; x) =

odd

(x ¡ct) + f

odd

(x + ct)] +

x¡ct

x+ct

odd

(s) ds: (3.1)

Here, f

odd

and g

odd

are the odd extensions of f and g, respectively, and are 2L-periodic

functions. Recall that the odd extension of a f unction f deﬁned on [0; L] is

odd

(x) =



f(x) 0 ≤x ≤L

¡f(¡x) ¡L ≤x ≤0

The solution given by the above formula is obtained by applying the D'Alembert formula to

the odd extension of the initial condition f(x) and the odd extension of the initial velocity

g(x). The ﬁrst term in solution represents the contrib ution of the initial displacement to

the solution, and the second term repres ents the contribution of the initial velocity to the

solution.

The boundary conditions u(t; 0)=u(t; L)=0 are automatically satisﬁed by this solution

since fodd(0)= f

odd

(L)=0.

Example 3.3. Consider the following problem

= u

u(0; t) = u(4; t) = 0

u(x; 0) = f (x)

(x; 0) = 0

;

16 Linear Second-Order Equations

where

f(x) =

x ¡1 1 ≤x ≤2

3 ¡x 2 ≤x ≤3

0 otherwise

The ﬁgure below illustrates the initial condition f(x), and u(x; 1):

0 1 2 3 4

0.2

0.4

0.6

0.8

For t > 1, two branches of waves hit the boundary points x = 0 and x = 4. According to

formula (3.1), the o dd extensio n of the wave moves back after it hits the end point. The

ﬁgure below de picts the solution u(x; t) fo r t = 0; 1; 3, and 4.

0 2 4

-1

-0.5

0.5

0 2 4

-1

-0.5

0.5

0 2 4

-1

-0.5

0.5

0 2 4

-1

-0.5

0.5

Exercise 3.17. Verify the formula ( 3.1).

3.2 Wave equation 17

Exercise 3.18. Consider the following equation

= u

¡1; 0 < x < 1

u(0; t) = u(1; t) = 0

u(x; 0) =

x(x ¡1);

(x; 0) = sin(πx)

: (3.2)

Verify that the following function is the solution of the above problem:

u(x; t) =

x(x ¡1) +

sin(πx) sin(πt):

Exercise 3.19. Consider the following wave problem

= c

u(0; t) = u(L; t) = 0

u(x; 0) = f (x)

(x; 0) = g(x)

The energy of the wave function is deﬁned by the following formula

E(t) :=

[ju

(x; t)j

] dx:

a) Suppose that f

(x) is integrable. Verify the following formula which is known as the conservation

of energy of the wave function

E(t) =

fjg(x)j

(x)j

gdx:

b) Conclude that the the given problem has a unique solution.

Exercise 3.20. Let Ω ⊂ R

be a bounded set with smooth boundary bnd(Ω). Consider the following

wave problem on Ω

= c

∆u

bnd(Ω)

u(x; y; 0) = f (x; y )

(x; y; 0) = g(x; y)

;

and assume that

Ω

jg j

jrf j

dS < 1:

a) The energy of the vibrating membrane is deﬁned through the following integral

E(t): =

Ω

[ju

(x; y; t)j

jru(x; y; t)j

] dS

Prove that E(t) is constant and conclude the following relation

E(t) =

Ω

[jg j

jrf j

] dS:

b) Prove that the above problem can not have multiple solutions.

18 Linear Second-Order Equations

Exercise 3.21. Consider the following wave problem

= u

¡u

u(0; t) = u(L; t) = 0

u(x; 0) = f (x)

(x; 0) = g(x)

i. Show the following conservation of energy for the above problem

(juj

+ ju

) dx =

(jf j

+ jf

) dx:

ii. If f (0) = f(L), show that

(juj

+ ju

) dx =

jf + f

dx:

Exercise 3.22. Consider the following damped equation

+ ξu

= u

u(0; t) = u(L; t) = 0

u(x; 0) = f (x)

(x; 0) = g(x)

;

where ξk ≥0 is a constant. For E(t), the energy deﬁned as

E(t): =

(ju

+ ju

) dx;

show the following relation

¡2ξt

E(0) ≤E(t) ≤E(0):

Exercise 3.23. Consider the following wave problem

= c

∆u

bnd(Ω)

u(x; 0) = f (x)

(x; 0) = g(x)

;

where x 2 Ω ⊂ R

, and Ω is a bounded domain wit h smooth boundary bnd(Ω). Verify the following

conservation of energy formula for the problem

Ω

[ju

(x; t)j

jru(x; t)j

] dV =

Ω

[jg(x)j

jrf(x)j

] dV :

Exercise 3.24. Consider the following wave problem

= ∆u ¡u

u(t; bnd(Ω)) = 0;

u(0; p) = f(p); @

u(0; p) = 0

;

for x 2 Ω ⊂ R

, and Ω is a bounded domain with the smooth boundary b nd( Ω). Show the following

conservation of energy for the above problem

Ω

(juj

+ ju

+ jruj

) dV =

Ω

(jf j

+ jrf j

) dV :

3.2 Wave equation 19

Exercise 3.25. Consider the following damped equation

+ ξu

= ∆u

bnd(Ω)

u(x; 0) = f (x)

(x; 0) = g(x)

;

where ξ > 0 is a constant, and Ω is the unit cube [0; 1]

in R

i. For E(t), the energy deﬁned by

E(t): =

Ω

jruj

;

show that E(t) ≤E(0) for t > 0.

ii. Show the following relation

E(t) ≤E(0) e

¡ξt

3.3 Poisson and Lapl ace equations

The Po isson and Laplace equations are two important partial diﬀerential equations that

arise in many areas of science and engineering. The Poisson equation describes the steady-

state behavior of a ﬁeld, such as the electric potential or temperature, in the presence

of a given source or charge distribution. The Laplace equation is a special case of the

Poi sson equation where the source is zero, and it describes the equilibrium state of the ﬁeld.

These equations have many applications in physics, engineering, and mathematics, and their

solutions can provide valuable i nsights into the behavior of physical systems. In this chapter,

we will discuss the Poisson and Laplace equations and their solutions in one, two, and three

dimensions.

3.3.1 Derivation of Poisson an d Laplace equations

Let Ω be an open subset of R

. The Poisson's equation has the general form ¡∆u = f . If Ω

is bounded with boundary bnd(Ω), the associated boundary condition can be

αu + β



bnd(Ω)

= g;

where α and β a re constants. I n particular, if f is identically zero, the partial diﬀerential

equation is called a Laplace equation.

The Poisson equation c a n be considered as the steady-state solution of a heat problem.

For example, the steady-state solution of the non-homogeneous heat equation

= ∆u + f on Ω

αu + β

= g on bnd(Ω)

u(x; 0) = u

(x)

;

is expressed as the equation

(

¡∆u = f on Ω

αu + β

= g on bnd(Ω)

;

20 Linear Second-Order Equations

which is independent of th e initial condition u

In the electrostatic context, the Poisson and Laplace equations arise as the potential ﬁeld

of electric charges. Consider a distribution of electric cha rge with density ρ(r) in R

. The

electric ﬁeld at r 2R

generated by this distribution is determined by the following integral:

E(r) =

4π"

r~ ¡r~

jr~ ¡r~

ρ(r

) dr

; (3.3)

where "

> 0 is the permittivity constant. We can apply Gauss's law, which relates the ﬂux of

electric ﬁeld through a surface to the charge enclosed by that surface, to a bounded domain

Ω ⊂R



bnd(Ω)

E(r) ·ν dS =

;

where Q i s the total charge enclosed by bn d(Ω):

Q =

ZZZ

Ω

ρdV :

Using the divergence theorem, this can be rewritten in terms of the dive rgen ce of the electric

ﬁeld:

ZZZ

div (E) dV =

ZZZ

Ω

ρdV :

It turns out that the electric ﬁeld E is potential in the sense that there is a scalar function

φ such that E = ¡rφ. Accordingly, we obtain the following diﬀerential equation for φ:

¡∆φ(r) =

ρ(r)

Exercise 3.26. Recall that if we p ut a pointwise electric charge q at the origin, it pro duces an electric

ﬁeld E that bu the Columb law is

E(r) =

4π"

jrj

r^;

for any r 2R

¡f0g, where r^ is the unit vector in the direction o f r.

a) Ver ify that for φ = ¡

4π"

jrj

, we have E = ¡rφ.

b) For any r 2R

¡f0g, verify the relation ∆φ(r) = 0.

Exercise 3.27. f the electrical ﬁeld E(r ) is known, its associated potential φ with the relation E =¡rφ

is derived by the relation

φ(r) =

E(r

) ·T

dl;

where T

is the unit tangent vector to the straight line connecting 1 to r.

Exercise 3.28. The most well know conservative ﬁeld in physics is the gravitation ﬁeld. If a particle

with the mass m is located at the origin in R

, then its potential U (x; y; x) is

U = ¡

jrj

;

for r 2R

¡f0g. Now, suppose that in R

; n ≥3, t he potential function U has the form

U(x

; :::; x

) = ¡

+ ···+ x

)

α/2

;

3.3 Poisson and Laplace equations 21

for some constant α. Find α such that ∆U = 0.

Exercise 3.29. Despite the heat and wave equations, the Poisson and Laplace equations with Neumann

bo undary conditions can be tricky and the existence of a solution may fail. Consider the following

equation:

(

¡∆u = f on Ω

= g on bnd(Ω)

a) Use the divergence theorem to show that the necessary condition for the existence of a solution

to the equation is:

bnd(Ω)

g =

Ω

b) Even if a solution u exists for the equation, show that the solution is not unique and there are

inﬁnitely many other solutions.

c) Let B ⊂ R

be the unit disk centered at the origin. Find α such that the following equation is

solvable.

(

∆u = 1 + x



= α

;

where S = bnd(B) is the unit circle.

Exercise 3.30. The Maxwell's equations are as follows

r·E =

; r·B = 0; r×E = ¡B

; r×B = µσE + µ"

;

where E ; B are respectively the electri c and magnetic ﬁelds in space, ρ is t he density of the electric

charge, µ is the magnetic inductive capacity and " is the electric inductive capacity and σ is the electrical

conductivity. Use the above equations and derive the following equation

∆E = µ"

+ µσE

rρ

3.3.2 Fundamental solution of Poisson's equation

Consider the equation ¡∆u = f(r), where r 2R

. We can write f as a convolution integral

using the Dirac delta function δ:

f(r) =

f(r

) δ(r ¡r

) dV :

This suggests that the solution to the above equation is related to the solution of th e equation

¡∆u = δ(r). Recall that the potential function φ = ¡

4π"

jrj

describes the electric ﬁeld of a

pointwise unit charge a t the origin. We can take u(r) = Φ(r) :=

4π jrj

for r 2 R

¡f0g, and

although a rigorous proof is beyond the scope of this book, we c a n verify the following:

¡∆Φ = 0 = δ(r);

for r 2R

¡f0g. Moreover, for any ball B

of radius a centered at 0, we have:

ZZZ

δ(r) dV = 1:

22 Linear Second-Order Equations

Using the divergence theorem, we can rewrite the left-hand side of the above equation:

ZZZ

∆ΦdV = ¡



bnd(B

)

rΦ ·νdS = ¡



bnd(B

)

4π j rj

·r^dS ;

where we used the fact that ν = r^ on the surface of the ball B

. Using the relation

4π j rj

= ¡

4π jrj

;

and the equality jrj

= a

on bnd(B

), we obtain



bnd(B

)

4π j rj

·r^dS =

4πa



bnd(B

)

dS = 1:

Therefore,

ZZZ

∆Φ(r) dV =

ZZZ

δ(r) dV ;

for any a > 0. We can now obtain the solution to the equation ¡∆u = f using the convolution

integral:

u(r) =

ZZZ

f(r

) Φ(r ¡r

) dV :

Example 3.4. Consider the equation ¡∆u(r) = jrje

¡jrj

. The solution of the equation is

u(r) =

4π

ZZZ

¡jr

jr ¡r

dV :

The value of u a t the origin r = 0 is equal to

u(0) =

4π

ZZZ

¡jr

dV =

4π

lim

R!1

¡π

¡r

sin(θ) dr

dθdφ = 2:

For r =/ 0, we can use a numerical integrator to ﬁnd the value of u(r) by the above integral.

Remark 3.2. While the convolution solution for the Poisson equation works well on R

, it

does not always satisfy the desired bou ndary conditions when applied to bounded domains Ω.

In fact, ﬁnding a closed-form solution for the Poisson equation on a general b o unded domain

with arbitrary boundary conditions is typically impossible. However, in later chapters, we

will develop tools to solve this equation on speciﬁc types of domains such as rectangles, disks,

cylinders, and spheres.

Exercise 3.31. Let f = e

¡jrj

. Use a numerical integrator and ﬁnd the solution of equation ¡∆u = f

on the sphere jrj= 1.

Exercise 3.32. Let us ﬁnd the fundamental solut ion o f the Poisson equation in R

¡∆Φ = δ(r ):

3.3 Poisson and Laplace equations 23

Assume Φ = Φ(jrj), where jrj= x

+ y

. Calculate ∆Φ and obtain th e following equation for r =/ 0:

djrj

jrj

dΦ

djrj

= 0;

Solve the equation and obtain Φ(r) = C lnjpj for some constant C. Find C.

3.3.3 Harmonic functions

Harmonic functions play an important role in many areas of mathematics and physics. A

smooth function u deﬁned on an open set Ω ⊂R

is c a lle d harmonic if it satisﬁes the Laplace

equation ∆u(x) = 0 for all x 2 Ω. There are many examples of harmonic functions in R

including some simple ones that can be obtained using basic calculus. The following exercise

gives a few such examples.

Exercise 3.33. Verify that the following functions are solutions t o the Laplace equation in R

i. f (x; y) = constant

ii. f (x; y) = ax + by

iii. f(x; y) = x

¡ y

iv. f (x; y ) = cos (nx) cosh (ny).

Harmonic fu nc tions po ssess several interesting properties. One of these properties is that

the net ﬂux of the gradient vector ﬁeld F = ru through an arbitrary closed surface S is zero.

Speciﬁcally, for any closed surface S that encloses the open se t Ω, we have

0 =

Ω

∆u =

ru ·ν; (3.4)

where ν is the unit outward norma l vector to S. This property also explains why the Laplace

equation ∆u = 0 corresponds to the steady state of a heat equation u

= ∆u: if we interpret

u as the temperature f unction, then the gradient ru represents the heat ﬂux, and since

u is s teady, the net heat ﬂux through any closed surface S must be zero, otherwise the

temperature o f points inside S would change in time. This observation leads to the following

important theorem

Theorem 3.2. A harmonic function u in an open set Ω can not have a proper maximum

and minimum inside Ω.

The theorem expresses a fundamental result in the theory of harmo nic functions. It

states that a smooth function u, which is harmonic in an open set Ω, cannot have a proper

maximum or minimum inside Ω. By proper maximum or minimum, we mean a point x

and

a neighborhood D containing x

such that u(x

) > u(x) or u(x

) < u(x) for all x 2D.

To understand why this is true, consider the associated vector ﬁeld F = r u. As we saw

in t he ﬁrst chapter, proper maximum points behave like sinks for the vector ﬁeld, while

minimum points behave like so urces. In either case, the net ﬂux passing through the level

curves arou nd the m a ximum or minimum points is not zero, which violates the relation (3.4).

24 Linear Second-Order Equations

For example, the function u = x y e

¡x

¡y

can not be a harmonic function. The ﬁgure

below illustrates th e associated vector ﬁeld of the function. Observe that the net ﬂux passing

through the level curves around the maximum and minimum points are non-zero.

The following fact is an immediate result of the above argument:

Theorem 3.3. (Maximum principle) If u is harmonic function in an bounded open

set Ω ⊂R

and if u is continuous on cl(Ω), then u attains its maximum and minimum on

bnd(Ω).

Exercise 3.34. Let Ω ⊂R

be a bounded open set with smooth boundary bnd(Ω). Use the maximum

principle and prove that the unique solution of the Laplace equation



∆u = 0 on Ω

u = a on bnd(Ω)

;

for a a constant, is the constant solution u = a.

Exercise 3.35. Let Ω ⊂R

be a bounded open set with smooth boundary bnd(Ω). Use the maximum

principle and show that if there is a s olution to the following problem



∆u = f on Ω

u = g on bnd(Ω)

;

then it is unique.

Exercise 3.36. Let Ω ⊂R

be a bounded open set with smooth boundary bnd(Ω). Use the maximum

principle and show that if there is a s olution to the following problem

(

¡∆u = f on Ω

u +

= g on bnd(Ω)

;

then it is unique.

Exercise 3.37. Let Ω be an open bounded set with smooth boundary bnd(Ω). Prove that the necessary

condition that the following problem has a non-trivial sol ution is λ > 0

(

∆u = ¡λu on Ω

u +

= 0 on bnd(Ω)

Exercise 3.38. Let us write the wave equation on a bounded domain Ω as follows

= ∆u:

3.3 Poisson and Laplace equations 25

When c !1, the equation approached to the Laplace equation ∆u = 0. Therefore, we can consider the

Laplace equation as a wave pr opagation equation with the wave speed inﬁnity. This justiﬁes t ha t if the

condition on bnd(Ω) changes, this change aﬀects immediately on the value of u inside Ω.

a) Ver ify that the functio n u = 1 + "r

sin(2θ) satisﬁes the Laplace equation



∆u = 0 on B

u = 1 + " sin(2θ) on bnd(B)

;

where B is the unit disk in R

b) For " = 0, the unique solution is u = 1. For the boundary condition uj

bnd(B)

=1 + "sin(2θ), the

solution inside B changes immediately to u = 1 + "r

sin(2θ). Run the following code in Matlab

and draw the solut ion for " = 0; 0.5; 1.

[th,r]=meshgrid(-pi:pi/20:pi,0:0.2:1);

up=@(e) 1+e*r.^2.*sin(2*th);

[x,y]=pol2cart(th,r,up(0.5));

axis([-1 1 -1 1 0 2]);

subplot(2,2,1)

surf(x,y,up(0));

axis([-1 1 -1 1 0 2]);

title('$\epsilon=0$','interpreter','latex','fontsize',12);

subplot(2,2,2)

surf(x,y,up(0.5));

axis([-1 1 -1 1 0 2]);

title('$\epsilon=0.5$','interpreter','latex','fontsize',12);

subplot(2,2,3)

surf(x,y,up(1));

axis([-1 1 -1 1 0 2]);

title('$\epsilon=1$','interpreter','latex','fontsize',12)

Another consequence of the relation (3.4) is the Mean Value Property of harmonic f unc-

tion. In its rigorous form, it states that i f u is a harmonic function on an open set Ω in R

and B

(x) denotes the open ball of radius r centered at x, and S

(x) denotes the sphere

bnd(B

(x)), then for any x 2Ω and r > 0 such that B

(x) ⊆Ω, we have:

u(x) =

(x)

u(z) dS;

where jS

j denotes the volume of the ball sphere S

. In other words, the value of a harmonic

function at any point is the average of its values over any sphere cente red at that point.

In the case of harmonic functions on Ω ⊂R

, the Mean Va lue Property can be expressed

as follows: if u is a harmonic function on an open set Ω, and B

) denotes the open disk

of radius r centered at z

= (x

; y

) and contained in Ω, then

u(z

) =

2πr

)

u(z) dz;

where C

) denotes the circle of radius r centered at z

, oriented counterclock wise. In other

words, the value of a harmonic function at any point is the average of its values over any

circle centered at that point.

To see why this is true, suppose u is harmonic in an open set Ω⊂R

and consider two

circles centered at z

=(x

; y

) with radii r and r + " inside Ω. For an arbit rary point z on

, the ci rcle of radius r, consider the a ssociated point z + " r^ on C

r+ "

26 Linear Second-Order Equations

z + "r^

r + "

We can calculate the limit:

lim

"!0



r + "

r+"

)

u(z) dl ¡

)

u(z) dl



Using polar coordinates for both integrals, we ob tain the expression in the bracket as:

r + "

r+"

)

u(z) dl ¡

r(z

)

u(z) dl = "

¡π

ru(jξ j; θ) ·r^dθ;

where ξ is a point in the line connecting z to z + " r^. T herefore, we obtain:

lim

"!0



r + "

r+"

)

u(z) dl ¡

)

u(z) dl



= lim

"!0

¡π

ru(jξj; θ) ·r^dθ;

and if we pass the limit inside the integral, we obtain

lim

"!0

¡π

ru(jξj; θ) ·r^dθ =

¡π

ru(r; θ) ·r^ dθ = 0

where the last equality is justiﬁed by the relation (3.4). This calculation justiﬁes the fact that:

)

u(z) dl = 0;

and then the integral

)

u(z) dl

is a constant. In particular, if r !0, we obtain

lim

r!0

)

u(z) dl = u(z

) lim

r!0

)

dl = 2πu(z

);

and then

u(x

; y

) =

2πr

)

u(z)dl:

This is the ﬁrst version of the mean value property of ha rmonic functions in 2D.

Example 3.5. Let us verify the mean value property for the harmonic function u = x

¡y

We have

2πr

¡π

(cos

θ ¡sin

θ) dθ =

2π

¡π

cos(2θ)dθ = 0;

3.3 Poisson and Laplace equations 27

and the obtained value is equal to u(0; 0). Now, consider the arbitrary point z

: (x

; y

). The

circle C of radius r centered at z

has the following representation in the polar coordinate

x = x

+ r cosθ; y = y

+ r sinθ:

Hence,

2πr

¡y

)dl =

2πr

¡π

[(x

+ r cosθ)

¡(y

+ r sinθ)

]dθ = x

¡y

= u(z

Exercise 3 .39. Verify that the function u = x

+ y

¡3x

y ¡3xy

is harmonic and calculate the following

integral

u(x; y) dl;

where C is the circle of radius 1 centered at the point (1; 1).

Exercise 3.40. We prove another version of Mean Value Property for harmonic functions in 2D

u(x

; y

) =

πr

)

u(x; y) dxdy:

For this, use the polar coordinate a nd write

)

u(x; y) dxdy =

¡π

u(r; θ) rdrdθ:

Then, use the ﬁrst version of the mean value property and conclude the claim.

Exercise 3.41. Use the spherical coordinate and prove the ﬁrst and second versions of the Mean Value

Proper ty for harmonic functions in 3D.

3.4 Cl as siﬁcation of linear second-order PD Es

The diﬀerential equations we have explored thus far are instances of linear second-order

equations. Understanding the behav ior of their solutions is a signiﬁcant aspect of studying

these equations, both qualitatively and quantitatively. In order to analyze them eﬀectively,

it is helpful to classify linear second-order partial diﬀerential equations into three main

categories: elliptic, parabolic, and hyperbolic equations. This classiﬁcation is based on the

properties of the principal part of the d iﬀerential operator near a speciﬁc point, which in

turn determines the characteristic beh avior of the solution in the vicinity of that point.

3.4.1 Operator form, li nearity and superposition

A useful perspective on linear diﬀerential equations is to view them as operator equations.

When dealing with a real-valued function u(x) deﬁned on a domain Ω⊂R

, the partial deriv-

ative

can be seen as a mapping @

that takes functions from the domain of continuously

diﬀerentiable functions and returns continuous functions:

: C

(Ω) !C(Ω):

28 Linear Second-Order Equations

Here, C

represents the set of continuously diﬀerentiable functions, and C represents the set

of continuous functions.

By considering

as an o perator @

, we can analyze how it acts on functions and how

it re lates to other operators in the diﬀerential equation. This viewpoint allows us to study

the properties of the operator, such a s linearity, compo sition, and its eﬀect on the solutions

of the diﬀerential equation. Similarly, we can extend the notion of diﬀerenti a l operators to

higher order derivatives. For instance, the se cond partial derivative operator can be denoted

as @

and maps functions from the domain of twice continuously diﬀerentiable functions:

: C

(Ω) !C(Ω):

Furthermore, we can deﬁne a complete linear second-order diﬀerential operator as:

L :=

i;j =1

(x) @

j=1

(x) @

+ c(x); (3.5)

where a

i; j

(x); b

(x) and c(x) are some continuous or smooth functions deﬁned on the d omain

Ω. This general form allows us to represent a wide range of second-order linear diﬀerential

equations.

Recall from linear algebra that a mapping T from R

to R

is called linear if it satisﬁes

the following relation for any α

; α

2R and u

; u

2U:

T (α

+ α

) = α

T (v~

) + α

T (v~

Similarly, we can apply this concept to the diﬀerential operator L. The operator L is linear

because it satisﬁes the following property:

L(α

(x) + α

(x)) = α

T (u

(x)) + α

T (u

(x));

where u

(x) and u

(x) are functions in the domain of the operator L, and α

and α

are

constants.

Although we wo n't delve into the full vector space formulation of partial diﬀerential

operators in this book, representing partial de rivatives in o perator notation can help us

understand the classiﬁcation of second-order linear PDEs. In general, a linear second-order

partial diﬀerential equation can be written as:

L[u]= f ;

where L is the ope rator deﬁned in (3.5). Here, L represents a linear mapping that takes a

function u as input and produces the corresponding function f as output. The goal of such

an equation is to determine the function u, given a known function f.

Interpreting a linear partial diﬀerential equation as an operator equation has several

advantages. I t allows us to make a connection between linear algebra, which studies the

properties of linear mappings in ﬁnite-dimensional vector spaces, and diﬀerential equations.

One beneﬁt of this approach is the superposition principle, which applies to homogeneous

equations of the form L[u] = 0.

3.4 Classification of linear second-order PD Es 29

Proposition 3.1. If u

; :::; u

are homogeneous solutions of linear equation L[u] = 0, then

the function

u =

j=1

;

for arbitrary c

2 R solves the equation. Moreover, if u

is a particular solution for the

equation L[u] = f, then the function

u = u

j=1

;

solves the equation if u

; ::: ; u

are in the kernel or null space of L.

Exercise 3.42. Let L be the linear operator L = @

+ @

. Consider the equation L[u] = x + y.

a) Find the general homogeneous solution of the equation.

b) Verify that u

= xy is a particular solution of the equation. Therefore, u = u

+ xy is the general

solution of the equation. It is simply veriﬁed that u

is also a particular solution.

Show that the latter particular solution can be derived from the general solution u = u

+ xy.

Exercise 3.43. Find the gener al solution of the linear equation L[u] = x for u = u(x; y), where

L = @

¡c@

for a constant c.

Exercise 3.44. Consider the following equation

¡c

= 0;

and let L be the oper ator L : = @

¡c

. Show that

L = L

◦L

;

where L

:= @

¡c@

, L

: =@

+ c@

. With the aid of the the decomposition, ﬁnd two linearly independent

solution of the given equation.

By taking some precautions, we can also extend the superposition principle to the inﬁnite

or continuous versions. For example, if u

; u

; ::: are in the null space of L, then L[u] = 0 for

the function

u =

j=1

Similarly, if u(x; k) solves the equation L[u] = 0 for k in some interval I, then the following

function also solves the equation:

u =

c(k) u(x; k) dk:

Example 3.6. Consider the partial diﬀerential equation ∆u =0 deﬁned on the half-plane x >

0. An example of a solution to this equation is give n by the functions u

(x; y) = e

¡kx

sin(ky)

for any k 2(¡1; 1). Now, suppose we deﬁne c(k) as follows:

c(k) =



0 k < 0

1 k > 0

Then, the function

u(x; y) :=

¡1

c(k) u

(x; y) dk =

+ y

;

30 Linear Second-Order Equations

solves the equation ∆u =0 on x > 0. However, as x approach es 0, the function u(x; y) becomes

unbounded and cannot be extended to the boundary of the domain. To remedy this, we can

choose c(k) as

c(k) =

(

0 k < 0

¡k

k > 0

This transforms the superposition integral to

u(x; y) =

(x + 1)

+ y

This solution can be extended to the boundary x = 0 in a proper manner.

Remark 3.3. The version o f inﬁnite superposition should be used with caution. For instance,

let us consider the equation L =

+ 1 and the functions u

(x) =

cos(nx) for any

o dd n, which solve the equation. The f unction

u =

n=1;n:odd

(x) = cos(x) +

cos(2x) +

cos(3x) + ···+

converges tof(x)=

jxj for x 2(¡π; π). However, this function is discontinuous and does

not have a derivative at x = 0. Therefore, we cannot deﬁne L [u ] as equal to zero at x = 0,

and we must be careful when using the inﬁnite version of superposition.

Exercise 3.45. Find a solution to equation ∆u = 0 on the half-plane x > 0 such that u on the boundary

x = 0 is equal to

4 + y

Exercise 3.46. Consider linear equation @

u = xy, for u = u(x; y).

a) Show that u

y solves the equation.

b) Find the homogeneous solution of the equation and write down the general solution.

Exercise 3.47. Find the gener al s olution of the following equation

L[u] = xy;

for a function u = u(x; y), where L = @

Exercise 3.48. We know that the linear equations enjoy superposition principle. This property does

not necessarily hold for non-linear equations. Verify that functions u

= x + c

and u

= y + c

are the

solutions to the nonlinear equation

+ ju

= 1;

however, u

+ u

is not a solution.

3.4.2 Classiﬁcation: c onstant coeﬃcient case

Linear seco nd-order partial diﬀerential equations can be classiﬁed based on their principal

part, which is deﬁned by the sum:

i;j

(x) @

3.4 Classification of linear second-order PD Es 31

This classiﬁcation is based on the idea that the behavior of a PDE near a given poi nt is

determined by the highest-order derivatives app earing in the equation at that point. To

illustrate the logic behind the classiﬁcation, we ﬁrst consider the case where a

are constants,

and the domain of the problem is open subsets of the plane.

It is simply seen that L

can be rewritten in the matrix form as follows:

:= (

) A





;

where A is the symmetric coeﬃcient matrix given by

A =

+ a

From linear algebra, the eigenvalues of A are real. Let λ

and λ

be the eigenvalues of A.

Based on these eigenvalues, we can classify the PDE as follows:

• Elliptic: λ

>0.

◦ Solutions to elliptic PDEs are smooth and have no s ingu larities. This type of

PDE typically arises in problems involving steady-state phenomena, such as

in electrostatics or ﬂu id mechanics. The most familiar example is the Poisson

equation ¡∆u = f, where ∆ operator can be written as

∆ = (

)



1 0

0 1





;

with the coeﬃcient matrix A =



1 0

0 1



having eigenvalues λ

= λ

= 1.

• Parabolic: λ

=0.

◦ Solutions to parabolic PDEs have a smooth ing eﬀect that can regularize initial

data. This type o f PDE often arises in problems involving diﬀusion, heat

ﬂow, or time-dependent phenomena that evolve in a certain direction. The

most fa miliar equation is the heat equation u

= ku

in the plane (x; t). The

coeﬃcient matrix is

A =



0 0

0 k



;

withe eigenvalues λ

= 0, λ

= k.

• Hyperbolic: λ

<0.

◦ Solutions to hyperbolic PDEs have wave-like behavior and can exhibit shock

waves and other d iscontinuities. This type of PDE typically arises in problems

involving wave propagation or vibrations. The most familiar example is a wave

equation u

= c

in the (x; t)-plane. The coeﬃcient matrix of this equation

A =

¡1 0

0 c

;

32 Linear Second-Order Equations

with eigenvalues λ

= ¡1, λ

= c

Remark 3.4. It is known from linear algebra that λ

= det(A), and thus the classiﬁcation

can be done based on the sign of det(A) too. Therefore, we can classify the P DE based on

the sign of det(A), as follows:

• Elliptic: det(A)>0

• Parabolic: det(A)=0

• Hyperbolic: det(A)<0

Example 3.7. Consider the following equation

x x

¡αu

+ 2u

+ 3u

= u + x:

We classify the type of the equation based on the value of α. The coeﬃcient matrix of the

principal part of the equation is

A =

We have det(A) = 4 ¡

, and thus the equation is elliptic if α

< 16 or ¡4 < α < 4. The

equation is hyperbolic if α > 4 or α < ¡4. The equation is parabolic if α = 4 or ¡4.

Exercise 3.49. Classify following equations

a) u

+ u

+ 2u

+ u

= 1 ¡u

b) u

+ 2u

¡3u

= u.

Exercise 3.50. Determine the type of the following equation based on the values of α

+ u

+ 2αu

= 0:

Deﬁnition 3.1. Let L

be the principal part of a linear second-order partial diﬀerential

equation (PDE), given by

i;j=1

;

where a

are constant coeﬃcients. The associated coeﬃcient matrix of the operator is deﬁned

by A = [a~

] where a~

+ a

. Let λ

; :::; λ

be the eigenvalues of A.

• The PDE is called elliptic if all eigenvalues λ

; :::; λ

have the same sign.

• The PDE is called hyperbolic if one eigenvalue has the opposite sign of the other

eigenvalues.

• The PDE is called parabolic if one eigenvalue is zero, and the other eigenvalues have

the same sign.

Exercise 3.51. For function u = u(x; y; t), classify the equations

a) u

= k∆u,

3.4 Classification of linear second-order PD Es 33

b) u

= c

∆u.

c) u

+ u

+ 3u

= 0.

Exercise 3.52. Find α such that the following equation is parabolic in the coordinate (t; x; y):

= 2u

+ αu

+ 2u

3.4.3 Canonic al form

For the principal part L

, the terms @

for i =/ j are called multiplicative terms. We can

remove these terms by performing a lin ear transformation of the co ordinate. Before intro-

ducing the transformation, let's explain the terminology of ellipticity, parabolicity, and

hyberbolicity by comparing them to the types of conic sections.

Recall that the equation of a conic section is given by:

+ 2bxy + cy

+ dx + ey = const:

The principal part of this section is ax

+ 2bxy + cy

, which ca n be represented in matrix

form as

+ 2bxy + cy

= (

x y

) A





;

where th e c oeﬃcient matrix A is symmetric and has real eigenvalues λ

; λ

, and two orthog-

onal eigenvectors. Let Q = [v~

jv~

] be the eigenvector matrix of A with jv~

j= jv~

j= 1.

By linear algebra, we know that:

¡1

AQ =



0 λ



Thus, the transformation





= Q





;

maps the old coordinate system to the new coordinate system (X ; Y ), and transforms the

conic section to the following form:

+ λ

= ﬁrst and zero order terms:

If λ

and λ

have the same sign, the conic section is an ellipse. If one of λ

or λ

is zero,

the conic section is a parabola. If λ

and λ

have opposite signs, then the conic section is a

hyperbola.

This argument shows how we can remove the multiplicative terms and why the termi-

nology is justiﬁed based on geometrical observation.

Example 3.8. Consider the equation

+ 2y

¡2xy = 1:

In the matrix form, we can represent the equation as

(

x y

)



2 ¡1

¡1 2





= 1:

34 Linear Second-Order Equations

The coeﬃ cient matrix A =



2 ¡1

¡1 2



has two eigenvalues λ = 1; 3 with associated eigenvectors

, v~

. Consider the coordinate transformation









; (3.6)

where Q is the matrix of eigenvectors. By performing this transformation to the equation,

we obtain

(

x y

)



2 ¡1

¡1 2





= (

X Y



2 ¡1

¡1 2







= (

X Y

)



1 0

0 3





Thus, th e given equation in the new coo rdinate system has the representation X

+ 3Y

= 1

which is a ellipse. Note that Q is a 45 degree rotation m a trix, and thus the coord inate change

is just a rotation of the old coordinate system.

-0.5 0.5

-0.5

0.5

-1 -0.5 0.5 1

We can employ the above technique for linear partial diﬀerential equation.

Example 3.9. Let us rewrite the following equation in the canonical form

+ 2u

¡2u

+ 3u

= u + x:

The coeﬃcient matrix of the principal part is A =



2 ¡1

¡1 2



, with the eigenvalues λ

= 1,

and λ

= 3. The eigenvectors of the coeﬃcient matrix A are v~





, v~



¡1



, and

thus Q =



1 ¡1

1 1



is the transformation matrix. Deﬁne the new coordinate system @

; @

through the following transformation





= Q





: (3.7)

By performing this transformation to the principal part of the equation, we reach

x x

+ 2u

¡2u

x y

= u

+ 3u

We need also to transform the ﬁrst and zero terms of the equation as well. From (3.7), we

have

;

3.4 Classification of linear second-order PD Es 35

and x =

X ¡

Y . Finally, the given PDE reduces to the following one in the new

coordinate system

+ 3u

3 2

= u +

X ¡

Y :

Exercise 3.53. Wr ite down the canonical form of the following equation

+ u

+ 4u

+ 2

= u

Exercise 3.54. Find values of α s uch that the equation

+ 3u

+ 2αu

+ 2u

¡3u

= 0;

is elliptic, parabo lic or hyperbolic. For α = 1 write down the equation in the new coordinate and

determine the associated linear transformation.

Exercise 3.55. Write the following equations in the canonical form and classify them into elliptic,

parabolic and hyperbolic equations

a) u

= 0

b) u

+ 2u

= u

c) u

¡u

¡2u

= u

+ u

d) 3u

+ 4u

+ 6u

= u

¡u

e) u

¡4u

+ 4u

¡u

= 1

f) ¡2u

¡6u

+ 6u

¡u

= 0

g) u

+ 4u

+ u

+ 2

= u

Exercise 3.56. Consider the following equation

+ u

+ 4u

= 2

a) Determine the type of equation in terms of elliptic, parabolic o r hyperbolic equation. Your answer

should be based on a calculation.

b) Use a change of coordinate from (t; x) to ( T ; X) and rewrite the equation in the normal form

without any multiplicative partial diﬀerentiation.

c) For the new equation in the coordinate (T ; X), ﬁnd a solution in the form u(T ; X) = αX

+ βT

for suitable constants α; a; β; b.

Exercise 3.57. Consider the following equation

+ 3u

+ 2 2

= 0

a) Determine the type of the equation in terms of ellipticity, parabolicity and hyberbolicity.

b) Write the canonical form of the equation with the aid of a transformation of the type









;

for an appropriate matrix Q.

c) In the new coordinate (X ; Y ), use the transformation ξ = Y + i2X, η = Y ¡i2X, and derive the

following equation in (ξ; η)

ξη

= 0:

Find the general solution of the equation in ( X ; Y ).

36 Linear Second-Order Equations

3.4.4 Classiﬁcation: variable coeﬃcients

The classiﬁcation of linear second-order PDEs with variable coeﬃcients is similar to that of

equations with constant coeﬃcients. However, for variable c oeﬃcient equations, the classi-

ﬁcation is loca l, meaning it holds only in an open neighborhood of the focal point, and not

for the entire domain. The following example illustrates this fact.

Example 3.10. Consider the diﬀerential operator

L := (1 + x

y)@

¡2xy@

+ y@

The coeﬃcient matrix of L

A(x; y) =

1 + x

y ¡xy

¡xy y

Note that det(A)= y. Thu s, L is elliptic in the half plane y >0 and hype rbolic in the half

plane y <0. For a ny point on th e x-axis, the e quation exhibits degenerate behavior.

Exercise 3.58. Consider the diﬀer ential operator

L = x@

+ xy @

+ y@

Rewrite L in the form

L =div ([a

(x; y)] r) + ﬁr st term derivatives

where [a

] is a symmetric matri x.

Problem 3. 1 . Find domains in which the following operators can be classiﬁed into elliptic, parabolic

or hyperbolic typ e

i. L := y@

¡y@

+ 2 (x ¡1)@

ii. L := cos(xy)@

+ cos(xy)@

+ 2 sin(xy)@

iii. L := x@

+ y@

+ 2xy@

iv. L := x

+ y

+ x

+ y

+ 2@

3.4 Classification of linear second-order PD Es 37