1

1 Partial Differential Equations (PDEs)

PDEs enable modelling physical phenomena. The main interset in this

course is obtaining explicit soln’s to PDEs which are very limited only
to special cases. The other alternative of obtaining numerical soln’s is
out of the scope of this course.

Another specific interest is the 2nd order PDEs that more often arises in
modelling physical phenomena. Some examples are:

Laplace eqn.: 𝑢𝑥𝑥 + 𝑢𝑦𝑦 = 0 ⇒ 𝑢 = 𝑢(𝑥, 𝑦)

Heat eqn.: 𝑢𝑡 − 𝜅𝑢𝑥𝑥 = 0 ⇒ 𝑢 = 𝑢(𝑥, 𝑡)

Wave eqn.: 𝑢𝑡𝑡 − 𝑐 2 𝑢𝑥𝑥 = 0 ⇒ 𝑢 = 𝑢(𝑥, 𝑡)

or in higher dimension

𝑢𝑥𝑥 + 𝑢𝑦𝑦 + 𝑢𝑧𝑧 = 0,

𝑢𝑡 − 𝜅(𝑢𝑥𝑥 + 𝑢𝑦𝑦 + 𝑢𝑧𝑧 ) = 0,

𝑢𝑡𝑡 − 𝑐 2 (𝑢𝑥𝑥 + 𝑢𝑦𝑦 + 𝑢𝑧𝑧 ) = 0,

respectively, or the (spatial) Laplace operator

ℒ[𝑢] ≡ 𝑢𝑥𝑥 + 𝑢𝑦𝑦 + 𝑢𝑧𝑧

in cartesian coordinates (𝑥, 𝑦, 𝑧) or

1 1
𝑢𝑟𝑟 + 𝑟 𝑢𝑟 + 𝑟 2 𝑢𝜃𝜃 + 𝑢𝑧𝑧 = 0,

in cylindrical coordinates (𝑟, 𝜃, 𝑧) or

 2

2 1 1
𝑢𝑟𝑟 + 𝑟 𝑢𝑟 + 𝑟 2 𝑠𝑖𝑛𝜙 (𝑠𝑖𝑛𝜙 𝑢𝜙 )𝜙 + 𝑟 2 𝑠𝑖𝑛2 𝜙 𝑢𝜃𝜃 = 0,

in spherical (𝑟, 𝜃, 𝜙) coordinates.

2nd Order Linear PDEs

A 2nd order linear PDE in two variables 𝒙 = (𝑥, 𝑦) can be written in

compact form as

ℒ[𝑢] = 𝑔(𝑥, 𝑦)

where ℒ is the differential operator defined by

ℒ [𝑢] ≡ 𝑎𝑢𝑥𝑥 + 2𝑏𝑢𝑥𝑦 + 𝑐𝑢𝑦𝑦 + 𝑑𝑢𝑥 + 𝑒𝑢𝑦 + 𝑓𝑢.

Here, the given coefficients (data) 𝑎, 𝑏, 𝑐, 𝑑, 𝑒, 𝑓 may be functions of

(𝑥, 𝑦). Note that the variable 𝑦 may be replaced by the variable 𝑡. The
PDE is said to be homogeneous if 𝑔 ≡ 0, otherwise nonhomogeneous.

The linearity of ℒ can be simply expressed as

ℒ [𝛼𝑢 + 𝛽𝑣] = 𝛼ℒ [𝑢] + 𝛽ℒ[𝑣]

for a pair of functions 𝑢, 𝑣 and constants 𝛼, 𝛽.

For a unique solution, boundary (and/or initial) conditions must be

specified as many as the degree of the highest derivative in space (and/or
time), respectively. The combination of a PDE with initial and boundary
conditions is referred to as an Initial-Boundary Value Problem (IBVP).
 3

Considering a 2nd order PDE on a spatial domain 𝛺,

the possible types of boundary conditions on the boundary 𝜕𝛺 are:

(i) Dirichlet Type: the value of u is specified

𝐵[𝑢] ≡ 𝑢(𝒙) = ℎ𝐷 (𝒙) for 𝒙 ∈ 𝜕𝛺𝐷 .

(ii) Neumann Type: the value of the (outward) normal derivative of 𝑢 is

specified
𝜕𝑢
𝐵[𝑢] ≡ 𝜕𝑛 (𝒙) = ℎ𝑁 (𝒙) for 𝒙 ∈ 𝜕𝛺𝑁 .

(iii) Robin Type: a linear combination of the value of u and its

(outward) normal derivative is specified

𝜕𝑢
𝐵[𝑢] ≡ 𝛼𝑢(𝒙) + 𝛽 𝜕𝑛 (𝒙) = ℎ𝑅 (𝒙) for 𝒙 ∈ 𝜕𝛺𝑅 .

(iv) Mixed Type: the three types may be applied on nonoverlapping

partitioning of 𝜕𝛺 = 𝜕𝛺𝐷 ∪ 𝜕𝛺𝑁 ∪ 𝜕𝛺𝑅 .

 4

Here, the operator 𝐵 stands for the boundary operator that has to be linear
together with the differential operator 𝐿 for the linearity of the IBVP.

e.g. An IBVP for the heat equation is completely specified, say, by the
nonhomogeneous differential operator

ℒ [𝑢(𝑥, 𝑡)] ≡ 𝑢𝑡 − 𝜅𝑢𝑥𝑥 = 𝑔(𝑥, 𝑡), (𝑥, 𝑡) ∈ (0, 𝐿) × (0, 𝑇),

by Dirichlet-Neumann boundary conditions

𝑢(0, 𝑡)
𝑩[𝑢] ≡ { } = 𝒉(𝑡), 𝑡 > 0,
𝑢𝑥 (𝐿, 𝑡)

and by the initial condition

𝑢(𝑥, 0) = 𝑠(𝑥 ), 𝑥 ∈ (0, 𝐿).

Then, let

ℒ [𝑢 ], (𝑥, 𝑡) ∈ (0, 𝐿) × (0, 𝑇) 𝑔(𝑥, 𝑡)

𝜞[𝑢] ≡ { 𝑩[𝑢], 𝑡>0 } = { 𝒉(𝑡) } ≡ 𝑭(𝑥, 𝑡)
𝑢(𝑥, 0), 𝑥 ∈ (0, 𝐿) 𝑠(𝑥 )

so that IBVP becomes

𝜞[𝑢(𝑥, 𝑡)] = 𝑭(𝑥, 𝑡).

Note that an IBVP having a unique solution that varies continuously

with the initial and boundary data is said to be well-posed. Otherwise
ill-posed. i.e. the small change 𝛿𝑭 in the data 𝑭 induces small change
𝛿𝑢 in the solution u:
 5

𝜞[𝑢 + 𝛿𝑢] = 𝑭 + 𝛿𝑭

in the sense of some function norms ‖ ∙ ‖𝑎 and ‖ ∙ ‖𝑏

‖𝛿𝑢‖𝑎 ≤ 𝐶 ‖𝛿𝑭‖𝑏 .

PDE Models

Classical PDEs are derived from physical laws with some modelling
assumptions. The reality of these assumptions determines the success of
the model in predicting the physical phenomena.

Wave Equation for a string

Consider a string with line density of 𝜌 stretched with a tension T.

Assume small angles 𝛼∓ such that 𝑐𝑜𝑠 𝛼∓ ≈ 1 and 𝑠𝑖𝑛 𝛼∓ ≈ tan 𝛼∓ .

Then, balancing horizontal motion

𝑇+ cos 𝛼+ = 𝑇− cos 𝛼− ⇒ 𝑇+ = 𝑇− = 𝑇

while using Newton’s 2nd Law in the vertical

 6

𝜌 𝛿𝑠 𝑢𝑡𝑡 (𝑥, 𝑡) = 𝑇+ sin 𝛼+ − 𝑇− sin 𝛼−

= 𝑇 (tan 𝛼+ − tan 𝛼− )
1 1
= 𝑇 (𝑢𝑥 (𝑥 + 2 𝛿𝑥 , 𝑡) − 𝑢𝑥 (𝑥 − 2 𝛿𝑥 , 𝑡))

= 𝑇 𝛿𝑥 𝑢𝑥𝑥 (𝑥, 𝑡) + 𝑂 (𝛿𝑥 2 )

leads to the wave equation

𝑢𝑡𝑡 = 𝑐 2 𝑢𝑥𝑥 as 𝛿𝑥 → 0

where 𝛿𝑥 ≈ 𝛿𝑠 and 𝑐 = √𝑇⁄𝜌.

As boundary conditions associated with the wave eqn.

(i) Dirichlet: 𝑢 may be specified as physical displacement of the end

of the string,
(ii) Neumann: Tensile force 𝑇𝑢𝑥 may be specified at the end of the
string,

(iii) Robin: Some kind of elastic attachement at the end of the string

𝑇𝑢𝑥 = 𝑘𝑢 where 𝑘 is an elasticity constant.

Conservation Law

Consider a quantity 𝑄(𝒙, 𝑡) per unit volume in a volume 𝑉 (≡ 𝛺)

bounded by a surface 𝑆 (≡ 𝜕𝛺). Then the rate of change of Q in 𝑉 is
𝑑
∫ 𝑄 𝑑𝑉 = ∫ 𝐹 𝑑𝑉 − ∫𝒒 ⋅ 𝒏 𝑑𝑆
𝑑𝑡 𝑉 𝑉 𝑆
 7

where F is the net production and 𝒒 ⋅ 𝒏 is the flux along the outward
normal 𝒏 to 𝑆. Using divergence theorem and arbitrariness of 𝑉, it
yields
𝜕𝑄
+ 𝛻 ⋅ 𝒒 = 𝐹.
𝜕𝑡
Heat equation

Let 𝑄 = 𝜌𝑐𝑇 denote heat per unit volume and 𝒒 = −𝜆𝛻𝑇 by Fourier
Law. Then
𝜕𝑇 𝜕𝑇
𝜌𝑐 + 𝜵 ⋅ (−𝜆𝜵𝑇) = 𝐹 ⇒ − 𝜅 𝛻 2𝑇 = 𝑓
𝜕𝑡 𝜕𝑡
where 𝜅 = 𝜆/𝜌𝑐 is thermal diffusivity and 𝑓 = 𝐹/𝜌𝑐 with c the specific
heat, 𝜆 thermal conductivity coefficient. In 1D and without the source
term, it becomes

𝜕𝑇 𝜕2𝑇
−𝜅 = 0.
𝜕𝑡 𝜕𝑥 2
As boundary conditions associated with the heat eqn.

(i) Dirichlet: 𝑇 may be specified on 𝑆, i.e. 𝑇|𝑆 = 𝑇0

𝜕𝑇
(ii) Neumann: 𝑆 may be insulated, i.e. 𝜕𝑛| = 0.
𝑆

(iii) Robin: Heat may be lost to the surrounding at 𝑇0 , i.e.

𝜕𝑇
| = −ℎ(𝑇|𝑆 − 𝑇0 ).
𝜕𝑛 𝑆
 8

Classification of PDEs

The classification of PDEs as Elliptic, Parabolic or Hyperbolic can be

obtained from a study of second-order linear PDEs.

Consider the general form of a 2nd order linear constant coefficients

homogeneous PDE of two independent variables x and y (or t):
Principal Part
⏞
𝑎𝑢𝑥𝑥 + 2𝑏𝑢𝑥𝑦 + 𝑐𝑢𝑦𝑦 + 𝑑𝑢𝑥 + 𝑒𝑢𝑦 + 𝑓𝑢 = 0

where 𝑎, 𝑏, … , 𝑓 are given (variable) coefficients.

The type of this equation is defined by its principal part which consists
of the highest-order terms, written in matrix form

𝜕𝑦 ) (𝑎 𝑏 𝜕𝑥
𝑎𝜕𝑥𝑥 𝑢 + 2𝑏𝜕𝑥𝑦 𝑢 + 𝑐𝜕𝑦𝑦 𝑢 = (𝜕𝑥 ) ( ) 𝑢,
⏟𝑏 𝑐 𝜕𝑦
𝚫

where the coefficients 𝑎, 𝑏, 𝑐 are frozen at some (𝑥, 𝑦) value. The

classification of PDEs is based on the relative sign of the eigenvalues of
the coefficient matrix 𝚫 as follows:

Recall that the determinant of the coefficient matrix 𝚫 (discriminant)

gives the product of the eigenvalues:

𝜆1 𝜆2 = 𝑑𝑒𝑡 𝚫 = 𝑎𝑐 − 𝑏2 .

The sign of 𝑎𝑐 − 𝑏2 determines the type of the PDE:

𝑎𝑐 − 𝑏2 > 0 : Elliptic, i.e. eigenvalues have the same sign.

 9

𝑎𝑐 − 𝑏2 = 0 : Parabolic, i.e. at least one eigenvalue is zero.

𝑎𝑐 − 𝑏2 < 0 : Hyperbolic, i.e. eigenvalues have the opposite sign.

e.g. These are the canonical forms of each type of PDEs:

[𝜕𝑡 𝜕𝑥 ] 0 0 𝜕𝑡
ℑ[𝑢] ≡ 𝜕𝑡 𝑢 − 𝜕𝑥𝑥 𝑢 ⇒ [ ][ ] : Parabolic
0 1 𝜕𝑥

[𝜕𝑡 𝜕𝑥 ] 1 0 𝜕𝑡
ℑ[𝑢] ≡ 𝜕𝑡𝑡 𝑢 − 𝜕𝑥𝑥 𝑢 ⇒ [ ] [ ] : Hyperbolic
0 −1 𝜕𝑥

[𝜕𝑥 𝜕𝑦 ] 1 0 𝜕𝑥
ℑ[𝑢] ≡ 𝜕𝑥𝑥 𝑢 + 𝜕𝑦𝑦 𝑢 ⇒ [ ][ ] : Elliptic
0 1 𝜕𝑦

Note: If the coefficients 𝑎, 𝑏, 𝑐 depend on the variables 𝑥, 𝑦 (or 𝑡), then

the PDE may be of different type at different points in the domain
evaluated by freezing the coefficients at single point (𝑥, 𝑦) (or 𝑡).

e.g. Consider the Tricomi operator

𝐿[𝑢] ≡ 𝜕𝑥𝑥 𝑢 + 𝑥 𝜕𝑦𝑦 𝑢,

whose discriminant is 𝑑𝑒𝑡 𝚫 = 𝑥. Hence, 𝐿 is elliptic in the half-plane

𝑥 > 0, hyperbolic in the half-plane 𝑥 < 0, and parabolic on the 𝑦-axis.
 10

Some Further Notes

 The classification scheme does not depend on the coordinate system

in which the PDE is expressed: Consider a change of independent
variables:

𝜉 = 𝜉(𝑥, 𝑦), 𝜂 = 𝜂(𝑥, 𝑦).

The transformation is nonsingular if the Jacobian determinant of the

transformation is nonzero, i.e.
𝜕(𝜉, 𝜂)
𝐽=| | ≠ 0.
𝜕(𝑥, 𝑦)

Let us denote the transformed dependent variable as

𝜔(𝜉, 𝜂) = 𝑢(𝑥, 𝑦)

then the PDE becomes

𝛼𝜔𝜉𝜉 + 2𝛽𝜔𝜉𝜂 + 𝛾𝜔𝜂𝜂 + 𝛿𝜔𝜉 + 𝜀𝜔𝜂 + 𝑓𝜔 = 0

where

𝛼 = 𝑎𝜉𝑥2 + 2𝑏𝜉𝑥 𝜉𝑦 + 𝑐𝜉𝑦2 ,

𝛽 = 𝑎𝜉𝑥 𝜂𝑥 + 𝑏(𝜉𝑥 𝜂𝑦 + 𝜉𝑦 𝜂𝑥 ) + 𝑐𝜉𝑦 𝜂𝑦 ,

𝛾 = 𝑎𝜂𝑥2 + 2𝑏𝜂𝑥 𝜂𝑦 + 𝑐𝜂𝑦2 ,
𝛿 = 𝑑𝜉𝑥 + 𝑒𝜉𝑦 , 𝜀 = 𝑑𝜂𝑥 + 𝑒𝜂𝑦 .

The determinant of the coefficient matrix 𝚫 is now:

𝜆1 𝜆2 = 𝑑𝑒𝑡 𝚫 = 𝛼𝛾 − 𝛽 2 = (𝑎𝑐 − 𝑏2 )𝐽2 .

 11

The sign of 𝛼𝛾 − 𝛽 2 is the same as that of 𝑎𝑐 − 𝑏2 . Hence, the

classification of PDEs is invariant under a change of coordinates.

e.g. It can be shown that the wave equation (hyperbolic)

𝑢𝑡𝑡 − 𝑐 2 𝑢𝑥𝑥 = 0, − ∞ < 𝑥 < ∞, 𝑡 > 0,

under the coordinate transformation

𝜉 = 𝑥 + 𝑐𝑡, 𝜂 = 𝑥 − 𝑐𝑡,

reduces to

𝜔𝜉𝜂 = 0,

that is still of an hyperbolic type (𝛼 = 𝛾 = 0, 𝛽 = 1⁄2).

The last form can actually be solved to get the general solution

𝜔(𝜉, 𝜂) = 𝐹 (𝜉 ) + 𝐺(𝜂),

or

𝜔(𝜉, 𝜂) = 𝑢(𝑥, 𝑡) = 𝐹 (𝑥 + 𝑐𝑡) + 𝐺 (𝑥 − 𝑐𝑡).

This solution represents two waves propagating without change of

shape with velocity 𝑐 in opposite directions along the 𝑥-axis.
 12

 In the general case of n independent variables 𝑥1 , 𝑥2 , . . . , 𝑥𝑛 ,

consider the 2nd order linear PDE
Principal Part
⏞𝑛 𝑛 𝑛
𝜕2𝑢 𝜕𝑢
∑ ∑ 𝑎𝑖𝑗 + ∑ 𝑏𝑖 + 𝑐𝑢 + 𝑑 = 0.
𝜕𝑥𝑖 𝜕𝑥𝑗 𝜕𝑥𝑖
𝑖=1 𝑗=1 𝑖=1

This equation can not be reduced to a simple canonical form over a

full region unless the coeff. 𝑎𝑖𝑗 are constants or we restrict to a
single point in (𝑥1 , . . . , 𝑥𝑛 )-space (freezing the coefficients).

The classification is then based on the eigenvalues 𝜆1 , . . . , 𝜆𝑛 of the

coefficients matrix 𝐀 = [𝑎𝑖𝑗 ]:

Def.: The n-dimensional PDE is called:

 Elliptic if all the eigenvalues 𝜆1 , . . . , 𝜆𝑛 are positive (or all

negative);
 Hyperbolic if none of 𝜆1 , . . . , 𝜆𝑛 vanish and one of them has the
opposite sign from the (𝑛 − 1) others;
 Parabolic if exactly one of 𝜆1 , . . . , 𝜆𝑛 is zero and all others have
the same sign.
 13

2 Separation of Variables Technique for Solving Linear PDEs

It is an important technique for solving linear PDEs provided that they

are separable. But, first some basic definitions towards solving PDEs.

Def : Superposition principle

 If 𝑢1 solves the linear PDE ℒ [𝑢1 ] = 𝑓1 and 𝑢2 solves ℒ [𝑢2 ] = 𝑓2 , then

𝑢 = 𝑐1 𝑢1 + 𝑐2 𝑢2 solves
ℒ[𝑢] = 𝑐1 𝑓1 + 𝑐2 𝑓2
for some functions 𝑓1 , 𝑓2 and constants 𝑐1 , 𝑐2 . In particular, if 𝑢1 and
𝑢2 both solve the same homogeneous (𝑓1 = 𝑓2 = 0) linear PDE, so
does 𝑢 = 𝑐1 𝑢1 + 𝑐2 𝑢2 .
 If 𝑢1 satisfies the linear boundary condition 𝐵[𝑢]|𝜕Ω = 𝑔1 and 𝑢2
satisfies 𝐵[𝑢]|𝜕Ω = 𝑔2 on the boundary 𝜕Ω, then 𝑢 = 𝑐1 𝑢1 + 𝑐2 𝑢2
satisfies
𝐵 [𝑢]|𝜕Ω = 𝑐1 𝑔1 + 𝑐2 𝑔2 .
In particular, if 𝑢1 and 𝑢2 both satisfy the same homogeneous (𝑔1 =
𝑔2 = 0) linear boundary condition, so does 𝑢 = 𝑐1 𝑢1 + 𝑐2 𝑢2 .

e.g. Let 𝑢1 (𝑥, 𝑡) and 𝑢2 (𝑥, 𝑡) satisfy 1D homogeneous wave equation

𝑢𝑡𝑡 = 𝑢𝑥𝑥

subject to the homogeneous boundary conditions

𝑢(0, 𝑡) = 𝑢(𝐿, 𝑡) = 0,
 14

then 𝑢 = 𝑐1 𝑢1 + 𝑐2 𝑢2 satisfies

𝑢𝑡𝑡 = (𝑐1 𝑢1 + 𝑐2 𝑢2 )𝑡𝑡 = 𝑐1 (𝑢1 )𝑡𝑡 + 𝑐2 (𝑢2 )𝑡𝑡 = 𝑐1 (𝑢1 )𝑥𝑥 + 𝑐2 (𝑢2 )𝑥𝑥
= (𝑐1 𝑢1 + 𝑐2 𝑢2 )𝑥𝑥 = 𝑢𝑥𝑥

and

𝑢(0, 𝑡) = 𝑐1 𝑢1 (0, 𝑡) + 𝑐2 𝑢2 (0, 𝑡) = 0,

𝑢(𝐿, 𝑡) = 𝑐1 𝑢1 (𝐿, 𝑡) + 𝑐2 𝑢2 (𝐿, 𝑡) = 0.

Thus, the superposition of 𝑢1 (𝑥, 𝑡) and 𝑢2 (𝑥, 𝑡);

𝑢(𝑥, 𝑡) = 𝑐1 𝑢1 + 𝑐2 𝑢2

satisfies the same linear, homogeneous PDE with homogeneous

boundary conditions as 𝑢1 and 𝑢2 .

e.g. Consider the (elliptic) Poisson problem

ℒ[𝑢] ≡ 𝑢𝑥𝑥 + 𝑢𝑦𝑦 = 𝐹(𝑥, 𝑦)

subject to the boundary conditions

𝑢(0, 𝑦) = 𝑔1 (𝑦), 𝑢(𝑎, 𝑦) = 𝑔2 (𝑦),

𝑢(𝑥, 0) = ℎ1 (𝑥), 𝑢(𝑥, 𝑏) = ℎ2 (𝑥).

 15

then the solution can be written as

𝑢 = 𝑢1 + 𝑢2 + 𝑢3

such that 𝑢1 (𝑥, 𝑦) satisfies

ℒ [𝑢1 ] ≡ (𝑢1 )𝑥𝑥 + (𝑢1 )𝑦𝑦 = 𝐹 (𝑥, 𝑦),

subject to the homogeneous boundary conditions

𝑢1 (0, 𝑦) = 𝑢1 (𝑎, 𝑦) = 𝑢1 (𝑥, 0) = 𝑢1 (𝑥, 𝑏) = 0,

and 𝑢2 (𝑥, 𝑦) satisfies the (homogeneous) Laplace equation

ℒ [𝑢2 ] ≡ (𝑢2 )𝑥𝑥 + (𝑢2 )𝑦𝑦 = 0,

subject to the (semi-) nonhomogeneous boundary conditions

𝑢2 (0, 𝑦) = 0, 𝑢2 (𝑎, 𝑦) = 0,

𝑢2 (𝑥, 0) = ℎ1 (𝑥), 𝑢2 (𝑥, 𝑏) = ℎ2 (𝑥).

and 𝑢3 (𝑥, 𝑦) satisfies the (homogeneous) Laplace equation

ℒ [𝑢3 ] ≡ (𝑢3 )𝑥𝑥 + (𝑢3 )𝑦𝑦 = 0,

subject to the (semi-) nonhomogeneous boundary conditions

𝑢3 (0, 𝑦) = 𝑔1 (𝑦), 𝑢3 (𝑎, 𝑦) = 𝑔2 (𝑦),

𝑢3 (𝑥, 0) = 0, 𝑢3 (𝑥, 𝑏) = 0.
 16

The Implementation of Separation of Variables Technique:

Recall that a conservation law is shown to lead to the spatial differential

operator involving a quantity 𝑢:

ℒ [𝑢] ≡ −𝛻 ⋅ (𝑝(𝒙)𝛻𝑢) + 𝑞(𝒙)𝑢

that embodies the internal effects proportional to 𝑢 acting in the spatial
domain 𝛺 and the flux through its boundary 𝜕𝛺.

In particular, over 1D domains 𝛺 = { 𝑥 | 𝑎 < 𝑥 < 𝑏 }, it becomes

𝜕 𝜕𝑢
ℒ [𝑢] ≡ − (𝑝(𝑥 ) ) + 𝑞(𝑥 )𝑢.
𝜕𝑥 𝜕𝑥
The balance between the time variation and the spatial variation of the
quantity 𝑢 may then lead to a parabolic PDE

𝜕𝑢
+ ℒ [𝑢] = 𝐹 (𝑥, 𝑡),
𝜕𝑡
or an hyperbolic PDE

𝜕2𝑢
2
+ ℒ[𝑢] = 𝐹 (𝑥, 𝑡),
𝜕𝑡
with the nonhomogeneous term 𝐹 (𝑥, 𝑡).
 17

Method of separation of variables is one of the most widely used

techniques to solve PDEs and is based on the assumption that the solution
of the equation is separable, that is, the final solution can be represented
as a product of several functions, each of which is only dependent upon
a single independent variable such as

𝑢(𝑥, 𝑡) = 𝑋(𝑥)𝑇(𝑡).

 When substituted into, say, a 1D homogeneous parabolic PDE

𝑇′ 𝑋 + 𝑇 ℒ [𝑋] = 0

and divided by 𝑇𝑋, it yields

𝑇′(𝑡) ℒ [𝑋(𝑥 )]
=− = −𝜆
𝑇(𝑡) 𝑋(𝑥 )

where 𝜆 is called separation constant. Since 𝑥 and 𝑡 are independent

variables, if 𝑥 is fixed and 𝑡 varies, then the left-hand side is a constant.
Analogously, if 𝑡 is fixed and 𝑥 varies, then the right-hand side is a
constant. Thus, the separation constant 𝜆 represents the common constant
value of both sides. As a result, we obtain two ODEs

ℒ[𝑋(𝑥 )] = 𝜆 𝑋(𝑥) and 𝑇′(𝑡) + 𝜆 𝑇(𝑡) = 0.

 18

 In 1D homogeneous hyperbolic case, one obtains

𝑇′′(𝑡) ℒ [𝑋(𝑥 )]
=− = −𝜆
𝑇(𝑡) 𝑋(𝑥)

and thus

ℒ [𝑋(𝑥 )] = 𝜆 𝑋(𝑥) and 𝑇′′(𝑡) + 𝜆 𝑇(𝑡) = 0.

 Let the parabolic and hyperbolic cases be considered for the

homogeneous boundary conditions

𝛼𝑢(𝑎, 𝑡) + 𝛽𝑢𝑥 (𝑎, 𝑡) = 0 & 𝛾𝑢(𝑏, 𝑡) + 𝛿𝑢𝑥 (𝑏, 𝑡) = 0.

that can be separated by substituting

𝑢(𝑥, 𝑡) = 𝑋(𝑥)𝑇(𝑡)

to get

{𝛼𝑋(𝑎) + 𝛽𝑋′(𝑎 )}𝑇(𝑡) = 0 & [𝛾𝑋(𝑏) + 𝛿𝑋′(𝑏)]𝑇(𝑡) = 0.

For, 𝑇(𝑡) ≠ 0 (otherwise trivial solution results), they reduce to

𝛼𝑋(𝑎) + 𝛽𝑋 ′ (𝑎) = 0 & 𝛾𝑋(𝑏) + 𝛿𝑋′(𝑏) = 0

respectively.

Note that in each case (parabolic & hyperbolic), the equation

ℒ [𝑋(𝑥 )] = 𝜆 𝑋(𝑥)

with the homogeneous boundary conditions

 19

𝛼𝑋(𝑎) + 𝛽𝑋 ′ (𝑎) = 0 & 𝛾𝑋(𝑏) + 𝛿𝑋′(𝑏) = 0

is a differential eigenproblem whose nonzero (nontrivial) soln’s 𝑋(𝑥) are

called the eigenfunction corresponding to the eigenvalue 𝜆. This problem
for ℒ [∙] self-adjoint is already studied as the Sturm-Liouville (SL)
problem.

SL-problem is guaranteed to have has countably ∞-many nonnegative

real eigenvalues 𝜆𝑘 ≥ 0, 𝑘 = 1,2. .. such that 𝜆𝑘 → ∞ as 𝑘 → ∞ with
corresponding real eigenfunctions 𝑋𝑘 (𝑥). Further, the eigenfunctions are
orthogonal based on the the inner product, i.e.
𝐿
⟨𝑋𝑘 , 𝑋𝑗 ⟩ ≡ ∫ 𝑋𝑘 (𝑥 ) 𝑋𝑗 (𝑥 ) 𝑑𝑥 = 0
0

for 𝑘 ≠ 𝑗.

Assuming 𝜆𝑘 = 𝛼𝑘2 > 0 for simplicity, for each eigenvalue 𝜆𝑘 we obtain

an equation for 𝑇𝑘 (𝑡) or 𝑌𝑘 (𝑦) that yield:

 hyperbolic case,

𝑇′′(𝑡) + 𝛼𝑘2 𝑇(𝑡) = 0 ⇒ 𝑇𝑘 (𝑡) = 𝐴𝑘 𝑐𝑜𝑠( 𝛼𝑘 𝑡) + 𝐵𝑘 𝑠𝑖𝑛( 𝛼𝑘 𝑡).

 parabolic case,

𝑇′(𝑡) + 𝛼𝑘2 𝑇(𝑡) = 0 ⇒ 𝑇𝑘 (𝑡) = 𝐴𝑘 𝑒𝑥𝑝( − 𝛼𝑘2 𝑡).

By construction, the product (kth mode)

 20

𝑢𝑘 (𝑥, 𝑡) = 𝑇𝑘 (𝑡)𝑋𝑘 (𝑥)

satisfies the corresponding PDEs for each 𝑘 together with the

homogeneous boundary conditions. However, for complete
representation of the solution, they are superposed to get
∞
𝑢(𝑥, 𝑡) = ∑ 𝑢𝑘 (𝑥, 𝑡).
𝑘=1

In order to determine the coefficients 𝐴𝑘 , 𝐵𝑘 , the additional conditions

must be supplied, such as initial conditions.

 In the hyperbolic case, the two initial conditions must be satisfied

∞ ∞
𝑢(𝑥, 0) = ∑ 𝑇𝑘 (0)𝑋𝑘 (𝑥 ) = ∑ 𝐴𝑘 𝑋𝑘 (𝑥 ) = 𝑓 (𝑥 ),
𝑘=1 𝑘=1
∞ ∞
𝑢𝑡 (𝑥, 0) = ∑ 𝑇𝑘′ (0) 𝑋𝑘 (𝑥 ) = ∑ 𝛼𝑘 𝐵𝑘 𝑋𝑘 (𝑥 ) = 𝑔(𝑥 ).
𝑘=1 𝑘=1

The orthogonality of the SL-eigenfunctions yields

⟨𝑓(𝑥), 𝑋𝑘 (𝑥)⟩ ⟨𝑔(𝑥), 𝑋𝑘 (𝑥)⟩
𝐴𝑘 = & 𝛼𝑘 𝐵𝑘 =
⟨𝑋𝑘 (𝑥), 𝑋𝑘 (𝑥)⟩ ⟨𝑋𝑘 (𝑥), 𝑋𝑘 (𝑥)⟩

and thus the solution is completely determined

∞
𝑢(𝑥, 𝑡) = ∑ [𝐴𝑘 𝑐𝑜𝑠(𝛼𝑘 𝑡) + 𝐵𝑘 𝑠𝑖𝑛(𝛼𝑘 𝑡)] 𝑋𝑘 (𝑥 ).
𝑘=1
 21

 In the parabolic case, the initial condition to be satisfied is

∞ ∞
𝑢(𝑥, 0) = ∑ 𝑇𝑘 (0)𝑋𝑘 (𝑥 ) = ∑ 𝐴𝑘 𝑋𝑘 (𝑥 ) = 𝑓(𝑥)
𝑘=1 𝑘=1

This yields
⟨𝑓(𝑥), 𝑋𝑘 (𝑥)⟩
𝐴𝑘 =
⟨𝑋𝑘 (𝑥), 𝑋𝑘 (𝑥)⟩
and thus the complete solution is
∞
𝑢(𝑥, 𝑡) = ∑ 𝐴𝑘 𝑒𝑥𝑝( − 𝛼𝑘2 𝑡) 𝑋𝑘 (𝑥).
𝑘=1

Some Specific Examples:

e.g. Consider the homogeneous Wave equation (hyperbolic type)

𝑢𝑡𝑡 = 𝑐 2 𝑢𝑥𝑥 , 0 < 𝑥 < 𝐿, 𝑡>0

subject to the homogeneous boundary

𝑢(0, 𝑡) = 𝑢(𝐿, 𝑡) = 0, 𝑡>0

and the initial conditions

𝑢(𝑥, 0) = 𝑓(𝑥 ), 𝑢𝑡 (𝑥, 0) = 0, 0 < 𝑥 < 𝐿.

Assume separable solution 𝑢(𝑥, 𝑡) = 𝑋(𝑥 )𝑇(𝑡) and substitute into the
PDE to get

𝑋 𝑇 ″ = 𝑐 2 𝑋 ″ 𝑇,
 22

and in turn,

𝑇″ 𝑋″
= =𝜆
𝑐2 𝑇 𝑋
where 𝛼 is the separation constant. It yields two ODEs:

(1) 𝑋 ″ − 𝜆𝑋 = 0,

(2) 𝑇 ″ − 𝜆𝑐 2 𝑇 = 0.

The homogeneous boundary conditions then yield

𝑢(0, 𝑡) = 𝑋(0)𝑇(𝑡) = 0 ⇒ 𝑋(0) = 0

𝑢(𝐿, 𝑡) = 𝑋(𝐿)𝑇(𝑡) = 0 ⇒ 𝑋(𝐿) = 0.

While the nonhomogeneous initial condition

𝑢(𝑥, 0) = 𝑓(𝑥)

will be handled later, the homogeneous initial condition yields

𝑢𝑡 (𝑥, 0) = 𝑋(𝑥)𝑇 ′ (0) = 0 ⇒ 𝑇 ′ (0) = 0.

Now, (1) forms a regular S-L problem:

𝑋 ″ − 𝜆𝑋 = 0, 𝑋(0) = 𝑋(𝐿) = 0

and for negative separation constant 𝜆 = −𝛼 2 it yields the eigensoln’s:

𝛼𝑘 = 𝑘𝜋/𝐿 & 𝑋𝑘 (𝑥) = 𝑠𝑖𝑛( 𝛼𝑘 𝑥), 𝑘 = 1,2, ….

Now, (2) follows with the solution

𝑇𝑘 (𝑡) = 𝐴𝑘 𝑐𝑜𝑠( 𝛼𝑘 𝑐𝑡) + 𝐵𝑘 𝑠𝑖𝑛( 𝛼𝑘 𝑐𝑡)

 23

and using 𝑇 ′ (0) = 0, it reduces to

𝑇𝑘 (𝑡) = 𝐴𝑘 𝑐𝑜𝑠( 𝛼𝑘 𝑐𝑡).

The form of the soln’s are then

𝑘𝜋 𝑘𝜋
𝑢𝑘 (𝑥, 𝑡) = 𝑋𝑘 (𝑥)𝑇𝑘 (𝑡) = 𝐴𝑘 𝑐𝑜𝑠( 𝑐𝑡) 𝑠𝑖𝑛( 𝑥),
𝐿 𝐿

that yields the solution

∞ ∞
𝑘𝜋 𝑘𝜋
𝑢(𝑥, 𝑡) = ∑ 𝑢𝑘 (𝑥, 𝑡) = ∑ 𝐴𝑘 𝑐𝑜𝑠( 𝐿
𝑐𝑡) 𝑠𝑖𝑛( 𝐿
𝑥)
𝑘=1 𝑘=1

by superposition principle.

Note that the nonhomogeneous initial condition 𝑢(𝑥, 0) = 𝑓(𝑥) is not

satisfied yet and it will be used to determine the expansion coefficients
𝐴𝑘 by the sine expansion:
∞
𝑘𝜋
𝑓(𝑥 ) = ∑ 𝐴𝑘 𝑠𝑖𝑛( 𝐿
𝑥).
𝑘=1

Consider the four cases for 𝑓(𝑥):

(i)
𝜋𝑥 ∞
𝑘𝜋
𝑓(𝑥) = 3 𝑠𝑖𝑛( ) = ∑ 𝐴𝑘 𝑠𝑖𝑛( 𝐿 𝑥).
𝐿 𝑘=1

Note that it is satisfied for 𝐴1 = 3 and 𝐴𝑘>1 = 0 that yields the solution
𝜋 𝜋
𝑢(𝑥, 𝑡) = 3
⏟𝑐𝑜𝑠( 𝐿 𝑐𝑡) 𝑠𝑖𝑛( 𝐿 𝑥),
Amp(𝑡)
 24

simply by comparison. Verify that the orthogonality of the S-L

eigenfunction yields the same solution.

This is the fundamental mode which is a standing wave with amplitude

𝜋
Amp(𝑡) = 3 𝑐𝑜𝑠( 𝐿 𝑐𝑡)

at any given x-value and frequency

𝜈 = 𝑐/2𝐿.

Using trigonometric identities, it can also be written as

3 𝜋 𝜋
𝑢(𝑥, 𝑡) = 2 [𝑠𝑖𝑛 ( 𝐿 (𝑥 + 𝑐𝑡)) + 𝑠𝑖𝑛 ( 𝐿 (𝑥 − 𝑐𝑡))]
1
= 2[𝑓(𝑥 + 𝑐𝑡) + 𝑓(𝑥 − 𝑐𝑡)].

That is two travelling waves,

3 𝜋
𝑠𝑖𝑛 ( 𝐿 (𝑥 + 𝑐𝑡)) (left-going)
2

and
3 𝜋
𝑠𝑖𝑛 ( 𝐿 (𝑥 − 𝑐𝑡)) (right-going)
2

superpose to yield the standing wave solution

𝜋 𝜋
𝑢(𝑥, 𝑡) = (3
⏟ 𝑐𝑜𝑠( 𝐿 𝑐𝑡)) 𝑠𝑖𝑛( 𝐿 𝑥)
Amp(𝑡)
 25

as shown below:

Note that the profile 𝑓(𝑥) is needed to be extended beyond the actual
interval 0 < 𝑥 < 𝐿 to visualize the travelling wave pattern. This is
already done as a result of the globally defined sine-expansion.

(ii)
𝜋𝑥 2𝜋𝑥 ∞
𝑘𝜋
𝑓(𝑥) = 3 𝑠𝑖𝑛( ) − 𝑠𝑖𝑛( )=∑ 𝐴𝑘 𝑠𝑖𝑛( 𝐿 𝑥).
𝐿 𝐿 𝑘=1

It is satisfied when 𝐴1 = 3, 𝐴2 = −1 and 𝐴𝑘>2 = 0 that yields

𝜋 𝜋 2𝜋 2𝜋
𝑢(𝑥, 𝑡) = 3
⏟𝑐𝑜𝑠( 𝐿 𝑐𝑡) 𝑠𝑖𝑛( 𝐿 𝑥) − 𝑐𝑜𝑠(
⏟ 𝑐𝑡) 𝑠𝑖𝑛( 𝑥).
𝐿 𝐿
Fundamental Mode The Second Harmonic

again simply by comparison. Verify also that the orthogonality of the S-

L eigenfunction yields the same solution.

The solution is the sum of two modes. The first one is the fundamental
mode while the second one is the second harmonic with amplitude
2𝜋
Amp(𝑡) = 𝑐𝑜𝑠( 𝑐𝑡),
𝐿
 26

at any given x-value and frequency

𝜈 = 𝑐/𝐿,

i.e. the second harmonic oscillates twice as fast.

(iii)
∞
𝑘𝜋
𝑓(𝑥) = 𝑥(𝐿 − 𝑥) = ∑ 𝐴𝑘 𝑠𝑖𝑛( 𝐿
𝑥).
𝑘=1

It clearly requires the use of the orthogonality of the SL-eigenfunctions

𝑋𝑘 (𝑥) that provides the coefficients by

2 𝐿 𝑘𝜋 4𝐿2 1
𝐴𝑘 = ∫ 𝑥(𝐿 − 𝑥) 𝑠𝑖𝑛( 𝑥)𝑑𝑥 = 3 3 [1 + (−1)𝑘+1 ] ∝
𝐿 0 𝐿 𝑘 𝜋 𝑘3

and, in turn, the solution

∞
𝑘𝜋 𝑘𝜋
𝑢(𝑥, 𝑡) = ∑ 𝐴𝑘 𝑐𝑜𝑠( 𝐿
𝑐𝑡) 𝑠𝑖𝑛( 𝐿
𝑥).
𝑘=1

It is formed by the contribution of infinitely many modes (harmonics),

each having an amplitude
 27

𝑘𝜋 1
Amp𝑘 (𝑡) = 𝐴𝑘 𝑐𝑜𝑠( 𝐿 𝑐𝑡) ∝
𝑘3
thus lower modes are contributing more significantly than higher modes
and each with frequency 𝜈𝑘 = 𝑘𝑐/2𝐿.

(iv)
∞
𝑥, 0 ≤ 𝑥 < 𝐿/2 𝑘𝜋
𝑓(𝑥) = { }=∑ 𝐴𝑘 𝑠𝑖𝑛( 𝐿 𝑥).
𝐿 − 𝑥, 𝐿/2 ≤ 𝑥 ≤ 𝐿 𝑘=1

In this case, the initial profile has a discontinuity in its derivative. By

using the orthogonality, the eigenfunction expansion coefficients 𝐴𝑘 can
be obtained as
𝐿/2 𝐿
2 𝑘𝜋 𝑘𝜋
𝐴𝑘 = { ∫ 𝑥 𝑠𝑖𝑛( 𝑥) 𝑑𝑥 + ∫ (𝐿 − 𝑥) 𝑠𝑖𝑛( 𝑥) 𝑑𝑥 }
𝐿 0 𝐿 𝐿/2 𝐿
(1 − cos(𝑛𝜋/2)) 1
= 4𝐿 sin(𝑛𝜋/2) ∝
𝑛2 𝜋 2 𝑘2
leading to the solution
∞ (1 − cos(𝑛𝜋/2)) 𝑘𝜋 𝑘𝜋
𝑢(𝑥, 𝑡) = ∑ 4𝐿 sin(𝑛𝜋/2) 𝑐𝑜𝑠( 𝑐𝑡) 𝑠𝑖𝑛( 𝑥).
𝑘=1 𝑛2 𝜋 2 𝐿 𝐿

Again, infinitely many modes (harmonics), each having an amplitude

proportional to
1
Amp𝑘 (𝑡) ∝
𝑘2
 28

in terms of the wave number 𝑘. As a consequence of the discontinuity in

the initial profile, the decay rate of the amplitude the wave number 𝑘 is
slower when compared with the amplitude in case (iii).

The superposition of the two travelling waves resulting in the standing

wave pattern in the interval 0 < 𝑥 < 𝐿 is shown in the sampled
sequence of plots covering half the oscillation period, below:

The Matlab codes that perform the simulations above, wavesup_i.m,

wavesup_ii.m and wavesup_iv.m, will be made available via
METUClass.
 29

e.g. Consider the Heat equation (parabolic type):

𝑢𝑡 = 𝜅𝑢𝑥𝑥 , 0 < 𝑥 < 𝐿, 𝑡>0

subject to the boundary

𝑢𝑥 (0, 𝑡) = 𝑢𝑥 (𝐿, 𝑡) = 0, 𝑡 > 0,

and initial conditions

𝑢(𝑥, 0) = 𝑓(𝑥 ), 0 < 𝑥 < 𝐿.

In order to implement the Separation of Variables technique, let

𝑢(𝑥, 𝑡) = 𝑋(𝑥)𝑇(𝑡)

and substitute into the PDE to get

𝑋 𝑇′ = 𝜅 𝑋′′ 𝑇.

After dividing by 𝜅𝑋𝑇, it becomes

𝑇′ 𝑋″
= =𝜆
𝜅𝑇 𝑋
where 𝛼 is the separation [Link] yields two ODEs:

(1) 𝑋 ″ − 𝜆𝑋 = 0,

(2) 𝑇′ − 𝜆𝜅𝑇 = 0.

The separation of the homogeneous boundary conditions lead to

𝑢𝑥 (0, 𝑡) = 𝑋′(0)𝑇(𝑡) = 0 ⇒ 𝑋′(0) = 0,

𝑢𝑥 (𝐿, 𝑡) = 𝑋′(𝐿)𝑇(𝑡) = 0 ⇒ 𝑋′(𝐿) = 0.

 30

The nonhomogeneous initial condition 𝑢(𝑥, 0) = 𝑓(𝑥) will be handled

later because it is not separable.

Now, (1) forms a regular S-L problem:

𝑋 ″ − 𝜆𝑋 = 0, 𝑋 ′ (0) = 𝑋 ′ (𝐿) = 0

and for 𝜆 = −𝛼 2 it yields the eigensolutions:

𝛼𝑘 = 𝑘𝜋⁄𝐿 & 𝑋𝑘 (𝑥) = 𝑐𝑜𝑠( 𝛼𝑘 𝑥), 𝑘 = 0,1, . ..

Note that 𝜆0 = 0 leads to a nontrivial eigenfunction 𝑋0 (𝑥 ) = 1.

Thus, (2) follows

(2) 𝑇′ + 𝛼𝑘2 𝜅𝑇 = 0 ⟹ 𝑇𝑘 (𝑥) = 𝐴𝑘 𝑒𝑥𝑝( − 𝛼𝑘2 𝜅𝑡).

The form of the solution for one mode is then

𝑘𝜋 𝑘𝜋
𝑢𝑘 (𝑥, 𝑡) = 𝑋𝑘 (𝑥)𝑇𝑘 (𝑡) = 𝐴𝑘 𝑒𝑥𝑝( − ( 𝐿 )2 𝜅𝑡) 𝑐𝑜𝑠( 𝑥),
𝐿

that provides the complete solution

∞ 𝑘𝜋 𝑘𝜋
𝑢(𝑥, 𝑡) = ∑ 𝐴𝑘 𝑒𝑥𝑝( − ( )2 𝜅𝑡) 𝑐𝑜𝑠( 𝑥).
𝑘=0 𝐿 𝐿

Let us now proceed to satisfy the initial condition

∞
𝑘𝜋
𝑢(𝑥, 0) = 𝑓(𝑥 ) = ∑ 𝐴𝑘 𝑐𝑜𝑠( 𝑥)
𝐿
𝑘=0

that is in the form of an eigenfunction expansion. Using the

orthogonality of the S-L eigenfunctions, one gets
 31

1 𝐿
𝐴0 = ∫ 𝑓(𝑥 )𝑑𝑥 ,
𝐿 0

2 𝐿 𝑘𝜋
𝐴𝑘 = ∫ 𝑓(𝑥) 𝑐𝑜𝑠( 𝑥)𝑑𝑥 for 𝑘 = 1,2, …
𝐿 0 𝐿

For the special choice of 𝑓(𝑥) = 𝑥, the coefficients become

𝐿 2𝐿
𝐴0 = , 𝐴𝑘 = 2 2
[(−1)𝑘 − 1].
2 𝑘 𝜋
Then the solution is
∞
𝐿 2𝐿 𝑘
𝑘𝜋 2 𝑘𝜋
( ) [ ]
𝑢 𝑥, 𝑡 = + ∑ 2 2 (−1) − 1 𝑒𝑥𝑝( − ( ) 𝜅𝑡) 𝑐𝑜𝑠( 𝑥).
2 𝑘 𝜋 𝐿 𝐿
𝑘=1

Note that as 𝑡 → ∞, 𝑢 → 𝐿⁄2 exponentally fast. This is the constant

equilibrium temperature. Since the ends of the domain are insulated, the
temperature is equalized in the domain to an average value (dotted-line).

Note also that the intermediate solutions satisfy vanishing derivative at

the ends even though the inital profile (dash-line) does not.
 32

e.g. Consider the Laplace equation (elliptic type):

𝑢𝑥𝑥 + 𝑢𝑦𝑦 = 0, 0 < 𝑥 < 𝑎, 0<𝑦<𝑏

subject to the boundary conditions

𝑢(0, 𝑦) = 𝑢(𝑎, 𝑦) = 0, 0 < 𝑦 < 𝑏,

𝑢(𝑥, 0) = 1, 𝑢(𝑥, 𝑏) = 0, 0 < 𝑥 < 𝑎.

Assume the solution to be separabel 𝑢(𝑥, 𝑦) = 𝑋(𝑥)𝑌(𝑦) and substitute

into the PDE to get

𝑋′′ 𝑌 + 𝑋 𝑌′′ = 0

that becomes

𝑋″ 𝑌″
=− =𝜆
𝑋 𝑌
after dividing by 𝑋𝑌. This yields two ODEs:

(1) 𝑋′′ − 𝜆𝑋 = 0,

(2) 𝑌′′ + 𝜆𝑌 = 0.

Along 𝑥, the boundary conditions are homogeneous, so separable

𝑢(0, 𝑦) = 𝑋(0)𝑌(𝑦) = 0 ⇒ 𝑋(0) = 0,

𝑢(𝑎, 𝑦) = 𝑋(𝑎)𝑌(𝑦) = 0 ⇒ 𝑋(𝑎) = 0,

while along 𝑦, the one homogeneous boundary condition is separable

𝑢(𝑥, 𝑏) = 𝑋(𝑥)𝑌(𝑏) = 0 ⇒ 𝑌(𝑏) = 0.

 33

The other boundary condition 𝑢(𝑥, 0) = 1 is nonhomogeneous and so it

is not separable. That is why it will be handled later.

Clearly, (1) forms a regular S-L problem:

𝑋 ″ − 𝛼𝑋 = 0, 𝑋(0) = 𝑋(𝑎) = 0

and for 𝜆 = −𝛼 2 , it yields the eigensolutions:

𝛼𝑘 = 𝑘𝜋/𝑎 & 𝑋𝑘 (𝑥) = 𝑠𝑖𝑛( 𝛼𝑘 𝑥) for 𝑘 = 1,2, …

It then follows that

(2) 𝑌 ′′ − 𝛼𝑘2 𝑌 = 0 ⟹ 𝑌𝑘 (𝑦) = 𝐴𝑘 𝑐𝑜𝑠ℎ( 𝛼𝑘 𝑦) + 𝐵𝑘 𝑠𝑖𝑛ℎ( 𝛼𝑘 𝑦).

The separated homogeneous condition 𝑌(𝑏) = 0 gives one relation

𝑐𝑜𝑠ℎ( 𝑘𝜋𝑏/𝑎)
0 = 𝐴𝑘 𝑐𝑜𝑠ℎ( 𝛼𝑘 𝑏) + 𝐵𝑘 𝑠𝑖𝑛ℎ( 𝛼𝑘 𝑏) ⟹ 𝐵𝑘 = −𝐴𝑘 .
𝑠𝑖𝑛ℎ( 𝑘𝜋𝑏/𝑎)

The satisfaction of the nonhomogeneous condition 𝑢(𝑥, 0) = 1 requires

superposition of ∞–many solutions, i.e.
∞ 𝑘𝜋 𝑘𝜋 𝑘𝜋
𝑢(𝑥, 𝑡) = ∑ [𝐴𝑘 𝑐𝑜𝑠ℎ( 𝑦) + 𝐵𝑘 𝑠𝑖𝑛ℎ( 𝑦)] 𝑠𝑖𝑛( 𝑥).
𝑘=1 𝑎 𝑎 𝑎

so that
∞
𝑘𝜋
𝑢(𝑥, 0) = 1 = ∑ 𝐴𝑘 𝑠𝑖𝑛( 𝑥).
𝑎
𝑘=1

Using the orthogonality of the SL-eigenfunctions, one gets

 34

2 𝑎 𝑘𝜋 2
𝐴𝑘 = ∫ 𝑠𝑖𝑛( 𝑥)𝑑𝑥 = [1 + (−1)𝑘+1 ]
𝑎 0 𝑎 𝑘𝜋

and it follows that

2 𝑐𝑜𝑠ℎ( 𝑘𝜋𝑏/𝑎)
𝐵𝑘 = − [1 + (−1)𝑘+1 ] .
𝑘𝜋 𝑠𝑖𝑛ℎ( 𝑘𝜋𝑏/𝑎)

The solution then becomes

2[1+(−1)𝑘+1 ] 𝑘𝜋 𝑐𝑜𝑠ℎ(𝑘𝜋𝑏/𝑎) 𝑘𝜋 𝑘𝜋
𝑢(𝑥, 𝑦) = ∑∞
𝑘=1 (𝑐𝑜𝑠ℎ( 𝑦) − 𝑠𝑖𝑛ℎ( 𝑦)) 𝑠𝑖𝑛( 𝑥).
𝑘𝜋 𝑎 𝑠𝑖𝑛ℎ(𝑘𝜋𝑏/𝑎) 𝑎 𝑎

Note that at the corners (0 , 0) and (𝑎 , 0), convergence in the mean

occurs.
 35

Nonhomogeneous Problems:

When nonhomogeneities are present in the PDE or in the boundary

conditions, the direct application of the Separation of Variables technique
is not possible. However, by specialized use of the superposition principle
and the eigenfunction expansion, this limitation can be overcome. In the
following examples, we will consider two types of nonhomogeneities:
time-dependent and time-independent.

Time-independent nonhomogeneities

e.g. Consider the vibration problem of a taut string under gravity:

𝑢𝑡𝑡 = 𝑐 2 𝑢𝑥𝑥 − 𝑔, 0 < 𝑥 < 𝐿, 𝑡 > 0,

subject to boundary

𝑢(0, 𝑡) = 𝑢(𝐿, 𝑡) = 0, 𝑡 > 0,

and initial conditions

𝑢(𝑥, 0) = 𝑓(𝑥 ), 𝑢𝑡 (𝑥, 0) = 0, 0 < 𝑥 < 𝐿.

Here 𝑔 is the gravitational constant and it forms the time-independent

nonhomogeneity.

When the assumed separable solution 𝑢(𝑥, 𝑡) = 𝑋(𝑥)𝑇(𝑡) is substituted

into the nonhomogeneous PDE, it yields
 36

𝑇 ′′ 𝑋 ′′ 𝑔
𝑋𝑇 ′′ = 𝑐 2 𝑋 ′′ 𝑇 − 𝑔 ⟹ = − .
𝑐2𝑇 𝑋 𝑐 2 𝑋𝑇
Clearly the resulting form can not be separated.

In order to induce the separability of the PDE, we introduce a new

dependent variable 𝜓(𝑥) such that

𝑢(𝑥, 𝑡) = 𝑣(𝑥, 𝑡) + 𝜓(𝑥 ).

When this is substituted into the PDE, it becomes

𝑣𝑡𝑡 = 𝑐 2 𝑣𝑥𝑥 + [𝑐 2 𝜓 ′′ − 𝑔]

together with the boundary and initial conditions:

𝑣(0, 𝑡) + 𝜓(0) = 𝑣(𝐿, 𝑡) + 𝜓(𝐿) = 0,

𝑣(𝑥, 0) + 𝜓(𝑥 ) = 𝑓(𝑥 ),

𝑣𝑡 (𝑥, 0) = 0.

In order to remove the nonhomogeneity in the resulting PDE, we set

𝑐 2 𝜓 ′′ − 𝑔 = 0,

subject to the boundary conditions

𝜓(0) = 𝜓(𝐿) = 0,

whose solution can easily be shown to be

𝑔
𝜓(𝑥 ) = 𝑥 (𝑥 − 𝐿).
2𝑐 2
 37

This solution corresponds to static deflection of the string motionless

under the gravity.

Now, the modified PDE for 𝑣 = 𝑣(𝑥, 𝑡) takes the form

𝑣𝑡𝑡 = 𝑐 2 𝑣𝑥𝑥 , 0 < 𝑥 < 𝐿, 𝑡 > 0,

𝑣(0, 𝑡) = 𝑣 (𝐿, 𝑡) = 0, 𝑡 > 0,

𝑣(𝑥, 0) = 𝑓 (𝑥 ) − 𝜓(𝑥 ), 𝑣𝑡 (𝑥, 0) = 0, 0 < 𝑥 < 𝐿.

It was already shown before that it admits the solution

∞
𝑘𝜋 𝑘𝜋
𝑣(𝑥, 𝑡) = ∑ 𝐴𝑘 𝑐𝑜𝑠( 𝑐𝑡) 𝑠𝑖𝑛( 𝑥)
𝐿 𝐿
𝑘=1

where

2 𝐿 𝑘𝜋
𝐴𝑘 = ∫ [𝑓(𝑥) − 𝜓(𝑥)] 𝑠𝑖𝑛( 𝑥)𝑑𝑥 .
𝐿 0 𝐿

Then, the solution to the original problem becomes

∞
𝑘𝜋 𝑘𝜋
𝑢(𝑥, 𝑡) = 𝜓(𝑥 ) + ∑ 𝐴𝑘 𝑐𝑜𝑠( 𝑐𝑡) 𝑠𝑖𝑛( 𝑥)
𝐿 𝐿
𝑘=1
∞
𝑔 𝑘𝜋 𝑘𝜋
= 2 𝑥 (𝑥 − 𝐿) + ∑ 𝐴𝑘 𝑐𝑜𝑠( 𝑐𝑡) 𝑠𝑖𝑛( 𝑥).
2𝑐 𝐿 𝐿
𝑘=1

This represents the oscillations around the static deflection profile 𝜓(𝑥).
 38

e.g. Consider the temperature distribution in a rod with constant nonzero

temperatures at the ends:

𝑢𝑡 = 𝜅𝑢𝑥𝑥 , 0 < 𝑥 < 𝐿, 𝑡>0

subject to boundary

𝑢(0, 𝑡) = 𝑈0 , 𝑢(𝐿, 𝑡) = 𝑈𝐿 , 𝑡 > 0,

and initial conditions

𝑢(𝑥, 0) = 𝑓 (𝑥 ), 0 < 𝑥 < 𝐿.

The time-independent nonhomogeneity in this case is in the boundary

conditions.

Before attempting the separation of variables, we again introduce a new

dependent variable

𝑢(𝑥, 𝑡) = 𝑣(𝑥, 𝑡) + 𝜓(𝑥)

and substitute into the PDE to get

𝑣𝑡 = 𝜅 𝑣𝑥𝑥 + [𝜅 𝜓′′ ]

and

𝑣(0, 𝑡) + 𝜓(0) = 𝑈0 , 𝑣(𝐿, 𝑡) + 𝜓(𝐿) = 𝑈𝐿 .

Thus, the selection

𝜓 ″ = 0, 𝜓(0) = 𝑈0 , 𝜓(𝐿) = 𝑈𝐿

with solution
𝑈𝐿 − 𝑈0
𝜓(𝑥 ) = 𝑈0 + 𝑥
𝐿
 39

help remove the nonhomogeneity resulting in

𝑣𝑡 = 𝜅 𝑣𝑥𝑥 , 0 < 𝑥 < 𝐿, 𝑡 > 0,

𝑣(0, 𝑡) = 𝑣 (𝐿, 𝑡) = 0, 𝑡 > 0,

𝑣(𝑥, 0) = 𝑓 (𝑥 ) − 𝜓(𝑥 ), 0 < 𝑥 < 𝐿.

This reduced separable problem can be solved to get

∞
𝑘𝜋 2 𝑘𝜋
( )
𝑣 𝑥, 𝑡 = ∑ 𝐴𝑘 𝑒𝑥𝑝( − ( ) 𝜅𝑡) 𝑠𝑖𝑛( 𝑥)
𝐿 𝐿
𝑘=1

where

2 𝐿 𝑘𝜋
𝐴𝑘 = ∫ [𝑓(𝑥) − 𝜓(𝑥)] 𝑠𝑖𝑛( 𝑥)𝑑𝑥 .
𝐿 0 𝐿

Subsequently, the solution to the original problem is obtained as

∞
𝑘𝜋 𝑘𝜋
𝑢(𝑥, 𝑡) = 𝜓(𝑥 ) + ∑ 𝐴𝑘 𝑒𝑥𝑝( − ( )2 𝜅𝑡) 𝑠𝑖𝑛( 𝑥)
𝐿 𝐿
𝑘=1
∞
𝑈𝐿 − 𝑈0 𝑘𝜋 𝑘𝜋
= 𝑈0 + 𝑥 + ∑ 𝐴𝑘 𝑒𝑥𝑝( − ( )2 𝜅𝑡) 𝑠𝑖𝑛( 𝑥).
𝐿 𝐿 𝐿
𝑘=1

Note that as 𝑡 → ∞, the solution converges to equilibrium sol’n 𝜓(𝑥 ).

Let as consider some choices for 𝑓(𝑥):

(i) Consider 0𝑜 𝐶 initial temperature: 𝑓(𝑥) = 0, then

2 𝐿 𝑘𝜋 2
𝐴𝑘 = ∫ −𝜓(𝑥 ) 𝑠𝑖𝑛( 𝑥)𝑑𝑥 = − [𝑈0 + (−1)𝑘+1 𝑈𝐿 ]
𝐿 0 𝐿 𝑘𝜋
 40

and the resulting temperature distribution is shown below:

The time evolution of 𝑢 = 𝑢(𝑥, 𝑡) starting from 0𝑜 𝐶 initial temperature

towards linear equilibrium profile looks smooth, however, when 𝑡 = 0
at 𝑥 = 0 and 𝑥 = 𝐿, the slope (i.e. heat flow) is infinite. This is still
handled by the eigenfunction expansion (weak sol’n) approach.

(ii) Consider a smooth initial profile:

𝑥 𝑥
𝑓(𝑥) = 𝑈0 (1 − ( )2 ) + 𝑈𝐿 ( )
𝐿 𝐿
such that 𝑓(0) = 𝑈0 and 𝑓(𝐿) = 𝑈𝐿 . Then

2 𝐿 𝑘𝜋 2 3
𝐴𝑘 = ∫ [𝑓(𝑥) − 𝜓(𝑥)] 𝑠𝑖𝑛( 𝑥)𝑑𝑥 = −( ) 𝑈0 [1 + (−1)𝑘+1 ],
𝐿 0 𝐿 𝑘𝜋

and the resulting temperature distribution looks like:

 41

The time evolution towards linear equilibrium profile, in this case, is

smooth for the initial profile 𝑓(𝑥) satisfying the boundary conditions.
Note the faster decay of the expansion coefficients 𝐴𝑘 ∝ 𝑘 −3 in
contrast to the previous case where 𝐴𝑘 ∝ 𝑘 −1 .

Time-dependent nonhomogeneities

Recall in ODEs that Variation of Parameters for a particular solution

involves writing

𝑦𝑝 (𝑡) = 𝑐(𝑡)𝑦ℎ (𝑡)

in terms of the homogeneous solution. A similar approach will be

employed for nonhomogeneous PDEs.

e.g. Consider the vibration problem of a taut string with time-dependent

forcing:
 42

2
1 −𝑡
𝑢𝑡𝑡 = 𝑐 𝑢𝑥𝑥 + 𝑒 , 0 < 𝑥 < 𝐿, 𝑡>0
𝜌

subject to boundary

𝑢(0, 𝑡) = 𝑢(𝐿, 𝑡) = 0, 𝑡 > 0,

and initial conditions

𝑢(𝑥, 0) = 𝑓(𝑥 ), 𝑢𝑡 (𝑥, 0) = 0, 0 < 𝑥 < 𝐿.

In this example, we assume that forcing 𝑒 −𝑡 /𝜌 does not depend on 𝑥

for clarity even though the approach to be presented might as well apply
to the general case of time- and space-dependent nonhomogeneities.

The procedure is as follows:

(I) The solution to the associated homogeneous PDE, i.e. as if the

forcing were absent,

𝑢𝑡𝑡 = 𝑐 2 𝑢𝑥𝑥 , 0 < 𝑥 < 𝐿, 𝑡>0

𝑢(0, 𝑡) = 𝑢(𝐿, 𝑡) = 0, 𝑡 > 0,

𝑢(𝑥, 0) = 𝑓 (𝑥 ), 𝑢𝑡 (𝑥, 0) = 0, 0 < 𝑥 < 𝐿,

would be
∞
𝑘𝜋𝑐 𝑘𝜋
𝑢(𝑥, 𝑡) = ∑ 𝐴𝑘 𝑐𝑜𝑠 𝑡 𝑠𝑖𝑛 𝑥
⏟ 𝐿 𝐿
𝑘=1
𝐶𝑘 (𝑡)

as obtained in a previous example.

 43

(II) In order to account for the time-dependent forcing, we assume a

general lumped time-dependence of the expansion coefficients
∞
𝑘𝜋
𝑢(𝑥, 𝑡) = ∑ 𝐶𝑘 (𝑡) 𝑠𝑖𝑛 𝑥
𝐿
𝑘=1

in the eigenfunction expansion.

Substituting this form of 𝑢(𝑥, 𝑡) into the original nonhomogeneous PDE

yields
∞ ∞
′′
𝑘𝜋 2
𝑘𝜋 2 𝑘𝜋 1 −𝑡
∑ 𝐶𝑘 (𝑡) 𝑠𝑖𝑛 𝑥 = −𝑐 ∑( ) 𝐶𝑘 (𝑡) 𝑠𝑖𝑛 𝑥+ 𝑒
𝐿 𝐿 𝐿 𝜌
𝑘=1 𝑘=1

In order to incorporate the nonhomogeneous term into the series

representations, it is also expanded in terms of the eigenfunctions:
∞
1 −𝑡 𝑘𝜋
𝑒 = ∑ 𝐵𝑘 (𝑡) 𝑠𝑖𝑛 𝑥
𝜌 𝐿
𝑘=1

where

2 𝐿 1 −𝑡 𝑘𝜋 2 −𝑡
𝐵𝑘 (𝑡) = ∫ 𝑒 𝑠𝑖𝑛( 𝑥)𝑑𝑥 = 𝑒 [1 + (−1)𝑘+1 ].
𝐿 0 𝜌 𝐿 𝑘𝜋𝜌

The projected form of the original PDE becomes

∞
𝑐𝑘𝜋 2 𝑘𝜋
∑ {𝐶𝑘′′ (𝑡) + ( 𝐿 )2 𝐶𝑘 (𝑡) − (𝑘𝜋𝜌) 𝑒 −𝑡 [1 + 𝑘+1
(−1) ]} 𝑠𝑖𝑛 𝑥 = 0.
𝐿
𝑘=1
 44

Thus, the original nonhomogeneous PDE reduces to a system of 2nd

order nonhomogeneous ODEs:
𝑐𝑘𝜋 2 2
𝐶𝑘′′ (𝑡) + (𝐿
) 𝐶𝑘 (𝑡) = (𝑘𝜋𝜌) 𝑒 −𝑡 [1 + (−1)𝑘+1 ], 𝑘 = 1,2, …

whose general sol’n is

𝑐𝑘𝜋 𝑐𝑘𝜋
𝐶𝑘 (𝑡) = 𝑎𝑘 𝑐𝑜𝑠( 𝑡) + 𝑏𝑘 𝑠𝑖𝑛( 𝑡)
𝐿 𝐿
2𝐿2 −𝑡 𝑘+1 ]
+ 𝑒 [1 + ( −1) .
𝑘𝜋𝜌(𝐿2 + (𝑐𝑘𝜋)2 )

The homogeneous initial condition yields

∞
𝑘𝜋
𝑢𝑡 (𝑥, 0) = ∑ 𝐶𝑘′ (0) 𝑠𝑖𝑛 𝑥=0 ⇒ 𝐶𝑘′ (0) = 0
𝐿
𝑘=1

that is satisfied when

2𝐿3 [1 + (−1)𝑘+1 ]
𝑏𝑘 =
(𝑘𝜋)2 𝜌𝑐 (𝐿2 + (𝑐𝑘𝜋)2 )

for all 𝑘, so that

𝑐𝑘𝜋
𝐶𝑘 (𝑡) = 𝑎𝑘 𝑐𝑜𝑠( 𝑡)
𝐿
2𝐿2 [1 + (−1)𝑘+1 ] −𝑡
𝑐𝑘𝜋
+ (𝑘𝜋𝑐𝑒 + 𝐿 𝑠𝑖𝑛( 𝑡)).
(𝑘𝜋)2 𝜌𝑐(𝐿2 + 𝑐 2 (𝑘𝜋)2 ) 𝐿

The solution to the original PDE is obtained after substitution

 45

∞
𝑐𝑘𝜋
𝑢(𝑥, 𝑡) = ∑ { 𝑎𝑘 𝑐𝑜𝑠( 𝑡)
𝐿
𝑘=1

2𝐿2 [1 + (−1)𝑘+1 ] −𝑡
𝑐𝑘𝜋 𝑘𝜋
+ (𝑘𝜋𝑐𝑒 + 𝐿 𝑠𝑖𝑛( 𝑡)) } 𝑠𝑖𝑛 𝑥.
(𝑘𝜋)2 𝜌𝑐(𝐿2 + 𝑐 2 (𝑘𝜋)2 ) 𝐿 𝐿

The nonhomogeneous initial condition, on the other hand,

∞
2𝐿2 [1 + (−1)𝑘+1 ] 𝑘𝜋
𝑢(𝑥, 0) = 𝑓(𝑥) = ∑ { 𝑎𝑘 + } 𝑠𝑖𝑛 𝑥
𝑘𝜋𝜌(𝐿2 + 𝑐 2 (𝑘𝜋)2 ) 𝐿
𝑘=1

yields

2𝐿2 [1 + (−1)𝑘+1 ] 2 𝐿 𝑘𝜋
𝑎𝑘 + = ∫ 𝑓(𝑥) 𝑠𝑖𝑛( 𝑥)𝑑𝑥 .
𝑘𝜋𝜌(𝐿2 + 𝑐 2 (𝑘𝜋)2 ) 𝐿 0 𝐿

Now that 𝑎𝑘 and 𝑏𝑘 are determined, the solution procedure is complete.

e.g. Now, consider the temperature distribution in a rod with time-

dependent temperatures at the ends as time-dependent nonhomogeneity
in the boundary conditions.

𝑢𝑡 = 𝜅𝑢𝑥𝑥 , 0 < 𝑥 < 𝐿, 𝑡>0

subject to boundary

𝑢(0, 𝑡) = 𝑈0 (𝑡), 𝑢(𝐿, 𝑡) = 𝑈𝐿 (𝑡), 𝑡 > 0,

and initial conditions

𝑢(𝑥, 0) = 𝑓 (𝑥 ), 0 < 𝑥 < 𝐿.

 46

In this case, we introduce a new dependent variables 𝑣 and 𝜓:

𝑢(𝑥, 𝑡) = 𝑣(𝑥, 𝑡) + 𝜓(𝑥, 𝑡)

in order to have homogeneous boundary conditions for 𝑣(𝑥, 𝑡)

𝑣(0, 𝑡) = 𝑈0 (𝑡) − 𝜓(0, 𝑡) = 0,

𝑣(𝐿, 𝑡) = 𝑈𝐿 (𝑡) − 𝜓(𝐿, 𝑡) = 0.

This is accomplished with the simplest choice of linear blending of

𝑈0 (𝑡) and 𝑈𝐿 (𝑡) in the domain 0 < 𝑥 < 𝐿:
1
𝜓(𝑥, 𝑡) = 𝑈0 (𝑡) + [𝑈𝐿 (𝑡) − 𝑈0 (𝑡)] 𝑥.
𝐿
Note that 𝜓(0, 𝑡) = 𝑈0 (𝑡) and 𝜓(𝐿, 𝑡) = 𝑈𝐿 (𝑡).

Linearity of 𝜓(𝑥, 𝑡) in 𝑥 reduces the original PDE to

=0
⏞
𝑣𝑡 = 𝜅 𝑣𝑥𝑥 − [ 𝜓𝑡 − 𝜅 𝜓𝑥𝑥 ]

′( ) 1 ′( ) ′( )
= 𝜅 𝑣𝑥𝑥 −𝑈
⏟ 0 𝑡 − [𝑈𝐿 𝑡 − 𝑈0 𝑡 ] 𝑥.
𝐿
−𝜓𝑡

Now using the earlier procedure, we proceed to solving 𝑣(𝑥, 𝑡) with

time-dependent nonhomogeneity:

𝑣𝑡 = 𝜅𝑣𝑥𝑥 + 𝐹 (𝑥, 𝑡), 0 < 𝑥 < 𝐿, 𝑡 > 0,

𝑣(0, 𝑡) = 𝑣(𝐿, 𝑡) = 0, 𝑡 > 0,

1
𝑣(𝑥, 0) = 𝑓(𝑥 ) − 𝑈0 (0) − 𝐿 [𝑈𝐿 (0) − 𝑈0 (0)]𝑥, 0 < 𝑥 < 𝐿.
 47

1
where 𝐹 (𝑥, 𝑡) = −𝑈0′ (𝑡) − 𝐿 [𝑈𝐿′ (𝑡) − 𝑈0′ (𝑡)] 𝑥.

Again, we assume a general time-dependence

∞
𝑘𝜋
𝑣(𝑥, 𝑡) = ∑ 𝐶𝑘 (𝑡) 𝑠𝑖𝑛 𝑥
𝐿
𝑘=1

in terms of the eigenfunctions of the associated homogeneous problem

𝑣𝑡 = 𝜅𝑣𝑥𝑥 + 𝐹 (𝑥, 𝑡), 𝑣(0, 𝑡) = 𝑣(𝐿, 𝑡) = 0.

This results in
∞
𝑘𝜋 𝑘𝜋
∑ {𝐶𝑘′ (𝑡) + 𝜅 ( 𝐿 )2 𝐶𝑘 (𝑡) − 𝐹𝑘 (𝑡) } 𝑠𝑖𝑛 𝑥=0
𝐿
𝑘=1

where 𝐹𝑘 (𝑡) are the eigenfunction expansion coefficients of 𝐹(𝑥, 𝑡):

2 𝐿 𝑘𝜋
𝐹𝑘 (𝑡) = ∫ 𝐹(𝑥, 𝑡) 𝑠𝑖𝑛( 𝑥) 𝑑𝑥 .
𝐿 0 𝐿

The general solution of the resulting system of 1st order ODEs

𝑘𝜋
𝐶𝑘′ (𝑡) + 𝜅 ( 𝐿 )2 𝐶𝑘 (𝑡) = 𝐹𝑘 (𝑡)

is then obtained
𝑡
𝑘𝜋 𝑘𝜋
𝐶𝑘 (𝑡) = 𝑎𝑘 𝑒𝑥𝑝( − ( 𝐿 )2 𝜅𝑡) + ∫ 𝐹𝑘 (𝜏) 𝑒𝑥𝑝( − ( 𝐿 )2 𝜅(𝜏 − 𝑡)) 𝑑𝜏
0

by using the techniques of solving 1st order nonhomogeneous linear

ODEs as introduced earlier.
 48

Finally, the initial condition is satisfied by the eigenfunction expansion

∞
1 𝑘𝜋
𝑣(𝑥, 0) = 𝑓(𝑥) − 𝑈0 (0) − [𝑈𝐿 (0) − 𝑈0 (0)]𝑥 = ∑ 𝑎𝑘 𝑠𝑖𝑛 𝑥
𝐿 𝐿
𝑘=1

where

2 𝐿 1 𝑘𝜋
𝑎𝑘 = ∫ {𝑓(𝑥) − 𝑈0 (0) − [𝑈𝐿 (0) − 𝑈0 (0)]𝑥} 𝑠𝑖𝑛( 𝑥)𝑑𝑥
𝐿 0 𝐿 𝐿

using the orthogonality of the eigenfunctiosn.

The assembled final form of the solution is therefore

1
𝑢(𝑥, 𝑡) = 𝑣(𝑥, 𝑡) + 𝑈0 (𝑡) + [𝑈𝐿 (𝑡) − 𝑈0 (𝑡)]𝑥.
𝐿
In summary:

The main objective here is to use techniques above to homogenize the

given problem to achieve separability.

(1) When time-independent nonhomogeneities occur (in PDE or

boundary conditions), separate the static (steady-state) sol’n and solve the
remanining homogeneous PDE with homogeneous boundary conditions.

(2) When time-dependent nonhomogeneities are in the boundary

conditions, transfer the nonhomogeneity into the PDE and proceed with
the eigenfunction expansion reducing the PDE to a system of
nonhomogeneous ODEs in time.
 49

e.g. Consider the (elliptic) Helmholtz eqution

𝑢𝑥𝑥 + 𝑢𝑦𝑦 − 𝑠𝑢 = 0, 0 < 𝑥 < 𝑎, 0<𝑦<𝑏

subject to the nonhomogeneous boundary conditions

𝑢(0, 𝑦) = ℎ(𝑦), 𝑢𝑥 (𝑎, 𝑦) = 𝑟(𝑦), 0 < 𝑦 < 𝑏,

𝑢(𝑥, 0) = 𝑓(𝑥 ), 𝑢(𝑥, 𝑏) = 𝑔(𝑥 ), 0 < 𝑥 < 𝑎.

Nonhomogeneities in this case are easily handled by the superposition

𝑢 =𝑣+𝑤

that shares the nonhomogeneities in the boundary conditions as follows:

𝑣(0, 𝑦) = 0, 𝑣𝑥 (𝑎, 𝑦) = 0, 0 < 𝑦 < 𝑏,

𝑣(𝑥, 0) = 𝑓 (𝑥 ), 𝑣(𝑥, 𝑏) = 𝑔(𝑥 ), 0 < 𝑥 < 𝑎.

and

𝑤(0, 𝑦) = ℎ(𝑦), 𝑤𝑥 (𝑎, 𝑦) = 𝑟(𝑦), 0 < 𝑦 < 𝑏,

𝑤(𝑥, 0) = 0, 𝑤(𝑥, 𝑏) = 0, 0 < 𝑥 < 𝑎.

 50

thus, leaving one set of homogeneous boundary conditions for each

dependent variable.

For solving 𝑣𝑥𝑥 + 𝑣𝑦𝑦 − 𝑠𝑣 = 0, we proceed with 𝑣(𝑥, 𝑦) = 𝑋(𝑥)𝑌(𝑦)

to get the ODEs:

(1) 𝑋 ″ − 𝜆𝑋 = 0 with 𝑋(0) = 𝑋 ′ (𝑎) = 0,

(2) 𝑌 ″ + (𝜆 − 𝑠)𝑌 = 0.

Now, (1) forms a regular S-L problem and for 𝜆 = −𝛼 2 it yields the
eigensolutions:

𝛼𝑘 = ((2𝑘 − 1)𝜋)/2𝑎 & 𝑋𝑘 (𝑥) = 𝑠𝑖𝑛( 𝛼𝑘 𝑥), 𝑘 = 1,2, …

(2) yields the solution:

𝑌𝑘 (𝑦) = 𝐴𝑘 𝑐𝑜𝑠ℎ (√(𝑠 + 𝛼𝑘2 ) 𝑦) + 𝐵𝑘 𝑠𝑖𝑛ℎ (√(𝑠 + 𝛼𝑘2 ) 𝑦).

The form of the solutions are then

𝑣𝑘 (𝑥, 𝑦) = 𝑋𝑘 (𝑥 )𝑌𝑘 (𝑦) =

= [ 𝐴𝑘 𝑐𝑜𝑠ℎ (√(𝑠 + 𝛼𝑘2 ) 𝑦) + 𝐵𝑘 𝑠𝑖𝑛ℎ (√(𝑠 + 𝛼𝑘2 ) 𝑦) ] 𝑠𝑖𝑛( 𝛼𝑘 𝑥).

The superposed solution

∞

𝑣(𝑥, 𝑦) = ∑ 𝑣𝑘 (𝑥, 𝑦)
𝑘=1
 51

is now used to satisfy the boundary conditions:

∞

𝑣(𝑥, 0) = 𝑓(𝑥 ) = ∑ 𝐴𝑘 𝑠𝑖𝑛( 𝛼𝑘 𝑥) ,

𝑘=1

and

𝑣(𝑥, 𝑏) = 𝑔(𝑥 ) =
∞

= ∑ [𝐴𝑘 𝑐𝑜𝑠ℎ (√(𝑠 + 𝛼𝑘2 ) 𝑏) + 𝐵𝑘 𝑠𝑖𝑛ℎ (√(𝑠 + 𝛼𝑘2 ) 𝑏)] 𝑠𝑖𝑛( 𝛼𝑘 𝑥) .

𝑘=1

Using the orthogonality of the SL-eigenfunctions, we get

2 𝑎
𝐴𝑘 = ∫ 𝑓(𝑥) 𝑠𝑖𝑛( 𝛼𝑘 𝑥)𝑑𝑥
𝑎 0

and

1 2 𝑎
𝐵𝑘 = ( ∫ 𝑔(𝑥) 𝑠𝑖𝑛( 𝛼𝑘 𝑥)𝑑𝑥 − 𝐴𝑘 𝑐𝑜𝑠ℎ (√(𝑠 + 𝛼2𝑘 ) 𝑏)).
𝑠𝑖𝑛ℎ(√(𝑠+𝛼2 ) 𝑏) 𝑎 0
𝑘

Verify that for the selection of 𝑓(𝑥) = 0 and 𝑔(𝑥) = 𝑥, it yields

𝐴𝑘 = 0

and
 52

1 2 𝑎
𝐵𝑘 = ( ∫ 𝑥 𝑠𝑖𝑛( 𝛼𝑘 𝑥)𝑑𝑥 )
𝑠𝑖𝑛ℎ(√(𝑠+𝛼𝑘2 ) 𝑏)
𝑎 0

1 4√2𝑎3 (−1)𝑘+1
= .
𝑠𝑖𝑛ℎ(√(𝑠+𝛼𝑘2 ) 𝑏)
(2𝑘 − 1)2 𝜋 2

Similarly, one can proceed to solve

𝑤𝑥𝑥 + 𝑤𝑦𝑦 − 𝑠𝑤 = 0

for 𝑤 = 𝑤(𝑥, 𝑦) with the associated boundary conditions to complete

the solution 𝑢(𝑥, 𝑦) = 𝑣(𝑥, 𝑦) + 𝑤(𝑥, 𝑦).

Problems with Multi-Space Variables

This is to demonstrate the use of Separation of Variables technique to

problems in two or more space variables.

e.g. Consider the vertical oscillations of a rectangular membrane model:

𝑢𝑡𝑡 = 𝑐 2 (𝑢𝑥𝑥 + 𝑢𝑦𝑦 ), 0 < 𝑥 < 𝑎, 0 < 𝑦 < 𝑏, 𝑡 > 0

subject to the boundary conditions:

𝑢(0, 𝑦, 𝑡) = 𝑢(𝑎, 𝑦, 𝑡) = 0, 0 < 𝑦 < 𝑏, 𝑡 > 0,

𝑢(𝑥, 0, 𝑡) = 𝑢(𝑥, 𝑏, 𝑡) = 0, 0 < 𝑥 < 𝑎, 𝑡 > 0,

and the initial conditions:

𝑢(𝑥, 𝑦, 0) = 𝑓(𝑥, 𝑦), 𝑢𝑡 (𝑥, 𝑦, 0) = 0, 0 < 𝑥 < 𝑎, 0 < 𝑦 < 𝑏.

 53

We proceed again with 𝑢(𝑥, 𝑦, 𝑡) = 𝑋(𝑥)𝑌(𝑦)𝑇(𝑡) to get

𝑋″ 𝑌″ 𝑇″
=− + 2 =𝜆
𝑋 𝑌 𝑐 𝑇
and identify

(1) 𝑋 ″ − 𝜆𝑋 = 0 with 𝑋(0) = 𝑋(𝑎) = 0

that for 𝜆 = −𝛼 2 yields the eigensolutions

𝛼𝑛 = 𝑛𝜋/𝑎 and 𝑋𝑛 (𝑥) = 𝑠𝑖𝑛( 𝛼𝑛 𝑥).

Now, continuing with

𝑌″ 𝑇″
= 2 + 𝛼𝑛2 = 𝜇
𝑌 𝑐 𝑇
to identify

(2) 𝑌 ″ − 𝜇𝑌 = 0 with 𝑌(0) = 𝑌(𝑏) = 0

that for 𝜇 = −𝛽 2 yields the eigensolutions

𝛽𝑚 = 𝑚𝜋/𝑏 and 𝑌𝑚 (𝑦) = 𝑠𝑖𝑛( 𝛽𝑚 𝑦).

Finally obtaining

(3) 𝑇 ″ + 𝑐 2 (𝛼𝑛2 + 𝛽𝑚
2
)𝑇 = 0 with 𝑇 ′ (0) = 0

yielding the solution

𝑇𝑚𝑛 (𝑡) = 𝐴𝑚𝑛 𝑐𝑜𝑠( 𝑐√(𝛼𝑛2 + 𝛽𝑚

2
) 𝑡).
 54

The superposition
∞ ∞

𝑢(𝑥, 𝑦, 𝑡) = ∑ ∑ 𝑢𝑚𝑛 (𝑥, 𝑦, 𝑡)

𝑚=1 𝑛=1
∞ ∞

= ∑ ∑ 𝐴𝑚𝑛 𝑐𝑜𝑠( 𝑐 √(𝛼𝑛2 + 𝛽𝑚

2
) 𝑡) 𝑠𝑖𝑛( 𝛼𝑛 𝑥) 𝑠𝑖𝑛( 𝛽𝑚 𝑦)
𝑚=1 𝑛=1

is required to satisfy the nonhomogeneous initial condition

∞ ∞

𝑢(𝑥, 𝑦, 0) = 𝑓(𝑥, 𝑦) = ∑ (∑ 𝐴𝑚𝑛 𝑠𝑖𝑛( 𝛼𝑛 𝑥)) 𝑠𝑖𝑛( 𝛽𝑚 𝑦)

𝑚=1 𝑛=1

where
∞
2 𝑏
∑ 𝐴𝑚𝑛 𝑠𝑖𝑛( 𝛼𝑛 𝑥) = ∫ 𝑓(𝑥, 𝑦) 𝑠𝑖𝑛( 𝛼𝑛 𝑥)𝑑𝑦
𝑏 0
𝑛=1

and thus

2 𝑎 2 𝑏
𝐴𝑚𝑛 = ∫ ( ∫ 𝑓(𝑥, 𝑦) 𝑠𝑖𝑛( 𝛽𝑚 𝑦)𝑑𝑦) 𝑠𝑖𝑛( 𝛼𝑛 𝑥)𝑑𝑥 .
𝑎 0 𝑏 0

For 𝑓(𝑥, 𝑦) = 𝑥𝑦(𝑎 − 𝑥)(𝑏 − 𝑦), verify that

16𝑎2 𝑏2
𝐴𝑚𝑛 = 3 3 6 [1 + (−1)𝑛+1 ][1 + (−1)𝑚+1 ].
𝑛 𝑚 𝜋
 55

Normal Modes of Vibration

𝑢𝑛=1,𝑚=1 (𝑥, 𝑦, 0)

𝑢𝑛=1,𝑚=3 (𝑥, 𝑦, 0)

ex. Verify that conduction temperature distribution on a slab:

𝑢𝑡 = 𝜅(𝑢𝑥𝑥 + 𝑢𝑦𝑦 ), 0 < 𝑥 < 𝑎, 0 < 𝑦 < 𝑏, 𝑡>0

subject to the boundary conditions:

𝑢(0, 𝑦, 𝑡) = 𝑢(𝑎, 𝑦, 𝑡) = 0, 0 < 𝑦 < 𝑏, 𝑡 > 0,

𝑢(𝑥, 0, 𝑡) = 𝑢(𝑥, 𝑏, 𝑡) = 0, 0 < 𝑦 < 𝑎, 𝑡 > 0,

and the initial conditions:

 56

𝑢(𝑥, 𝑦, 0) = 𝑓 (𝑥, 𝑦), 0 < 𝑥 < 𝑎, 0 < 𝑦 < 𝑏,

is given by
∞ ∞

𝑢(𝑥, 𝑦, 𝑡) = ∑ ∑ 𝐴𝑚𝑛 𝑒𝑥𝑝( − 𝜅(𝛼𝑛2 + 𝛽𝑚

2
)𝑡) 𝑠𝑖𝑛( 𝛼𝑛 𝑥) 𝑠𝑖𝑛( 𝛽𝑚 𝑦)
𝑚=1 𝑛=1

where 𝛼𝑛 = 𝑛𝜋/𝑎, 𝛽𝑚 = 𝑚𝜋/𝑏, and

2 𝑎 2 𝑏
𝐴𝑚𝑛 = ∫ ( ∫ 𝑓(𝑥, 𝑦) 𝑠𝑖𝑛( 𝛽𝑚 𝑦)𝑑𝑦) 𝑠𝑖𝑛( 𝛼𝑛 𝑥)𝑑𝑥 .
𝑎 0 𝑏 0

ex. Verify that the solution to the Laplace (equilibrium) problem:

𝑢𝑥𝑥 + 𝑢𝑦𝑦 + 𝑢𝑧𝑧 = 0, 0 < 𝑥 < 𝑎, 0 < 𝑦 < 𝑏, 0 < 𝑧 < 𝑐,

subject to the boundary conditions

𝑢𝑥 (0, 𝑦, 𝑧) = 𝑢(𝑎, 𝑦, 𝑧) = 0, 0 < 𝑦 < 𝑏, 0 < 𝑧 < 𝑐,

𝑢(𝑥, 0, 𝑧) = 𝑢𝑦 (𝑥, 𝑏, 𝑧) = 0, 0 < 𝑥 < 𝑎, 0 < 𝑧 < 𝑐,

𝑢(𝑥, 𝑦, 0) = 𝑓(𝑥, 𝑦), 0 < 𝑥 < 𝑎, 0 < 𝑦 < 𝑏,

 57

𝑢(𝑥, 𝑦, 𝑐 ) = 𝑔(𝑥, 𝑦), 0 < 𝑥 < 𝑎, 0 < 𝑦 < 𝑏,

is given by
∞ ∞

𝑢(𝑥, 𝑦, 𝑧) = ∑ ∑ {𝐴𝑚𝑛 𝑐𝑜𝑠ℎ (√(𝛼𝑛2 + 𝛽𝑚

2
) 𝑧)
𝑚=0 𝑛=0

+ 𝐵𝑚𝑛 𝑠𝑖𝑛ℎ (√(𝛼𝑛2 + 𝛽𝑚

2
) 𝑧)} 𝑐𝑜𝑠( 𝛼𝑛 𝑥) 𝑠𝑖𝑛( 𝛽𝑚 𝑦),

where 𝛼𝑛 = (2𝑛 + 1)𝜋/2𝑎, 𝛽𝑚 = (2𝑚 + 1)𝜋/2𝑏,

2 𝑎 2 𝑏
𝐴𝑚𝑛 = ∫ ( ∫ 𝑓(𝑥, 𝑦) 𝑠𝑖𝑛( 𝛽𝑚 𝑦)𝑑𝑦) 𝑐𝑜𝑠( 𝛼𝑛 𝑥)𝑑𝑥
𝑎 0 𝑏 0

and

𝐵𝑚𝑛 =
1 2 𝑎 2 𝑏
= [ ∫ ( ∫ 𝑔(𝑥, 𝑦) 𝑠𝑖𝑛( 𝛽𝑚 𝑦)𝑑𝑦) 𝑐𝑜𝑠( 𝛼𝑛 𝑥)𝑑𝑥
𝑠𝑖𝑛ℎ (√(𝛼 + 𝛽 ) 𝑐) 𝑎 0 𝑏 0
2 2
𝑛 𝑚

− 𝐴𝑚𝑛 𝑐𝑜𝑠ℎ (√(𝛼𝑛2 + 𝛽𝑚

2
) 𝑐) ].