Partial Differential UCL Engineering

Teaching material
Equations Dr Vasos Pavlika

ENGF0004
Mathematical Modelling and
Analysis II

Dr Vasos Pavlika

This booklet was created within the Teaching and Learning Support Summer Studentship, by
courtesy of University College London’s Engineering Faculties, in 2021.

Supervising staff: Student:

Dr Vasos Pavlika, Darius-Marian Filip,
Associate Professor (Teaching): Undergraduate, 2nd Year,
Engineering Mathematics and Statistics. MEng Civil Engineering.

Page | 1
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Table of Contents
Introduction .................................................................................................................................................................3
1. Important Equations Involving PDEs ............................................................................................... 4
2. PDEs vs. ODEs ............................................................................................................................................ 4
3. Order of PDEs ............................................................................................................................................. 4
4. Linear PDEs ................................................................................................................................................. 5
Classification of Second Order Linear PDEs .....................................................................................5

1. Worksheet Questions – Tutorial 1 ..................................................................................................... 8

Boundary Conditions ........................................................................................................................................10

1. Dirichlet Conditions ...............................................................................................................................11

2. Neumann (or Flux) Conditions ..........................................................................................................11
3. Mixed (or Robin) Conditions ..............................................................................................................12
4. Cauchy Conditions ..................................................................................................................................12
5. Worksheet Questions – Tutorial 2 ...................................................................................................13
Separation of Variables ..................................................................................................................................14

1. Worksheet Questions – Tutorial 3 ....................................................................................................20

Numerical Solutions to PDEs ................................................................................................. 24

1. Solution Technique ................................................................................................................................26

2. Worksheet Questions – Tutorial 4 ...................................................................................................31
Worked Examples of PDEs ...........................................................................................................................35

1. ENGF0004 Coursework I: Series and Transform .......................................................................35

2. ENGF0004 Coursework, 2017: Block 2—PDEs ..........................................................................51
Appendix: Fourier Series ....................................................................................................................59

Page | 2
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Page | 3
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Introduction
In the following sections, we will be covering the following topics:
❖ Classification of Partial Differential ❖ Numerical techniques
equations (PDEs) ❖ Crank-Nicolson Technique
❖ Separation of variables

The general form of a PDE is:

𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕 2 𝑢 𝜕 2 𝑢 𝜕3𝑢
𝐹 ቆ𝑥, 𝑦, 𝑧, 𝑡, … , , , , , , ,…, 3 ,… ቇ = 0
𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑡 𝜕𝑥𝜕𝑦 𝜕𝑥𝜕𝑧 𝜕𝑥

Equation 1

What we will be looking at:

❖ First and second order PDEs ❖ Fluid flow

❖ Equations involving PDEs ❖ Elasticity
❖ Electrostatics (Maxwell’s Equations) ❖ Other areas of engineering

You should avoid using these formulations:

𝜕𝑢 𝜕2 𝑢 𝜕2 𝑢
𝑢𝑥 = 𝜕𝑥 ; 𝑢𝑥𝑦 = 𝜕𝑥𝜕𝑦 ; 𝑢𝑦𝑦 = 𝜕𝑣 2 .

Introducing the Laplacian (or the Laplace operator):

𝜕2 𝑢 𝜕2 𝑢 𝜕2 𝑢
∇2 𝑢 = ቀ𝜕𝑥 2 + 𝜕𝑦2 + 𝜕𝑧2 , ቁ; where u is the dependent variable

Equation 2

Important Equations involving PDEs

❖ The 3D Heat Equation ❖ The Laplace Equation

𝜕𝑢 𝜕2 𝑢 𝜕2 𝑢 𝜕2 𝑢
𝜕𝑡
= 𝑅 ቀ𝜕𝑥 2 + 𝜕𝑦2 + 𝜕𝑧2 , ቁ = R∇2 𝑢 ∇2 𝑢 = 0

Equation 3 Equation 4
Page | 4
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

❖ The Poisson Equation ❖ Linear 3D Wave Equation

∇2 𝑢 = 𝜎 ∂2 u
= 𝑐 2 ∇2 𝑢; 𝑤ℎ𝑒𝑟𝑒 𝑐 = 𝑠𝑝𝑒𝑒𝑑 𝑜𝑓 𝑙𝑖𝑔ℎ𝑡
∂t 2
Equation 5 Equation 6

❖ Non-linear Burgers’ Equation

∂u ∂u ∂2 u The term 𝑢 is the dependent variable and its

+𝑢 =𝜇 2 presence in the term 𝜕𝑢/𝜕𝑥 makes the Burgers’
∂t ∂x ∂x
Equation a non-linear equation.
Equation 7

PDEs vs. ODEs

In the case of Ordinary Differential Equations (ODEs), if the derivative of a function u (which is only a function
of x) with respect to x is equal to zero, u must be a constant. If the partial derivative of a function u (this time
a function of x and y) with respect to x is equal to zero, u could be a function of y.
This means that the integration of a partial derivative results in an arbitrary function of all other variables.

PDE ODE

∂𝑢(𝑥,𝑦) 𝑑𝑢
=0 =0
∂𝑥 𝑑𝑥

u(x,y) = f(y) u(x) = c

𝜕𝑓(𝑦)
∴ =0
𝜕𝑥

Similarly,
𝜕𝑢(𝑥,𝑦,𝑧)
=0
𝜕𝑥

u(x,y,z) = f(y,z)

Example

Order of PDEs
The order of a PDE is the highest order derivative of the dependent variable (here we have used u).

∂u ∂3 u
∂t
= 𝑢 ∂x3 + sin (𝑥) - 3rd order PDE

Example

Page | 5
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Linear Partial Differential Equations (PDEs)

A PDE is said to be linear if the dependent variable only occurs in a linear form, i.e.:

𝜕𝑢 𝜕𝑢
𝑎(𝑥, 𝑦) 𝜕𝑥 + 𝑏(𝑥, 𝑦) 𝜕𝑦 + 𝑐(𝑥, 𝑦)𝑢 = 𝑑(𝑥, 𝑦);

a≡ 𝑎(𝑥, 𝑦); 𝑏 ≡ 𝑏(𝑥, 𝑦); 𝑐 ≡ 𝑐(𝑥, 𝑦).

Equation 8
Almost linear PDEs:

𝜕𝑢 𝜕𝑢
𝑃(𝑥, 𝑦) + 𝑄(𝑥, 𝑦) = 𝑅(𝑥, 𝑦, 𝑢)
𝜕𝑥 𝜕𝑦

Equation 9
The right-hand side of the equation can be a function of u.

Quasi-linear PDEs:

𝜕𝑢 𝜕𝑢
𝑃(𝑥, 𝑦, 𝑢) + 𝑄(𝑥, 𝑦, 𝑢) = 𝑅(𝑥, 𝑦, 𝑢)
𝜕𝑥 𝜕𝑦

Equation 10

Classification of Second Order Linear PDEs

Consider the general 2nd order linear PDE:

𝜕 2𝑢 𝜕 2𝑢 𝜕 2𝑢 𝜕𝑢 𝜕𝑢
𝐴 2
+ 𝑩 + 𝐶 2
+𝐷 +𝐸 + 𝐹𝑢 = 𝐺(𝑥, 𝑦)
𝜕𝑥 𝜕𝑥𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑦
Some authors write this as:
𝜕 2𝑢 𝜕 2𝑢 𝜕 2𝑢 𝜕𝑢 𝜕𝑢
𝐴 2
+ 𝟐𝑩 + 𝐶 2
+𝐷 +𝐸 + 𝐹𝑢 = 𝐺(𝑥, 𝑦)
𝜕𝑥 𝜕𝑥𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑦

Equation 11

The classification is determined by the so-

called Discriminant (reminding us of the 𝐷 = 𝑩𝟐 − 4𝐴𝐶 𝑜𝑟 𝐷 = ൫𝟐𝑩)𝟐 − 4𝐴𝐶 = 4(𝐵2 − 𝐴𝐶൯
quadratic inversion formula):
𝑜𝑟 𝐷 = 𝐵2 − 𝐴𝐶
Equation 12
Page | 6
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

PDEs can be classified into 3 categories, depending on the value of the discriminant:

< 𝟎 , 𝑡ℎ𝑒 𝑃𝐷𝐸 𝑖𝑠 𝑬𝒍𝒍𝒊𝒑𝒕𝒊𝒄

If the discriminant is D = 𝟎 , 𝑡ℎ𝑒 𝑃𝐷𝐸 𝑖𝑠 𝑷𝒂𝒓𝒂𝒃𝒐𝒍𝒊𝒄

{ > 𝟎 , 𝑡ℎ𝑒 𝑃𝐷𝐸 𝑖𝑠 𝑯𝒚𝒑𝒆𝒓𝒃𝒐𝒍𝒊𝒄

These three geometric plane curves represent conic sections,

and they can be visualized in Figure 1, by “slicing” a cone with
various planes.
You can also associate these three curves with their Greek
meaning:

❖ Ellipse – to be less than

❖ Para – to be equal to
❖ Hyper – to be more than

Figure 1 – Conic sections

Now let us look at the following examples:

𝜕 2𝑢 𝜕 2𝑢 𝜕 2 𝑢 𝜕𝑢
+ 2𝑦 + 𝑥 − +𝑢 =0
𝜕𝑥 2 𝜕𝑥𝜕𝑦 𝜕𝑦 2 𝜕𝑥
D = 𝐵 2 − 4𝐴𝐶; 𝐴 = 1; 𝐵 = 2𝑦; 𝐶 = 𝑥;
∴ 𝐷 = 4𝑦 2 − 4𝑥
= 4(𝑦 2 − 𝑥)
The value of the discriminant depends on the value of y, so there are 3 possible cases:

𝒚𝟐 < 𝒙 , 𝑡ℎ𝑒 𝑃𝐷𝐸 𝑖𝑠 𝑒𝑙𝑙𝑖𝑝𝑡𝑖𝑐

if 𝒚𝟐 = 𝒙 , 𝑡ℎ𝑒 𝑃𝐷𝐸 𝑖𝑠 𝑝𝑎𝑟𝑎𝑏𝑜𝑙𝑖𝑐

{ 𝒚𝟐 > 𝒙 , 𝑡ℎ𝑒 𝑃𝐷𝐸 𝑖𝑠 ℎ𝑦𝑝𝑒𝑟𝑏𝑜𝑙𝑖𝑐

Example

Page | 7
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

𝜕 2𝑢 𝜕 2𝑢 𝜕𝑢
2𝑥𝑦 +𝑥 +𝑦 +𝑢 =0
𝜕𝑥𝜕𝑦 𝜕𝑦 𝜕𝑥
D = 𝐵 2 − 4𝐴𝐶; 𝐵 = 2𝑥𝑦; 𝐴 = 𝐶 = 0;
∴ 𝐷 = 4𝑥 2 𝑦 2
∴ 𝐷 ≥ 0, therefore, the PDE is hyperbolic for x, y ≠ 0.
At the origin or x = 0, y = 0 the equation is of parabolic type.

Example

𝜕 2𝑢 𝜕 2𝑢 𝜕 2𝑢 𝜕𝑢
3 2
+ 2 + 5 2
+𝑥 =0
𝜕𝑥 𝜕𝑥𝜕𝑦 𝜕𝑦 𝜕𝑦
D = 𝐵 2 − 4𝐴𝐶; 𝐴 = 3; 𝐵 = 2; 𝐶 = 5;
∴ 𝐷 = 22 − 4(3)(5)
∴ 𝐷 = −56 < 0, ℎ𝑒𝑛𝑐𝑒 𝑒𝑙𝑙𝑖𝑝𝑡𝑖𝑐

Example

It is crucial to state the form of the PDE and the discriminant that you will be using from the
beginning, and stick to it, in order to avoid calculation mistakes. In the previous examples we
have used the first formula:

𝜕 2𝑢 𝜕 2𝑢 𝜕 2𝑢 𝜕𝑢 𝜕𝑢
𝐴 + 𝐵 + 𝐶 + 𝐷 + 𝐸 + 𝐹𝑢 = 𝐺(𝑥, 𝑦)
𝜕𝑥 2 𝜕𝑥𝜕𝑦 𝜕𝑦 2 𝜕𝑥 𝜕𝑦
𝐷 = 𝐵𝟐 − 4𝐴𝐶

Page | 8
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Worksheet Questions – Tutorial 1

In all questions take as the standard form of a PDE:

𝝏𝟐 𝒖 𝝏𝟐 𝒖 𝝏𝟐 𝒖 𝝏𝒖 𝝏𝒖
𝑨 𝟐
+ 𝟐𝑩 +𝑪 = 𝒇(𝒙, 𝒚, 𝒖, , ) , where f denotes a continuous function of the
𝝏𝒙 𝝏𝒙𝝏𝒚 𝝏𝒚𝟐 𝝏𝒙 𝝏𝒚
arguments shown.

1. Classify the following partial differential equations as to whether they are elliptic, parabolic or
hyperbolic, also state their order and whether they are linear or non-linear, where:

2
𝜕2 𝜕2
∇ ≡ 2+ 2
𝜕𝑥 𝜕𝑦

is the Laplacian operator in 2 dimensions.

i) The wave equation:

𝜕 2𝑢 1 𝜕 2𝑢
=
𝜕𝑥 2 𝑐 2 𝜕𝑡 2

ii) The Laplace equation

∇2 𝑢 = 0

iii) The Poisson equation

∇2 𝑢 = 𝑓(𝑥, 𝑦)

iv) The Heat (or diffusion equation), where t is time

𝜕 2 𝑢 1 𝜕𝑢
=
𝜕𝑥 2 𝑘 𝜕𝑡

v) The Helmholtz’ equation, where k is a constant

∇2 𝑢 + 𝑘𝑢 = 0

vi) Fick’s 2nd law of diffusion, where D is a constant and C is the concentration.

𝜕𝐶 𝜕 2𝐶
=𝐷 2
𝜕𝑡 𝜕𝑥

Page | 9
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

𝜕2𝑢 𝜕2𝑢 𝜕2𝑢 𝜕𝑢

vii) 3 𝜕𝑥 2 + 𝜕𝑥𝜕𝑦 + 5 𝜕𝑦 2 + 𝑥 𝜕𝑦 = 0

viii) The Tricomi equation for transonic flow

𝜕2𝑢 𝜕2𝑢
+ 𝑦 𝜕𝑦 2 =0
𝜕𝑥 2

2. Classify the following partial differential equations (as elliptic, parabolic or hyperbolic etc.)

𝜕2𝑢 𝜕2𝑢 𝜕2𝑢

i. − 𝜕𝑥𝜕𝑦 + 𝜕𝑦 2 = 0
𝜕𝑥 2
𝜕2𝑢 𝜕2𝑢 𝜕𝑢
ii. + 𝜕𝑦 2 + 4 𝜕𝑥 = 0
𝜕𝑥 2
𝜕2𝑢 𝜕2𝑢 𝜕𝑢
iii. − 𝜕𝑦 2 + 2 𝜕𝑥 = 0
𝜕𝑥 2
𝜕2𝑢 𝜕𝑢 𝜕𝑢
iv. + +3 =0
𝜕𝑥 2 𝜕𝑥 𝜕𝑦
𝜕2 𝑢 𝜕2𝑢
v. 𝑦 𝜕𝑥 2 + 𝑥 𝜕𝑦 2 = 0
𝜕2𝑢 𝜕2𝑢 𝜕2𝑢 𝜕𝑢
vi. 3 𝜕𝑥 2 + 2 𝜕𝑥𝜕𝑦 + 5 𝜕𝑦 2 + 2 𝜕𝑦 = 0

3. State whether the following equations are linear, non-linear or quasi-linear

𝜕𝑢 𝜕𝑢
i. + 𝜕𝑦 = 0
𝜕𝑥
𝜕𝑢 𝜕𝑢 2
ii. +ቀ ቁ =0
𝜕𝑥 𝜕𝑦
𝜕𝑢 𝜕𝑢
iii. + 𝜕𝑦 + 𝑢2 = 0
𝜕𝑥
𝜕2𝑢 𝜕2𝑢
iv. + 𝜕𝑦 2 = 𝑥
𝜕𝑥 2
𝜕2𝑢 𝜕2 𝑢
v. + 𝑢 𝜕𝑦 2 = 0
𝜕𝑥 2

4. What is the order and degree of the of the following partial differential equations?
i.
3
𝜕 2𝑢 𝜕𝑢 𝜕𝑢
ቆ 2ቇ + 3 − =𝑥
𝜕𝑥 𝜕𝑥 𝜕𝑦

ii.

𝜕 2𝑢 𝜕𝑢 2 𝜕𝑢
+ ( ) + 𝑢 = 𝑦sin (𝑥)
𝜕𝑥 2 𝜕𝑥 𝜕𝑦

iii.

Page | 10
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

5
𝜕2𝑢 𝜕𝑢 6 𝜕𝑢
cos (𝑥) ቀ𝜕𝑥 2 ቁ + 2 ቀ𝜕𝑥 ቁ + 𝜕𝑦 = tan (𝑥)

Page | 11
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Boundary Conditions (BCs)

In this section you will be introduced to the four key boundary conditions for PDEs:
❖ Dirichlet Conditions ❖ Mixed (or Robin) Conditions
❖ Neumann (or Flux) Conditions ❖ Cauchy Conditions
Let us consider the following linear second order operator, L:

𝜕2 𝜕2 𝜕2 𝜕 𝜕
𝐿 =𝐴 2+𝐵 +𝐶 2+𝐷 +𝐸 +𝐹
𝜕𝑥 𝜕𝑥𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑦
𝑜𝑟 𝐿[𝑢] = 0

Equation 13

We will use the example below to explain this concept in more detail.

L[u] = 0 on some domain Ω. Find the solution to:

𝛻 2𝑢 = 0 on Ω = {0 ≤ 𝑥 ≤ 𝑎, 0 ≤ 𝑦 ≤ 𝑏 }

Figure 2 –
Boundaries
of a PDE

We refer to dΩ (or 𝝏𝛀) as the boundary of the region Ω. This notation does NOT denote a derivative!

Example

Page | 12
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Dirichlet Conditions
If u(x,y) is given on the boundary, we call such conditions Dirichlet Conditions.

u(x,y) = g(x,y) on dΩ

u(x,y) = 2x + 5y
0≤𝑦≤𝑎
x=0

u(x,y) = 3𝑦 2
x=a Figure 3 –
Dirichlet
0≤𝑦≤𝑏 Conditions

Example

Neumann (or Flux) Conditions

The normal derivative of u is given on the boundary i.e. on dΩ. These conditions are known as Neumann.

𝜕𝑢
= 𝑓(𝑥, 𝑦) 𝑜𝑛 𝑑Ω; where n is the normal derivative.
𝜕𝑛

The normal derivative on each one of the four boundaries of the region Ω are given below:

Figure 4 –
Neumann
Conditions

Example
Page | 13
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Mixed (or Robin) Conditions

This is a linear combination of the Dirichlet and Neumann conditions on the boundary 𝒅𝜴.

𝜕𝑢
𝛼𝑢 + 𝛽 = 𝑓(𝑥, 𝑦) 𝑜𝑛 𝑑Ω; where n is the normal derivative.
𝜕𝑛

Figure 5 –
Mixed
Conditions

Example

Cauchy Conditions

𝝏𝒖
𝝏𝒏
are prescribed on the boundary 𝒅𝛀.
𝝏𝒖 y
It should be noted that are some
𝝏𝒏
b
functions on the boundary.

Figure 6 –
Cauchy 0 a x
Conditions

Example

Page | 14
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Worksheet Questions – Tutorial 2

1. Write the appropriate auxiliary conditions that can be used to solve the problem described below.
Provide detailed explanations to support your answers and classify the auxiliary conditions you write
according to their type (Dirichlet/Neumann/Cauchy/Robin – choose appropriately).

A metal rod of length L (shown below) has an initial temperature throughout its whole length equal to
f(x). At time 0 (t=0), the left end is exposed to a constant temperature of 0°𝐶 while the right end is kept
insulated (i.e. no temperature gradient exists). The PDE that describes the change of temperature (T) as
a function of time (t) and position along the length of the rod (x) is given below. State the boundary
conditions and type that appropriately describe this problem.

 2T T
k =
x 2 t

2. A thin rectangular plate (shown below) rests in the region defined by 0 ≤ x ≤ 4 and 0 ≤ y ≤ 2. The left end
and the bottom of the plate are insulated. The top of the plate is held at temperature zero and the right
end of the plate is held at temperature f(y). The steady state PDE that describes how temperature (T)
changes along the position of the plate (x,y) is given below. State the boundary conditions and type that
appropriately describe the above problem.

 2T  2T
+ =0
x 2 y 2

Page | 15
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Separation of Variables
We will use the Heat Equation to explain the procedure of implementing the separation of variables. The
first step is to express a PDE as a pair of ODEs.
Consider the Heat Equation:

𝜕𝑢 𝐾𝜌 𝜕 2 𝑢
= × 2
𝜕𝑡 𝐶 𝜕𝑥
𝐾 = thermal conductivity; 𝜌 = density; 𝐶 = specific heat capacity;
𝐾𝜌
If we define 𝛼 2 = , this can be rewritten as:
𝐶

𝝏𝒖 𝝏𝟐 𝒖
∴ = 𝜶𝟐 𝟐
𝝏𝒕 𝝏𝒙
In 3D this would become:
𝜕𝑢
∴ = 𝛼2𝛻 2𝑢
𝜕𝑡
Equation 14

We are going to look at the following example: a one-dimensional heat conduction equation, in the x direction
(there is no y dimension). This is represented by the thin wire in the figure below. In order for this to be a
well-posed problem, we also need some boundary conditions to get arbitrary constants and functions
calculated, as they arise.

Ω = {0 ≤ 𝑥 ≤ L, t < ∞ }

Boundary Conditions
u(0,t) = u(L, t) = 0 (Dirichlet Conditions)

Alternatively,
𝜕(0,𝑡) 𝜕(𝐿,𝑡) Figure 7 – Electric wire
= = 0 (Neumann Conditions)
𝜕𝑥 𝜕𝑥

Furthermore, suppose also that we know the temperature distribution at 𝑡 = 0 such that

𝑢(𝑥, 0) = 𝑓(𝑥)

Example
Page | 16
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

The method of separation of variables is based on an ansatz. In this case, we will assume that
u(x,t)=X(x)T(t).

We are given the following PDE:

𝜕𝑢 𝜕 2𝑢
= 𝛼2 2 𝑎𝑠𝑠𝑢𝑚𝑒 𝑢(𝑥, 𝑡) = 𝑋(𝑥)𝑇(𝑡) (𝑜𝑟 𝑢 = 𝑋𝑇 𝑓𝑜𝑟 𝑠𝑖𝑚𝑝𝑙𝑖𝑐𝑖𝑡𝑦)
𝜕𝑡 𝜕𝑥

𝜕𝑋𝑇 𝜕 2 𝑋𝑇 𝜕𝑇 𝜕 2𝑋
= 𝛼2 𝑒𝑞𝑢𝑖𝑣𝑎𝑒𝑛𝑡 𝑡𝑜 𝑋 = 𝛼 2𝑇 2
𝜕𝑡 𝜕𝑥 2 𝜕𝑡 𝜕𝑥

Dividing both sides by XT:

1 𝜕𝑇 1 𝜕 2𝑋 1 𝜕𝑇 1 𝜕 2 𝑋
= 𝛼2 𝑜𝑟 = = −𝝀
𝑇 𝜕𝑡 𝑋 𝜕𝑥 2 𝛼 2 𝑇 𝜕𝑡 𝑋 𝜕𝑥 2

Since the LHS (function of t) and RHS (function of x) of the previous equation must be equal, the
only possibility is that both are equal to a constant 𝝀 (for convenience -𝝀). We then get two ODEs:

𝜕 2𝑋 𝜕𝑇
∴ + 𝜆𝑋 = 0 (1) ∴ + 𝛼 2 𝜆𝑋 = 0 (2)
𝜕𝑥 2 𝜕𝑡

𝑋(𝑥) = 𝐴 𝑠𝑖𝑛ξ𝜆𝑥 + 𝐵 𝑐𝑜𝑠ξ𝜆𝑥

𝑋′(𝑥) = 𝐴 ξ𝜆𝑐𝑜𝑠ξ𝜆𝑥 − 𝐵ξ𝜆 𝑠𝑖𝑛ξ𝜆𝑥
𝑋′′(𝑥) = −𝐴𝜆 𝑠𝑖𝑛ξ𝜆𝑥 − 𝐵 𝜆𝑐𝑜𝑠ξ𝜆𝑥
General solution:
𝑿(𝒙) = 𝑨 𝒔𝒊𝒏ξ𝝀𝒙 + 𝑩 𝒄𝒐𝒔ξ𝝀𝒙

To continue, we need to use the boundary condition:

𝑢(0,𝑡) = 0 𝑢(𝐿,𝑡) = 0
⇒ 𝑋(0)𝑇(𝑡) = 0
⇒ 𝑋(𝐿) = 𝐴 sinξ𝜆𝑥 = 0
We require: 𝑋(0) = 0
⇒ ξ𝜆𝑥 = 𝑛𝜋; 𝑛 = 1, 2, 3…
𝑋(0) = 𝐴 sinξ𝜆𝑥 + 𝐵 cosξ𝜆𝑥
⇒𝐵=0
𝒏𝟐 𝝅 𝟐
∴ 𝝀=
∴ X(x) = A sinξ𝝀𝒙 𝑳𝟐

Example
Page | 17
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Now we apply the principle of superposition to the initial PDE: if 𝒖𝟏 and 𝒖𝟐 are solutions to the PDE, so are
, 𝑪𝟏 𝒖𝟏 and 𝑪𝟐 𝒖𝟐 .

𝜕 𝐶1 𝑢 1 𝜕 2 𝐶1 𝑢 1 𝜕 𝑢1 𝜕 2 𝑢1 𝜕 𝑢1 𝜕 2 𝑢1
= 𝛼2 ⟺ 𝐶1 = 𝛼 2 𝐶1 ⟺ = 𝛼2
𝜕𝑡 𝜕𝑥 2 𝜕𝑡 𝜕𝑥 2 𝜕𝑡 𝜕𝑥 2

Similarly, 𝑪𝟏 𝒖𝟏 + 𝑪𝟐 𝒖𝟐 will satisfy the PDE:

𝑈 = 𝐶1 𝑢1 + 𝐶2 𝑢2 + ⋯ + 𝐶𝑛 𝑢𝑛

𝑛2 𝜋 2
Previously we had the eigenvalue 𝝀 = and the eigenfunction 𝑋(𝑥) = 𝐴 sinξ𝜆𝑥
𝐿2

⟹ 𝑋(𝑥) = 𝐴 sinξ𝜆𝑥, 𝑛 = 1,2,3 …

Linear algebra:

𝐴𝑋 = 𝜆𝑋
𝐴 = 𝑛 × 𝑛 matrix
𝑋 = 𝑛 × 1 column vector
𝜆 = scalar
𝐴𝑋 − 𝜆𝑋 = 0
(𝐴 − 𝜆𝐼) 𝑋 = 0

We have:

det(𝐴 − 𝜆𝐼) = 0
|𝐴 − 𝜆𝐼| = 0
𝝀 is called an eigenvalue
X is called an eigenvector

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

𝑑𝑇 𝑑𝑇
+ 𝛼 2 𝜆𝑇 = 0 ⟺ = −𝛼 2 𝜆𝑇
𝑑𝑡 𝑑𝑡

Example
Page | 18
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

The previous ODE can be solved using the principle of Separation of variables:

𝑑𝑇 𝑇
න = − 𝛼 2 𝜆 න 𝑑𝑡 ⟺ 𝑙𝑛𝑇 = −𝛼 2 𝜆𝑡 + 𝑙𝑛𝐵 ⟺ 𝑙𝑛 = −𝛼 2 𝜆𝑡
𝑇 𝐵
𝑇 𝑛2 𝜋 2 𝒏𝟐 𝝅𝟐
2 −𝜶𝟐 𝒕
= 𝑒 −𝛼 𝜆𝑡 ; 𝜆= 2 ; ∴ 𝑻 = 𝑩𝒆 𝑳𝟐
𝐵 𝐿

∞ 𝟐 𝟐
𝒏𝝅 𝒏 𝝅
−𝜶𝟐 𝟐 𝒕
𝒖𝒏 (𝒙, 𝒕) = ෍ 𝒃𝒏 𝐬𝐢𝐧 ቀ 𝒙ቁ 𝒆 𝑳
𝑳
𝒏=𝟏

The principle of
superposition

We know that:
𝑢(𝑥, 0) = 𝑓(𝑥)

∞
𝑛𝜋
∴ 𝑓(𝑥) = ෍ 𝑏𝑛 sin ቀ 𝑥ቁ
𝐿
𝑛=1

It can be observed that f(x) is a Fourier Series, in which case we know the formula for 𝒃𝒏 :

𝟐 𝑳 𝒏𝝅 𝑳
𝒃𝒏 = න 𝒇(𝒙)𝐬𝐢𝐧 ቀ 𝒙ቁ 𝒅𝒙 = 𝟐 න 𝒇(𝒙)𝐬𝐢𝐧(𝒏𝝅𝒙)𝒅𝒙
𝑳 𝟎 𝑳 𝟎

Example

Next, we will be looking at an application of the method of separation of variables to solve a PDE applicable
to a real-life scenario.

Page | 19
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

We have an insulated wire of length one unit, and the ends of the wire are held at 0℃ and let 𝛼 2 = 0.003.
The initial distribution is 𝑢(𝑥, 0) = 50𝑥(1 − 𝑥).

𝑇ℎ𝑒𝑟𝑒𝑓𝑜𝑟𝑒 𝑤𝑒 ℎ𝑎𝑣𝑒:

𝜕𝑇 𝜕 2𝑢
= 0.003 2
𝜕𝑡 𝜕𝑥

𝑢(0, 𝑡) = 𝑢(1, 𝑡) = 0; 𝐿 = 1;

𝑢(𝑥, 0) = 50𝑥(1 − 𝑥); 0 < 𝑥 < 1;

We had:
∞
𝑛𝜋
𝑓(𝑥) = ෍ 𝑏𝑛 sin ቀ 𝑥ቁ
𝐿
𝑛=1

This gives rise to a Fourier Series, which allows us to solve for 𝑏𝑛 :

2 𝐿 𝑛𝜋
𝑏𝑛 = න 𝑓(𝑥)sin ቀ 𝑥ቁ 𝑑𝑥; 𝑛 = 1,2,3 …
𝐿 0 𝐿

𝐿
𝑏𝑛 = 100 න 𝑥(1 − 𝑥)sin(𝑛𝜋𝑥)𝑑𝑥
0

𝐿
𝑏𝑛 = 100 න (𝑥 − 𝑥 2 )sin(𝑛𝜋𝑥)𝑑𝑥
0

Using LIATE:

𝑑𝑣 𝑑𝑢
න𝑢 𝑑𝑥 = 𝑢𝑣 − න 𝑣 𝑑𝑥
𝑑𝑥 𝑑𝑥

Let:

𝑑𝑢
𝑢 = 𝑥 − 𝑥2; = 1 − 2𝑥;
𝑑𝑥

𝑑𝑣 cos(𝑛𝜋𝑥)
= sin(𝑛𝜋𝑥); 𝑣=−
𝑑𝑥 𝑛𝜋

Example
Page | 20
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

cos(𝑛𝜋𝑥) 𝐿 𝐿
cos(𝑛𝜋𝑥)
∴ 𝑏𝑛 = 100 ቊ൤(𝑥 − 𝑥 2 ) ൨ + න (1 − 2𝑥)𝑑𝑥 ቋ
𝑛𝜋 0 0 𝑛𝜋

100 𝐿
𝑏𝑛 = න cos(𝑛𝜋𝑥) (1 − 2𝑥)𝑑𝑥
𝑛𝜋 0

𝑑𝑢
𝑢 = 1 − 2𝑥; = −2;
𝑑𝑥

𝑑𝑣 sin(𝑛𝜋𝑥)
= cos(𝑛𝜋𝑥); 𝑣=
𝑑𝑥 𝑛𝜋

100 sin(𝑛𝜋𝑥) 𝐿 𝐿
sin(𝑛𝜋𝑥)
∴ 𝑏𝑛 = ቊ൤(1 − 2𝑥) ൨ + 2න 𝑑𝑥 ቋ
𝑛𝜋 𝑛𝜋 0 0 𝑛𝜋

200 𝐿
𝑏𝑛 = න sin(𝑛𝜋𝑥) 𝑑𝑥
𝑛2 𝜋 2 0

−200 cos(𝑛𝜋𝑥) 𝐿
𝑏𝑛 = ൤ ൨
𝑛2 𝜋 2 𝑛𝜋 0

−200
𝑏𝑛 = [cos(𝑛𝜋𝐿) − cos(0)]
𝑛3 𝜋 3

𝑒𝑣𝑒𝑛, 𝑡ℎ𝑒𝑛: 𝑏𝑛 = 0
If n is ቊ 400
𝑜𝑑𝑑, 𝑡ℎ𝑒𝑛: 𝑏𝑛 = 𝑛3 𝜋3

∞
𝟒𝟎𝟎 𝟐 𝟐
∴ 𝒖(𝒙, 𝒕) = ෍ 𝟑 𝟑
𝐬𝐢𝐧(𝒏𝝅𝒙) 𝒆−𝒏 𝝅 𝟎.𝟎𝟎𝟑𝒕
𝒏=𝟏
𝒏 𝝅
𝒏 𝒐𝒅𝒅

Example

For further practice, try to plot the function in MATLAB and determine whether the maximum of u(x,t) will
occur at x=0.5.

Page | 21
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Worksheet Questions – Tutorial 3

1. Find solutions 𝑢(𝑥, 𝑦) of the following equations by separating variables:

i. 𝑢𝑥 + 𝑢𝑦 = 0
ii. 𝑢𝑥 − 𝑢𝑦 = 0
iii. 𝑦 2 𝑢𝑥 − 𝑥 2 𝑢𝑦 = 0
iv. 𝑢𝑥 + 𝑢𝑦 = (𝑥 + 𝑦)𝑢
v. 𝑢𝑥𝑥 + 𝑢𝑦𝑦 = 0
vi. 𝑢𝑥𝑦 − 𝑢 = 0
vii. 𝑢𝑥𝑥 − 𝑢𝑦𝑦 = 0
viii. 𝑥𝑢𝑥𝑦 + 2𝑦𝑢 = 0

2. Using the method of Separation of Variables and Fourier series find the solution to the Heat
Equation (Boundary Value Problem):

𝜕𝑢 𝜕 2𝑢
= 𝑘 2,
𝜕𝑡 𝜕𝑥

On 𝛀 = {𝟎 < 𝐱 < 𝐥, 𝐭 > 𝟎}

Given the flowing Dirichlet conditions on 𝜕Ω

𝑢(𝑥, 0) = 𝑓(𝑥), 0 < 𝑥 < 𝑙

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

3. Using the method of Separation of Variables and Fourier series find the solution to the Laplace
equation (Boundary Value Problem):

𝜕 2𝑇 𝜕 2𝑇
+ =0
𝜕𝑥 2 𝜕𝑦 2

For the temperature distribution along a bar as shown in the following figure:

Page | 22
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Subject to:

𝜕𝑇(𝑥,∞)
= 0 (right)
𝜕𝑥

𝜕𝑇(0,𝑦)
= 0 (bottom)
𝜕𝑦

𝜕𝑇(10,𝑦)
= 0 (top)
𝜕𝑦

𝜕𝑇(𝑥,0)
= 100 (left)
𝜕𝑥

4. Using the method of Separation of Variables and Fourier series find the solution to the one-
dimensional wave equation:

𝜕 2𝑢 2
𝜕 2𝑢
= 𝑐
𝜕𝑡 2 𝜕𝑥 2

Where u(x,t) is the deflection of the string. The string is fixed at both ends at x=0 and x=L, thus
the following boundary conditions apply.

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

𝑢(𝑥, 0) = 𝑓(𝑥),

and

𝜕𝑢(𝑥, 0)
= 𝑔(𝑥)
𝜕𝑡

Page | 23
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

5. Using the method of Separation of Variables and Fourier series find the solution to the Heat
equation (Boundary Value Problem):

𝜕𝑢 𝜕 2𝑢
= 𝑘 2,
𝜕𝑡 𝜕𝑥

On 𝛀 = {𝟎 < 𝐱 < 𝐥, 𝐭 > 𝟎}

Given the flowing Neumann and Dirichlet conditions on 𝜕Ω

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

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

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

6. Using the method of Separation of Variables series find the solution to the Boundary Value
Problem:

𝜕𝑢 𝜕 2𝑢
=𝑘 2 −𝑢
𝜕𝑡 𝜕𝑥

𝑢(𝑥, 0) = 𝑓(𝑥),

𝑢(0, 𝑡) = 0

𝜕𝑢(𝐿, 𝑡)
= −𝑢(𝐿, 𝑡)
𝜕𝑥

7. Find solutions 𝑢(𝑥, 𝑦) of the following equations by separating variables:

i. 𝑢𝑡 = 𝑘𝑢𝑥𝑥 0 < 𝑥 < 𝑙, 𝑡 > 0

𝑢(𝑥, 0) = 𝑓(𝑥) 0 < 𝑥 < 𝑙
𝑢(0, 𝑡) = 𝑢(𝑙, 𝑡) = 0 𝑡 > 0

Page | 24
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

𝜕2 𝑢 𝜕2𝑢 𝜕2𝑢
ii. ∇2 𝑢 = 0 𝑤ℎ𝑒𝑟𝑒 ∇2 𝑢 ≡ 𝜕𝑥 2 + 𝜕𝑦 2 + 𝜕𝑧 2 = 0

With 0 < 𝑥 < 𝐿, 0 < 𝑦 < 𝐿, 0 < 𝑧 < 𝐿

𝑎𝑛𝑑

𝑢 = 0 𝑜𝑛 𝑥 = 0, 𝐿
𝑢 = 0 𝑜𝑛 𝑧 = 0
𝑢 = 1 𝑜𝑛 𝑧 = 𝐿

In this case assume there is no dependency on y so that the separation solution to be

assumed is:

𝑢(𝑥, 𝑧) = 𝑋(𝑥)𝑍(𝑧) 𝑎𝑛𝑑 𝑡ℎ𝑎𝑡

𝜕 2𝑢 𝜕 2𝑢
+ =0
𝜕𝑥 2 𝜕𝑧 2

iii. 𝑢𝑡𝑡 = 𝑐 2 𝑢𝑥𝑥

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

Page | 25
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Numerical Solutions to PDEs

In most cases, the PDEs we are trying to solve, with the given boundary conditions, are unsolvable in closed
form. In this case, the best approach is to come up with a numerical approximation. We will first look at some
examples of Taylor Series to obtain some equations for the derivatives that we will be using.
Let us consider the following formulations of Taylor Series. They can be visualized in Figure 9.

ℎ2 ℎ3
𝑓(𝑥 + ℎ) = 𝑓(𝑥) + ℎ𝑓 ′ (𝑥) + 𝑓 ′′ (𝑥) + 𝑓 ′′ ′(𝑥)+ …
2! 3!

ℎ2 ℎ3
𝑓(𝑥 − ℎ) = 𝑓(𝑥) − ℎ𝑓 ′ (𝑥) + 𝑓 ′′ (𝑥) − 𝑓 ′′ ′(𝑥)+ …
2! 3!

Equation 14

By subtracting the two Taylor series, we can

determine an equation for the finite difference
approximation of the 1st derivative of a function.

∴ 𝑓(𝑥 + ℎ) − 𝑓(𝑥 − ℎ) = 2ℎ𝑓 ′ (𝑥) + 𝛰(ℎ3 )

𝒇(𝒙 + 𝒉) − 𝒇(𝒙 − 𝒉)
∴ = 𝒇′ (𝒙) + 𝜪൫𝒉𝟐 ൯
𝟐𝒉

the error term

Figure 8 – Taylor Series Equation 15

Now, we can perform a test to see if the formula we determined yields correct results.

𝐹𝑜𝑟 𝑒𝑥𝑎𝑚𝑝𝑙𝑒, 𝑓(𝑥) = 𝑥 2

𝑓 ′ (𝑥) = 2𝑥; @𝑥 = 2, 𝑓 ′ (2) = 4

32 − 12
𝑥 + ℎ → 𝑓(3) ↗ ↖ 𝑓(1) ⟵ 𝑥 − ℎ
2
Figure 9 –
9−1
= 4 ⟹ 𝑛𝑜 𝑒𝑟𝑟𝑜𝑟 quadratic
2 function

Example

Page | 26
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Why are we not surprised?

𝑓(𝑥) = 𝑥 2 𝑓 ′ (𝑥) = 2𝑥 𝑓 ′′ (𝑥) = 2 𝑓 ′′′ (𝑥) = 0

𝑓(𝑥 + ℎ) − 𝑓(𝑥 − ℎ)
= 𝑓 ′ (𝑥) + 𝛰(ℎ2 )
2ℎ

It’s because the 3rd derivative of the function, and therefore the error term, is equal to zero.
However, for 𝑓(𝑥) = 𝑥 3 then the result would not be in complete agreement as 𝑓 ′′′ (𝑥) = 6.

Sometimes we also use:

𝑓(𝑥 + ℎ) − 𝑓(𝑥)
≈ 𝑓 ′ (𝑥) Figure 10 –
2ℎ slope of a
function

Example

If we sum the two Taylor series instead, we can obtain a finite difference approximation formula for the 2nd
derivative of a function.

2 ′′ (𝑥)
ℎ4 (𝑖𝑣)
𝑓(𝑥 − ℎ) − 2𝑓(𝑥) + 𝑓(𝑥 + ℎ) = 2𝑓(𝑥) + ℎ 𝑓 + 2 𝑓 (𝑥) − 2𝑓(𝑥)
4!
𝑓(𝑥 − ℎ) − 2𝑓(𝑥) + 𝑓(𝑥 + ℎ)
∴ = 𝑓 ′′ (𝑥) + 𝛰(ℎ2 )
ℎ2

𝒇(𝒙 − 𝒉) − 𝟐𝒇(𝒙) + 𝒇(𝒙 + 𝒉)

∴ 𝒇′′ (𝒙) = + 𝜪(𝒉𝟐 )
𝒉𝟐

Equation 16

Page | 27
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Now let us consider the following example to test our expression.

𝑓(𝑥) = 𝑥 3

𝑓 ′ (𝑥) = 3𝑥 2 ; 𝑓 ′′ (𝑥) = 6𝑥;

@𝑥 = 1, 𝑓 ′′ (1) = 6;

@𝑥 = 2, 𝑓 ′′ (2) = 12;

13 − 2(23 ) + 33
𝑓 ′′ (2) = =
12

1 − 16 + 27
= = 12 𝑄. 𝐸. 𝐷. Figure 11 – Cubic equation
1
Example

To test your knowledge, find an expression for 𝑓 ′′′ (𝑥) as a finite difference form and compare your
numerical result with the function 𝑓(𝑥) = 𝑥 4 .

Solution Technique
For numerical solution of elliptic PDEs, the PDE is
transformed into an algebraic difference equation.
A grid system is used to divide the region of
interest. Since the PDE is satisfied at each point in
the area, it must be satisfied at each point of the
grid. A finite difference approximation is obtained
at each grid point.

𝒖𝒊+𝟏,𝒋 + 𝒖𝒊−𝟏,𝒋 + 𝒖𝒊,𝒋+𝟏 + 𝒖𝒊,𝒋−𝟏 − 𝟒𝒖𝒊,𝒋 = 𝟎

Figure 12 – Grid system

Page | 28
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Here you can see a worked example of Laplace’s Equation using the finite difference approximation on a
given grid system.

u(x,y)=400
xy

Figure 13 – Mesh system with boundary conditions

𝛻 2 u=0

𝜕2𝑢 𝜕2𝑢
+ = 0 𝑖𝑛 Ω = {(x, y)ȁ0 < x < 0.5, 0 < y < 0.5) 𝑜𝑛 𝜕Ω.
𝜕𝑥 2 𝜕𝑦 2

We know that u(x,y)=400xy is the exact solution and previously we have demonstrated that:

𝑓(𝑥 − ℎ) − 2𝑓(𝑥) + 𝑓(𝑥 + ℎ)

𝑓 ′′ (𝑥) =
ℎ2

Uniform mesh:

Δ𝑥 = Δ𝑦 = ℎ = 0.125

Dirichlet Conditions:

u(x,0)=0 u(0,y)=0

u(x, 0.5)=200x u(0.5, y)=200y

Example

Page | 29
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

, +

,
, + ,

Figure 14 – Determining the computational molecule

Each point on the grid is expressed in terms of its coordinates on the x and y axes:

𝑗
𝑢𝑖,𝑗 (𝑖Δ𝑥, 𝑗Δ𝑦) 𝑜𝑟 𝑢𝑖 (𝑖Δ𝑥, 𝑗Δ𝑦)

We know that

𝜕2 𝑢 𝜕2 𝑢
𝜕𝑥 2
+ 𝜕𝑦2 = 0 𝑖𝑛 Ω and 𝜕Ω.

By using the finite difference formula for the 2nd derivative, we obtain:

𝑢𝑖−1,𝑗 − 2𝑢𝑖,𝑗 + 𝑢𝑖+1,𝑗 𝑢𝑖,𝑗−1 − 2𝑢𝑖,𝑗 + 𝑢𝑖,𝑗+1

+ =0
ℎ2 ℎ2

𝑢𝑖−1,𝑗 − 4𝑢𝑖,𝑗 + 𝑢𝑖+1,𝑗 + 𝑢𝑖,𝑗−1 + 𝑢𝑖,𝑗+1 = 0

1
𝑢𝑖,𝑗 = [𝑢𝑖−1,𝑗 + 𝑢𝑖+1,𝑗 + 𝑢𝑖,𝑗−1 + 𝑢𝑖,𝑗+1 ]
4

This is the general formula for any given point 𝒖𝒊,𝒋 on the grid, known as the computational molecule.
Now, we have to calculate all nine points following this expression.

Example

Page | 30
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Figure 15 – Mesh system with boundary conditions and mesh points

1 1
𝑢1 = [𝑢2 + 0 + 𝑢4 + 25] 𝑢6 = [𝑢3 + 𝑢5 + 𝑢9 + 200 × 0.25]
4 4
1 1
𝑢1 = [𝑢2 + 𝑢4 + 25] 𝑢6 = [𝑢3 + 𝑢5 + 𝑢9 + 50]
4 4

1 1
𝑢2 = [𝑢1 + 𝑢3 + 𝑢5 + 200 × 0.25] 𝑢7 = [𝑢8 + 𝑢4 + 0 + 0]
4 4
1 1
𝑢2 = [𝑢1 + 𝑢3 + 𝑢5 + 50] 𝑢7 = [𝑢8 + 𝑢4 ]
4 4

1 1
𝑢3 = [𝑢2 + 𝑢6 + 200 × 0.375 + 200 × 0.375] 𝑢8 = [𝑢5 + 𝑢7 + +0 + 𝑢9 ]
4 4
1 1
𝑢3 = [𝑢2 + 𝑢6 + 150] 𝑢8 = [𝑢5 + 𝑢7 + 𝑢9 ]
4 4

1 1 200
𝑢4 = [𝑢5 + 𝑢1 + 𝑢7 + 0] 𝑢9 = ൤𝑢6 + 𝑢8 + 0 + ൨
4 4 8
1 1
𝑢4 = [𝑢5 + 𝑢1 + 𝑢7 ] 𝑢9 = [𝑢6 + 𝑢8 + 25]
4 4

1
𝑢5 = [𝑢2 + 𝑢6 + 𝑢4 + 𝑢8 ]
4
Example

Page | 31
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

This can be solved using a matrix. First the equations need to be rearranged:

4𝑢1 −𝑢2 − 𝑢4 = 25

4𝑢2 − 𝑢1 − 𝑢3 − 𝑢5 = 50

4𝑢3 − 𝑢2 − 𝑢6 = 150

4𝑢4 − 𝑢5 − 𝑢1 − 𝑢7 = 0

4𝑢5 − 𝑢2 − 𝑢6 − 𝑢4 − 𝑢8 = 0

4𝑢6 − 𝑢3 − 𝑢5 − 𝑢9 = 50

4𝑢7 − 𝑢8 − 𝑢4 = 0

4𝑢8 − 𝑢5 − 𝑢7 − 𝑢9 = 0

4𝑢9 − 𝑢6 − 𝑢8 = 25

4 −1 0 −1 0 0 0 0 0 𝑢1 25
−1 4 −1 0 −1 0 0 0 0 𝑢 2 50
‫ ۇ‬0 −1 4 0 0 −1 0 0 0 ‫ۇ ۊ 𝑢 ۇ ۊ‬ ‫ۊ‬
3 150
‫ۈ‬ ‫ۈ ۋ 𝑢ۈۋ‬ ‫ۋ‬
‫ۈ‬−1 0 0 4 −1 0 −1 0 0 ‫ ۈ ۋ‬4 ‫ ۈ ۋ‬0 ‫ۋ‬
‫ ۈ‬0 −1 0 −1 4 −10 0 −1 0 ‫𝑢 ۈ ۋ‬5 ‫ ۈ = ۋ‬0 ‫ۋ‬
‫ ۈ‬0 0 −1 0 −1 4 0 0 −1‫𝑢 ۈ ۋ‬6 ‫ ۈ ۋ‬50 ‫ۋ‬
‫ ۈ‬0 0 0 −1 0 0 4 −1 0 ‫𝑢 ۈ ۋ‬7 ‫ ۈ ۋ‬0 ‫ۋ‬
‫ ۈ‬0 0 0 0 −1 0 −1 4 −1‫𝑢 ۈ ۋ‬8 ‫ ۈ ۋ‬0 ‫ۋ‬
0 0 0 0 0 −1 0 −1 4 𝑢9 25
‫ۉ‬ ‫ۉ ی ۉی‬ ‫ی‬

Solution:

75 75
𝑢1 = 𝑢5 = 25 𝑢9 =
4 4

75 75
𝑢2 = 𝑢6 = Verify:
2 2

225 25 𝑢9 = 400𝑥𝑦 =
𝑢3 = 𝑢7 =
4 4
= 400×0.125× 0.375 =
25 25 75
𝑢4 = 𝑢8 = = 18.75 = 4 ⟹ exact solution!
2 2

Example

Page | 32
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Worksheet Questions – Tutorial 4

1. Defining

𝛀 ={(𝐱,𝐲)|𝟎<𝐱<𝟏,𝟎<𝐲<𝟏}

Find the numerical solution to the 2-dimensional Poisson Equation

𝜕 2𝑢 𝜕 2𝑢
+ = 2(𝑦 − 𝑥)
𝜕𝑥 2 𝜕𝑦 2

Given the flowing Dirichlet conditions on 𝜕Ω

𝑢(0, 𝑦) = 0, 0 < 𝑦 < 1

𝑢(1, 𝑦) = 𝑦 − 𝑦 2 , 0 < 𝑦 < 1

𝑢(𝑥, 0) = 0, 0 < 𝑥 < 1

𝑢(𝑥, 1) = 𝑥 2 − 𝑥, 0 < 𝑥 < 1

Compare your numerical solution with the exact solution given by:

𝑢(𝑥, 𝑦) = 𝑥 2 𝑦 − 𝑦 2 𝑥,

Use a uniform mesh spacing of 0.1.

2. Defining

Ω = {(x, y)ȁ0 < x < 1,0 < y < 1}

Find the numerical solution to the 2-dimensional Laplace Equation

𝜕 2𝑢 𝜕 2𝑢
+ =0
𝜕𝑥 2 𝜕𝑦 2

Page | 33
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Given the flowing Dirichlet conditions on 𝜕Ω

𝑢(0, 𝑦) = 𝑦, 0 < 𝑦 < 1

𝑢(1, 𝑦) = 301𝑦 + 1, 0 < 𝑦 < 1

𝑢(𝑥, 0) = 𝑥, 0 < 𝑥 < 1

𝑢(𝑥, 1) = 301𝑥 + 1, 0 < 𝑥 < 1

Compare your numerical solution with the exact solution given by:

𝑢(𝑥, 𝑦) = 300𝑥𝑦 + 𝑥 + 𝑦

Use a uniform mesh spacing of 0.1.

3. Defining

Ω = {(x, y)ȁ0 < x < 1.5,0 < y < 1.5}

Find the numerical solution to the 2-dimensional PDE:

𝜕 2𝑢 𝜕 2𝑢
− = (𝑦 2 − 𝑥 2 )sinh(𝑥𝑦)
𝜕𝑥 2 𝜕𝑦 2

Given the flowing Dirichlet conditions on 𝜕Ω

𝑢(0, 𝑦) = 0, 0 < 𝑦 < 1.5

𝑢(1.5, 𝑦) = sinh (1.5𝑦), 0 < 𝑦 < 1.5

𝑢(𝑥, 0) = 0, 0 < 𝑥 < 1.5

𝑢(𝑥, 1.5) = sinh (1.5𝑥), 0 < 𝑥 < 1.5

Page | 34
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Compare your numerical solution with the exact solution given by:

𝑢(𝑥, 𝑦) = sinh(𝑥𝑦)

Use a uniform mesh spacing of 0.1.

4. Defining

Ω = {(x, y)ȁ0 < x < 2.0,0 < y < 2.0}

Find the numerical solution to the 2-dimensional PDE:

𝜕 2 𝑢 𝜕 2 𝑢 𝜕𝑢 𝜕𝑢
+ + + = 𝑒 𝑥 (𝑐𝑜𝑠𝑦 − 𝑠𝑖𝑛𝑦)
𝜕𝑥 2 𝜕𝑦 2 𝜕𝑥 𝜕𝑦

Given the flowing Dirichlet conditions on 𝜕Ω

𝑢(0, 𝑦) = 𝑐𝑜𝑠𝑦, 0 < 𝑦 < 2

𝑢(2, 𝑦) = 𝑒 2 cosy, 0 < 𝑦 < 2.0

𝑢(𝑥, 0) = 𝑒 𝑥 , 0 < 𝑥 < 2.0

𝑢(𝑥, 2) = 𝑒 𝑥 cos(2) , 0 < 𝑥 < 2.0

Compare your numerical solution with the exact solution given by:

𝑢(𝑥, 𝑦) = 𝑒 𝑥 cos(𝑦)

Note that all trigonometric calculations must be in radians.

Use a uniform mesh spacing of 0.2.

Page | 35
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

5. Numerical solution to PDEs

𝜕 2𝑢 𝜕 2𝑢
+ = 4 0 < 𝑥 < 1; 0<𝑦<2
𝜕𝑥 2 𝜕𝑦 2

𝑢(𝑥, 0) = 𝑥 2 , 𝑢(𝑥, 2) = (𝑥 − 2)2 0 ≤ 𝑥 ≤ 1

𝑢(0, 𝑦) = 𝑦 2 , 𝑢(1, 𝑦) = (𝑦 − 1)2 0 ≤ 𝑦 ≤ 2

Use
ℎ = 0.25, 𝑘 = 0.5 𝑎𝑛𝑑 𝑐𝑜𝑚𝑝𝑎𝑟𝑒 𝑦𝑜𝑢𝑟 𝑟𝑒𝑠𝑢𝑙𝑡𝑠 𝑡𝑜 𝑡ℎ𝑒 𝑒𝑥𝑎𝑐𝑡 𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛
𝑢(𝑥, 𝑦) = (𝑥 − 𝑦)2 ,

6. Numerical solution to PDEs

𝜕 2𝑢 𝜕 2𝑢
+ = (𝑥 2 + 𝑦 2 )𝑒 𝑥𝑦 0 < 𝑥 < 2; 0<𝑦<1
𝜕𝑥 2 𝜕𝑦 2

𝑢(0, 𝑦) = 1, 𝑢(2, 𝑦) = 𝑒 2𝑦 0 ≤ 𝑦 ≤ 1
𝑢(𝑥, 0) = 1, 𝑢(𝑥, 1) = 𝑒 𝑥 0 ≤ 𝑥 ≤ 2

Use
ℎ = 0.2, 𝑘 = 0.1 𝑎𝑛𝑑 𝑐𝑜𝑚𝑝𝑎𝑟𝑒 𝑦𝑜𝑢𝑟 𝑟𝑒𝑠𝑢𝑙𝑡𝑠 𝑡𝑜 𝑡ℎ𝑒 𝑒𝑥𝑎𝑐𝑡 𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛
𝑢(𝑥, 𝑦) = 𝑒 𝑥𝑦 ,

7. Numerical solution to PDEs

𝜕 2𝑢 𝜕 2𝑢 𝑥 𝑦
+ = + 1 < 𝑥 < 2; 1<𝑦<1
𝜕𝑥 2 𝜕𝑦 2 𝑦 𝑥

𝑢(𝑥, 1) = 𝑥𝑙𝑛(𝑥), 𝑢(𝑥, 2) = 𝑥𝑙𝑛(4𝑥 2 ) 1 ≤ 𝑥 ≤ 2

𝑢(1, 𝑦) = 𝑦𝑙𝑛(𝑦), 𝑢(2, 𝑦) = 2𝑦𝑙𝑛(2𝑦) 1 ≤ 𝑦 ≤ 2

Use
ℎ = 𝑘 = 0.1, 𝑎𝑛𝑑 𝑐𝑜𝑚𝑝𝑎𝑟𝑒 𝑦𝑜𝑢𝑟 𝑟𝑒𝑠𝑢𝑙𝑡𝑠 𝑡𝑜 𝑡ℎ𝑒 𝑒𝑥𝑎𝑐𝑡 𝑠𝑜𝑙𝑢𝑡𝑖𝑜𝑛
𝑢(𝑥, 𝑦) = 𝑥𝑦𝑙𝑛(𝑥𝑦),

Page | 36
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Worked Examples of PDEs

ENGF0004 Worked Examples.

Question One: Partial Differential Equations.

a) Use the Separation of Variables technique to find the solution of the following boundary value problem

𝜕𝑢 𝜕2𝑢
=𝑘 2
𝜕𝑡 𝜕𝑥

with boundary conditions

𝑢(𝑥, 0) = 𝑓(𝑥) 0 ≤ 𝑥 ≤ 𝑙

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

Take 𝑓(𝑥) = 𝑥 2 on 0 ≤ 𝑥 ≤ 𝑙 so that you can determine the Fourier coefficients.

i) Show using the method of separation of variables that the following 2nd order ordinary
differential equation is obtained
−𝑋 ′′ (𝑥) = 𝜇𝑋(𝑥)
As well as
𝑇 ′ (𝑡) = −𝜇𝑘𝑇(𝑡)

where 𝑢(𝑥, 𝑡) = 𝑋(𝑥)𝑇(𝑡) [10]

ii) By finding the solution using the conditions given, show that

∞
𝑛𝜋𝑥
𝑓(𝑥) = ෍ 𝑐𝑛 sin ቀ ቁ
𝑙
𝑛=1

where the 𝑐𝑛 are constants. [10]

Page | 37
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

iii) Complete the solution for 𝑢(𝑥, 𝑡) by determining the 𝑐𝑛 taking 𝑓(𝑥) = 𝑥 2 [20]

b) Consider the following Poisson equation:

𝜕2𝑢 𝜕2𝑢
+ = 𝑥𝑒 𝑦 for 𝜕𝜔 = {(𝑥, 𝑦)ȁ0 ≤ 𝑥 ≤ 2 and 0 ≤ 𝑦 ≤ 1}
𝜕𝑥 2 𝜕𝑦 2

with the boundary conditions

𝑢(0, 𝑦) = 0, 𝑢(2, 𝑦) = 2𝑒 𝑦 for 0 ≤ 𝑦 ≤ 1

𝑢(𝑥, 0) = 𝑥, 𝑢(𝑥, 1) = 𝑥𝑒 for 0 ≤ 𝑥 ≤ 2 where 𝒆 means 𝑒 1

Use a mesh spacing of ℎ = 𝛿𝑥 = 0.4 and 𝑘 = 𝛿𝑦 = 0.2

𝜕2 𝑢 𝜕2 𝑢
i) Using a Taylor series, derive the central difference approximations for 𝜕𝑥 2 and 𝜕𝑦2 at the point (𝑖ℎ, 𝑗𝑘)

and quote the order of the error terms. [5]

ii) Give the difference equation that must be satisfied by 𝑢(𝑖ℎ, 𝑗𝑘) at each point (𝑖ℎ, 𝑗𝑘) after discretizing
𝜕𝜔. [5]

iii) Obtain the 16 difference equations that must be solved in order to obtain the numerical solution to
the given boundary value problem. [10]

iv) Write MATLAB code that solves this boundary value problem on a mesh with ℎ = 𝛿𝑥 = 0.4 and 𝑘 =
𝛿𝑦 = 0.2. You may use an iterative scheme of your choice such as the Gauss-Seidel technique with
the stopping criteria:

(𝑛) (𝑛−1)
|𝑢𝑖,𝑗 − 𝑢𝑖,𝑗 | < 10−5

where the superscript denotes the iteration stage.

Alternatively, you may use Cramer’s rule, or you may choose to invert the matrix obtained from the
algebraic formulation of this problem (i.e., when you write the PDE as a system of equations):

𝐴 𝑤𝑖 = 𝑏𝑖

where 𝐴 is a 16 × 16 matrix, the 𝑤𝑖 is a 16 × 1 column vector of unknown interior points and the 𝑏𝑖 is
a 16 × 1 column vector of known value. [35]

Page | 38
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

v) Within your MATLAB code, compare your numerical solutions (the 16 interior points) with the exact
solution given by 𝑢(𝑥, 𝑦) = 𝑥𝑒 𝑦 . [5]

SOLUTION:
Question Two:

a)

𝜕𝑢 𝜕2𝑢
=𝑘 2
𝜕𝑡 𝜕𝑥

𝑢(𝑥, 0) = 𝑓(𝑥) 0 ≤ 𝑥 ≤ 𝑙
𝑢(0, 𝑡) = 𝑢(𝑙, 𝑡) = 0 𝑡 > 0

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

Substituting 𝑢(𝑥, 𝑡) = 𝑋(𝑥)𝑇(𝑡) in the given PDE:

𝑋(𝑥)𝑇 ′ (𝑡) = 𝑘𝑋 ′′ (𝑥)𝑇(𝑡)
𝑇′(𝑡) 𝑋′′(𝑥)
=
𝑘𝑇(𝑡) 𝑋(𝑥)
[5]

𝑇 ′ (𝑡) 𝑋 ′′ (𝑥)
− =− =𝜇
𝑘𝑇(𝑡) 𝑋(𝑥)
where 𝜇 is an arbitrary constant.

Therefore: −𝑿′′ (𝒙) = 𝝁𝑿(𝒙) and 𝑻′ (𝒕) = −𝝁𝒌𝑻(𝒕) [5]

(ii) For negative constant (𝜆 = −𝜇 < 𝟎)

−𝑋 ′′ (𝑥) = 𝜇𝑋(𝑥)

𝑋 ′′ (𝑥) + 𝜇𝑋(𝑥) = 0

The solution to the auxiliary equation is thus: 𝑚 = ±𝑖ξ𝜇

The general form of the solution for 𝑋(𝑥) is: 𝑿(𝒙) = 𝑨𝐜𝐨𝐬(ξ𝝁𝒙) + 𝑩𝐬𝐢𝐧(ξ𝝁𝒙) [1]

Page | 39
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Applying 𝑢(0, 𝑡) = 0:
𝑋(0)𝑇(𝑡) = 0 𝐴𝑇(𝑡) = 0

𝑇(𝑡) = 0 or 𝐴 = 0. 𝑇(𝑡) = 0 will lead to a trivial solution of 𝑢(𝑥, 𝑡) = 0. So: 𝐴 = 0

𝑿(𝒙) = 𝑩𝐬𝐢𝐧൫√𝝁𝒙൯

[2]
Applying 𝑢(𝑙, 𝑡) = 0:
𝑋(𝑙)𝑇(𝑡) = 0 𝐵sin(ξ𝝁𝑙)𝑇(𝑡) = 0
𝑇(𝑡) = 0 or 𝐵 = 0 or sin(ξ𝜇𝑙) = 0
𝑇(𝑡) = 0 and 𝐵 = 0 will lead to a trivial solution of 𝑢(𝑥, 𝑡) = 0. So:

sin൫√𝜇𝑙൯ = 0

√𝜇𝑙 = 𝑛𝜋 𝑛 = 1,2,3 …

𝑛𝜋 𝑛𝜋 2
√𝜇𝑛 = 𝑙
−𝜇𝑛 = − ቀ 𝑙 ቁ

𝑛𝜋𝑥
Therefore: 𝑋(𝑥) = 𝐵𝑛 sin ቀ ቁ 𝑛 = 1,2,3, …
𝑙

[3]
Now for the second ODE in 𝑇(𝑡):
𝑇 ′ (𝑡) = −𝜇𝑛 𝑘𝑇(𝑡)
𝑇 ′ (𝑡)
= − 𝜇𝑛 𝑘
𝑇(𝑡)
𝑇 ′ (𝑡)
න 𝑑𝑡 = න −𝜇𝑛 𝑘 𝑑𝑡
𝑇(𝑡)
ln 𝑇(𝑡) = −𝜇𝑛 𝑘𝑡 + 𝑐

𝑇(𝑡) = 𝑒 (−𝜇𝑛𝑘𝑡 + 𝑐)

𝑛𝜋 2
−𝜇𝑛 = − ቀ ቁ
𝑙

Page | 40
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

nπ 2
T(t) = Dn exp ൤− ቀ l ቁ kt൨ where Dn = ec
[2]
nπ 2 nπx
X(x)T(t) = cn exp ൤− ቀ ቁ kt൨ sin ቀ ቁ
l l

where 𝒄𝒏 = 𝑫𝒏 𝑩𝒏

By the principle of superposition:

∞
𝑛𝜋 2 𝑛𝜋𝑥
𝑢(𝑥, 𝑡) = ෍ 𝑐𝑛 𝑒𝑥𝑝 ൤− ቀ ቁ 𝑘𝑡൨ sin ቀ ቁ
𝑙 𝑙
𝑛=1

[1]
Applying 𝑢(𝑥, 0) = 𝑓(𝑥)
∞
𝒏𝝅𝒙
𝒇(𝒙) = ෍ 𝒄𝒏 𝐬𝐢𝐧 ቀ ቁ
𝒍
𝒏=𝟏

[1]
(iii) 𝑐𝑛 values follow from the Fourier sine series for 𝑓(𝑥).
2 𝑙 𝑛𝜋𝑥
𝑐𝑛 = න 𝑓(𝑥) sin ቀ ቁ 𝑑𝑥
𝑙 0 𝑙

2 𝑙 𝑛𝜋𝑥
𝑐𝑛 = න 𝑥2 sin ቀ ቁ 𝑑𝑥
𝑙 0 𝑙

Using Integration by Parts (or from Standard Tables of Integrals):

2𝑥 𝑎2 𝑥 2 − 2
න 𝑥 2 sin(𝑎𝑥) 𝑑𝑥 = sin(𝑎𝑥) − cos(𝑎𝑥)
𝑎2 𝑎3
[5]
𝑙
𝑛𝜋 2 2
2 𝑙
𝑛𝜋𝑥 2 2𝑥 𝑛𝜋𝑥 ቀ
ቁ 𝑥 −2 𝑛𝜋𝑥
𝑐𝑛 = න 𝑥2 sin ቀ ቁ 𝑑𝑥 = [ sin ቀ ቁ −
𝑙 cos ቀ ቁ]
2 3
𝑙 0 𝑙 𝑙 𝑛𝜋 𝑙 𝑛𝜋 𝑙
ቀ ቁ ቀ ቁ
𝑙 𝑙 0

Page | 41
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

2 2𝑙 3 𝑙 3 (𝑛2 𝜋 2 − 2) 2𝑙 3
𝑐𝑛 = [ 2 2 sin(𝑛𝜋) − cos(𝑛𝜋) − 0 − 3 3 ]
𝑙 𝑛 𝜋 𝑛3 𝜋 3 𝑛 𝜋

[10]
sin(𝑛𝜋) = 0 cos(𝑛𝜋) = (−1)𝑛 𝑛 = 1,2,3, … [1]

3 2 2
2 𝑛 𝑙 (𝑛 𝜋 − 2) 2𝑙 3
𝑐𝑛 = [−(−1) − 3 3]
𝑙 𝑛3 𝜋 3 𝑛 𝜋

𝟐𝒍𝟐
𝒄𝒏 = [(𝟐 − 𝒏𝟐 𝝅𝟐 )(−𝟏)𝒏 − 𝟐]
𝒏𝟑 𝝅 𝟑
[2]
Therefore:
∞
𝒏𝝅 𝟐 𝒏𝝅𝒙
𝒖(𝒙, 𝒕) = ෍ 𝒄𝒏 𝒆𝒙𝒑 [− ቀ ቁ 𝒌𝒕] 𝐬𝐢𝐧 ቀ ቁ
𝒍 𝒍
𝒏=𝟏

[2]
𝟐𝒍𝟐
where 𝑐𝑛 = 𝒏𝟑 𝝅𝟑 [(𝟐 − 𝒏𝟐 𝝅𝟐 )(−𝟏)𝒏 − 𝟐]

Question Two:
b)
(i) From Taylor Series:
′ (𝑥)
ℎ2 ′′ ℎ3 ′′′ ℎ4 𝑖𝑣
𝑢(𝑥 + ℎ) = 𝑢(𝑥) + ℎ𝑢 + 𝑢 (𝑥) + 𝑢 (𝑥) + 𝑢 (𝑥) + ⋯
2! 3! 4!

′ (𝑥)
ℎ2 ′′ ℎ3 ′′′ ℎ4 𝑖𝑣
𝑢(𝑥 − ℎ) = 𝑢(𝑥) − ℎ𝑢 + 𝑢 (𝑥) − 𝑢 (𝑥) + 𝑢 (𝑥) + ⋯
2! 3! 4!
[2]
Therefore:
ℎ2 ′′ ℎ4
𝑢(𝑥 − ℎ) − 2𝑢(𝑥) + 𝑢(𝑥 + ℎ) = 2 𝑢 (𝑥) + 2 𝑢𝑖𝑣 (𝑥) + ⋯
2! 4!

𝑢(𝑥 − ℎ) − 2𝑢(𝑥) + 𝑢(𝑥 + ℎ) ℎ2 𝑖𝑣

𝑢′′ (𝑥) = − 2 𝑢 (𝑥) + ⋯
ℎ2 4!

Page | 42
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

𝝏𝟐 𝒖 ′′ (𝒙)
𝒖(𝒙 − 𝒉) − 𝟐𝒖(𝒙) + 𝒖(𝒙 + 𝒉)
= 𝒖 ≈
𝝏𝒙𝟐 𝒉𝟐
[1]
Similarly,
′ (𝑦)
𝑘 2 ′′ 𝑘 3 ′′′ 𝑘 4 𝑖𝑣
𝑢(𝑦 + 𝑘) = 𝑢(𝑦) + 𝑘𝑢 + 𝑢 (𝑦) + 𝑢 (𝑦) + 𝑢 (𝑦) + ⋯
2! 3! 4!

𝑘 2 ′′ 𝑘3 𝑘4
𝑢(𝑦 − 𝑘) = 𝑢(𝑦) − 𝑘𝑢′ (𝑦) + 𝑢 (𝑦) − 𝑢′′′ (𝑦) + 𝑢𝑖𝑣 (𝑦) + ⋯
2! 3! 4!

Therefore:
𝝏𝟐 𝒖 𝒖(𝒚 − 𝒌) − 𝟐𝒖(𝒚) + 𝒖(𝒚 + 𝒌)
𝟐
= 𝒖′′ (𝒚) ≈
𝝏𝒚 𝒌𝟐
[1]

𝜕2𝑢 𝜕2𝑢
The error term for 𝜕𝑥 2 and 𝜕𝑦 2 is of the order of ℎ2 and 𝑘 2 , respectively and the
coefficient can be obtained from the fourth derivative.
[1]

(ii) The discretization is as shown below:

Figure 16 - Discretization of the domain of the PDE

Page | 43
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

𝜕2 𝑢 𝜕2 𝑢
Given 𝜕𝑥2 + 𝜕𝑦2 = 𝑥𝑒 𝑦 , from the Central Finite Difference formulation for Second Order derivatives
𝜕2𝑢 𝜕2𝑢
and in 2bi) above:
𝜕𝑥 2 𝜕𝑦 2

𝑢𝑖−1,𝑗 − 2𝑢𝑖,𝑗 + 𝑢𝑖+1,𝑗 𝑢𝑖,𝑗−1 − 2𝑢𝑖,𝑗 + 𝑢𝑖,𝑗+1

+ = 𝑥𝑖,𝑗 𝑒 𝑦𝑖,𝑗
ℎ2 𝑘2

[2]

1 1 1 1
2 ൫𝑢𝑖−1,𝑗 + 𝑢𝑖+1,𝑗 ൯ + 2 ൫𝑢𝑖,𝑗−1 + 𝑢𝑖,𝑗+1 ൯ − 2 ( 2 + 2 ) 𝑢𝑖,𝑗 = 𝑥𝑖,𝑗 𝑒 𝑦𝑖,𝑗
ℎ 𝑘 ℎ 𝑘

𝟏 𝟏 𝟏 𝟏
𝟐( 𝟐
+ 𝟐 ) 𝒖𝒊,𝒋 = 𝟐 ൫𝒖𝒊−𝟏,𝒋 + 𝒖𝒊+𝟏,𝒋 ൯ + 𝟐 ൫𝒖𝒊,𝒋−𝟏 + 𝒖𝒊,𝒋+𝟏 ൯ − 𝒙𝒊,𝒋 𝒆𝒚𝒊,𝒋
𝒉 𝒌 𝒉 𝒌

[1]

1 25 1
ℎ = 0.4 , 𝑘 = 0.2, 2
= , = 25
ℎ 4 𝑘2
𝟏 𝟐 𝟐
𝒖𝒊,𝒋 = ൫𝒖𝒊−𝟏,𝒋 + 𝒖𝒊+𝟏,𝒋 ൯ + ൫𝒖𝒊,𝒋−𝟏 + 𝒖𝒊,𝒋+𝟏 ൯ − 𝒙 𝒆𝒚𝒊,𝒋
𝟏𝟎 𝟓 𝟏𝟐𝟓 𝒊,𝒋

[2]

(iii) The difference equations are:

1 2 2
At (𝑖, 𝑗) = (1,1): 𝑢1 = 10 (0 + 𝑢2 ) + 5 (0.4 + 𝑢5 ) − 125 (0.4)𝑒 0.2
1 2 2
At (𝑖, 𝑗) = (2,1): 𝑢2 = 10 (𝑢1 + 𝑢3 ) + 5 (0.8 + 𝑢6 ) − 125 (0.8)𝑒 0.2
1 2 2
At (𝑖, 𝑗) = (3,1): 𝑢3 = 10 (𝑢2 + 𝑢4 ) + 5 (1.2 + 𝑢7 ) − 125 (1.2)𝑒 0.2
1 2 2
At (𝑖, 𝑗) = (4,1): 𝑢4 = 10 (𝑢3 + 2𝑒 0.2 ) + 5 (1.6 + 𝑢8 ) − 125 (1.6)𝑒 0.2
1 2 2
At (𝑖, 𝑗) = (1,2): 𝑢5 = 10 (0 + 𝑢6 ) + 5 (𝑢1 + 𝑢9 ) − 125 (0.4)𝑒 0.4
1 2 2
At (𝑖, 𝑗) = (2,2): 𝑢6 = 10 (𝑢5 + 𝑢7 ) + 5 (𝑢2 + 𝑢10 ) − 125 (0.8)𝑒 0.4
1 2 2
At (𝑖, 𝑗) = (3,2): 𝑢7 = 10 (𝑢6 + 𝑢8 ) + 5 (𝑢3 + 𝑢11 ) − 125 (1.2)𝑒 0.4
1 2 2
At (𝑖, 𝑗) = (4,2): 𝑢8 = 10 (𝑢7 + 2𝑒 0.4 ) + 5 (𝑢4 + 𝑢12 ) − 125 (1.6)𝑒 0.4
1 2 2
At (𝑖, 𝑗) = (1,3): 𝑢9 = (0 + 𝑢10 ) + (𝑢5 + 𝑢13 ) − (0.4)𝑒 0.6
10 5 125
1 2 2
At (𝑖, 𝑗) = (2,3): 𝑢10 = 10 (𝑢9 + 𝑢11 ) + 5 (𝑢6 + 𝑢14 ) − 125 (0.8)𝑒 0.6
1 2 2
At (𝑖, 𝑗) = (3,3): 𝑢11 = 10 (𝑢10 + 𝑢12 ) + 5 (𝑢7 + 𝑢15 ) − 125 (1.2)𝑒 0.6
1 2 2
At (𝑖, 𝑗) = (4,3): 𝑢12 = 10 (𝑢11 + 2𝑒 0.6 ) + 5 (𝑢8 + 𝑢16 ) − 125 (1.6)𝑒 0.6
1 2 2
At (𝑖, 𝑗) = (1,4): 𝑢13 = 10 (0 + 𝑢14 ) + 5 (𝑢9 + 0.4𝑒) − 125 (0.4)𝑒 0.8
1 2 2
At (𝑖, 𝑗) = (2,4): 𝑢14 = 10 (𝑢13 + 𝑢15 ) + 5 (𝑢10 + 0.8𝑒) − 125 (0.8)𝑒 0.8

Page | 44
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

1 2 2
At (𝑖, 𝑗) = (3,4): 𝑢15 = (𝑢14 + 𝑢16 ) + (𝑢11 + 1.2𝑒) − (1.2)𝑒 0.8
10 5 125
1 2 2
At (𝑖, 𝑗) = (4,4): 𝑢16 = 10
(𝑢15 + 2𝑒 0.8 ) + 5 (𝑢12 + 1.6𝑒) − 125 (1.6)𝑒 0.8
[10]

(iv) Various approaches can be used for this question.

For Approach 1 (Matrix Inversion or Cramer’s rule)

To use Matrix inversion method to solve the system of difference equations above, they need to be put in
the form 𝐴 𝑤𝑖 = 𝑏𝑖 .

A typical MATLAB code to solve the matrix equation above is:

% Typical MATLAB Code for solving the PDE using Matrix inversion method
%% Part 2biv)
h = 0.4;
k = 0.2;
x = 0:h:2; %discrete values of x
y = 0:k:1; %discrete values of y
%Meshing
[x,y] = meshgrid(x,y); % Mesh including boundary points.

% Coefficient matrix (A) of u1 to u16

A = [1 -1/10 0 0 -2/5 0 0 0 0 0 0 0 0 0 0 0;...
-1/10 1 -1/10 0 0 -2/5 0 0 0 0 0 0 0 0 0 0;...
0 -1/10 1 -1/10 0 0 -2/5 0 0 0 0 0 0 0 0 0;...
0 0 -1/10 1 0 0 0 -2/5 0 0 0 0 0 0 0 0;...
-2/5 0 0 0 1 -1/10 0 0 -2/5 0 0 0 0 0 0 0;...
0 -2/5 0 0 -1/10 1 -1/10 0 0 -2/5 0 0 0 0 0 0;...
0 0 -2/5 0 0 -1/10 1 -1/10 0 0 -2/5 0 0 0 0 0;...
0 0 0 -2/5 0 0 -1/10 1 0 0 0 -2/5 0 0 0 0;...
0 0 0 0 -2/5 0 0 0 1 -1/10 0 0 -2/5 0 0 0;...
0 0 0 0 0 -2/5 0 0 -1/10 1 -1/10 0 0 -2/5 0 0;...
0 0 0 0 0 0 -2/5 0 0 -1/10 1 -1/10 0 0 -2/5 0;...
0 0 0 0 0 0 0 -2/5 0 0 -1/10 1 0 0 0 -2/5;...
0 0 0 0 0 0 0 0 -2/5 0 0 0 1 -1/10 0 0;...
0 0 0 0 0 0 0 0 0 -2/5 0 0 -1/10 1 -1/10 0;...
0 0 0 0 0 0 0 0 0 0 -2/5 0 0 -1/10 1 -1/10;...
0 0 0 0 0 0 0 0 0 0 0 -2/5 0 0 -1/10 1];

% Constants matrix, b
b = [(1/10)*0+(2/5)*0.4-(2/125)*0.4*exp(0.2); (2/5)*0.8-(2/125)*0.8*exp(0.2);...
(2/5)*1.2-(2/125)*1.2*exp(0.2); (1/10)*2*exp(0.2)+(2/5)*1.6-
(2/125)*1.6*exp(0.2);...
(1/10)*0-(2/125)*0.4*exp(0.4); -(2/125)*0.8*exp(0.4);...
-(2/125)*1.2*exp(0.4); (1/10)*2*exp(0.4)-(2/125)*1.6*exp(0.4);...
(1/10)*0-(2/125)*0.4*exp(0.6); -(2/125)*0.8*exp(0.6);...
-(2/125)*1.2*exp(0.6); (1/10)*2*exp(0.6)-(2/125)*1.6*exp(0.6);...
(1/10)*0+(2/5)*0.4*exp(1)-(2/125)*0.4*exp(0.8); (2/5)*0.8*exp(1)-
(2/125)*0.8*exp(0.8);...
(2/5)*1.2*exp(1)-(2/125)*1.2*exp(0.8); (1/10)*2*exp(0.8)+(2/5)*1.6*exp(1)-
(2/125)*1.6*exp(0.8)];

Page | 45
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

%Matrix inversion and solution

u = A\b; %Inverse matrix method, u is u1 to u16
disp(u)

%% Part 2bv)
%Comparison of numerical and analytical solutions
n_x_int = length(x) - 2; %number of internal points along x
n_y_int = length(y) - 2; %number of internal points along y
%Exact Analytical solution
u_exact = @(x,y) x.*exp(y);
%Initialization of table
u_int = zeros(n_x_int*n_y_int,7);

u_int(:,1) =repmat((1:n_x_int)',n_y_int,1); % Values of i

r=repmat((1:n_y_int)',1,n_x_int)';
u_int(:,2) = r(:); % Values of j
u_int(:,3) =repmat((x(1,2:n_x_int+1))',n_y_int,1); % values of x
r_2=repmat(y(2:n_y_int+1,1),1,n_x_int)';
u_int(:,4) = r_2(:); % values of y
u_int(:,5) = u; % numerical results

u_int(:,6) = feval(u_exact,u_int(:,3),u_int(:,4)); % exact analytical results

%percentage error between analytical and numerical solutions

u_int(:,7) = ((u_int(:,5)-u_int(:,6))./u_int(:,6))*100; %percentage error
disp(u_int)
%%

This gives 𝒖𝟏 to 𝒖𝟏𝟔 as shown below:

0.4887
0.9774
1.4661
1.9546
0.5970
1.1939
1.7908
2.3875
0.7291
1.4582
2.1872
2.9161
0.8904
1.7808
2.6712
3.5614

Page | 46
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

For Approach 2 (Non-iterative approach)

This is the case if the coded solution of the PDE and its boundary conditions is attempted from
scratch in MATLAB using a non-iterative scheme.

One may write non-iterative MATLAB code that solves the PDE and its boundary conditions from
scratch. Typical MATLAB code for this is as shown below. It produces the same results for 𝒖𝟏 to
𝒖𝟏𝟔 as those in the previous approach (as copied below from top-left to bottom-right).
0.4887 0.9774 1.4661 1.9546
0.5970 1.1939 1.7908 2.3875
0.7291 1.4582 2.1872 2.9161
0.8904 1.7808 2.6712 3.5614

% Typical MATLAB Code to solve boundary-value PDEs using the Finite Difference
...Method using Non-Iterative Scheme
% Numerical scheme used is an implicit second order central difference
...in space(5-point difference)
%Domain is [x_start,x_end] x [y_start, y_end]

% u is the numerical solution of the PDE at the mesh points.

%% 2biv)
% Inputs
%Specifying parameters
h=0.4; %Width of step along x axis (mesh spacing)
k=0.2; %Width of step along y axis (mesh spacing)
x_start = 0; %Start value of x
x_end = 2; %End value of x
y_start = 0; %Start value of y
y_end = 1; %End value of y

%Right hand side function of the PDE equation (D2u/Dx2 + D2u/Dy2 = f)

%If the RHS is a constant, write it as constant*x.^0.*y.^0 to create a
%vector later
f = @(x,y)x.*exp(y); %The RHS of PDE equation

%Dirichlet Boundary conditions, g_b, g_t, g_l and g_r

%If the Boundary Condition is a constant, write it as constant*x.^0.*y.^0 to create
a
%vector later
g_b=@(x,y)x; %boundary condition at y=y_start, BOTTOM
g_t=@(x,y)x*exp(1); %boundary condition at y=y_end, TOP
g_l=@(x,y)0*y; %boundary condition at x=x_start, LEFT
g_r=@(x,y)2*exp(y); %boundary condition at x=x_end, RIGHT

%% EXACT ANALYTICAL RESULT

%This is optional. It is to compare the numerical result with the exact
%analytical results
u_exact_available = 1; %Set to 1 if exact analytical result is available
% Set to 0 if exact analytical result is not available and comment off next

Page | 47
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

% line
u_exact = @(x,y)x.*exp(y);
%%
%Mesh generation
nx = (x_end - x_start)/h + 1; %Number of grid points along x axis
ny = (y_end - y_start)/k + 1; %Number of grid points along y axis
n_x_int = nx-2; %Number of interior grid points along x axis
n_y_int = ny-2; %Number of interior grid points along y axis
x=x_start:h:x_end; %x values for the grid points
y=y_start:k:y_end; %y values for the grid points
%Meshing
[x,y] = meshgrid(x,y); % Mesh including boundary points.

%%
% Computation of u(x,y) on the interior boundary grid points from the Dirichlet
boundary condition
ub = zeros(n_y_int,n_x_int); %Initialization n=n_x_int
idx = 2:n_x_int+1; %Indices of interior nodes along x-axis

idy = 2:n_y_int+1; %Indices of interior nodes along y-axis

%% Boundaries
% Components of the Central Difference terms for D2u/Dx2 contributed by
% the values of the boundary conditions on the left and right boundary nodes
% which will be transferred to the right hand side of the equations later

% Left and Right boundaries

ub(:,1) = (1/h^2)*feval(g_l,x(2,1),y(idy,1));% Left Boundary
ub(:,n_x_int) = (1/h^2)*feval(g_r,x(2,n_x_int+2),y(idy,1));% Right Boundary

% Top and Bottom boundaries

ub(1,1:n_x_int) = ub(1,1:n_x_int) + (1/k^2)*feval(g_b,x(1,idx),y(1,2)); %Bottom
Boundary
ub(n_y_int,1:n_x_int) = ub(n_y_int,1:n_x_int) +
(1/k^2)*feval(g_t,x(1,idx),y(n_y_int+2,2)); %Top Boundary

%%
% Convert ub to a vector using column reordering (reshape)
% for it to be subtracted from the column vector of the RHS of the PDE
ub = reshape(ub,n_x_int*n_y_int,1);
% Evaluate the RHS of the PDE at the interior points.
uf = feval(f,x(idy,idx),y(idy,idx));
% Convert f to a vector using column reordering
uf = reshape(uf,n_x_int*n_y_int,1);

% Create the D2x and D2y matrices

z_x = [-2;1;zeros(n_x_int-2,1)]; %coefficients of u(x,y) at the internal nodes
z_y = [-2;1;zeros(n_y_int-2,1)]; %coefficients of u(x,y) at the internal nodes

D2x = (1/h^2)*kron(toeplitz(z_x,z_x),eye(n_y_int)); %kron creates a block diagonal

matrix
D2y = (1/k^2)*kron(eye(n_x_int),toeplitz(z_y,z_y));
%if A is m*n and B is p*q, kron(A,B) creates m*p-by-n*q diagonal block
%matrix
%eye(n) creates a n*n Identity matrix

Page | 48
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

%toeplitz(c,r) returns a nonsymmetric Toeplitz matrix with c as its first column

%and r as its first row. Toeplitz matrix or diagonal-constant matrix is a
%matrix in which each descending diagonal from left to right is constant
%%
% Solve the system
u = (D2x + D2y)\(uf-ub); %Inverse matrix method to calculate the unknown u(i,j)
values
% Convert u from a column vector to a matrix to make it easier towork with
% for plotting.
u = reshape(u,n_y_int,n_x_int);
disp(u)

%% 2bv)
%%Create the Table for the results
if u_exact_available == 1
%Exact analytical solution is available for comparison
u_int = zeros(n_x_int*n_y_int,7); %Initialization of table

u_int(:,1) =repmat((1:n_x_int)',n_y_int,1);
r=repmat((1:n_y_int)',1,n_x_int)';
u_int(:,2) = r(:);
u_int(:,3) =repmat((x(1,2:n_x_int+1))',n_y_int,1);
r_2=repmat(y(2:n_y_int+1,1),1,n_x_int)';
u_int(:,4) = r_2(:);
r_3 = u';

u_int(:,5) = r_3(:);
u_int(:,6) = feval(u_exact,u_int(:,3),u_int(:,4));
%percentage error between analytical and numerical solutions
u_int(:,7) = ((u_int(:,5)-u_int(:,6))./u_int(:,6))*100; %percentage error
disp(u_int)

else
%Exact analytical solution is NOT available for comparison
u_int = zeros(n_x_int*n_y_int,5); %Initialization of table

u_int(:,1) =repmat((1:n_x_int)',n_y_int,1);
r=repmat((1:n_y_int)',1,n_x_int)';
u_int(:,2) = r(:);
u_int(:,3) =repmat((x(1,2:n_x_int+1))',n_y_int,1);
r_2=repmat(y(2:n_y_int+1,1),1,n_x_int)';
u_int(:,4) = r_2(:);
u_int(:,5) = u(:);
disp(u_int)
end
%%
Alternative Approach 3 (Iterative approach)
From 2bii,
𝑢𝑖−1,𝑗 − 2𝑢𝑖,𝑗 + 𝑢𝑖+1,𝑗 𝑢𝑖,𝑗−1 − 2𝑢𝑖,𝑗 + 𝑢𝑖,𝑗+1
+ = 𝑥𝑖,𝑗 𝑒 𝑦𝑖,𝑗
ℎ2 𝑘2
This gives:

Page | 49
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

1
𝑢𝑖,𝑗 = [𝑘 2 ൫𝑢𝑖−1,𝑗 + 𝑢𝑖+1,𝑗 ൯ + ℎ2 ൫𝑢𝑖,𝑗−1 + 𝑢𝑖,𝑗+1 ൯ − ℎ2 𝑘 2 𝑥𝑖,𝑗 𝑒 𝑦𝑖,𝑗 ]
2(ℎ2 + 𝑘2)

One might write iterative MATLAB code that solves the PDE and its boundary conditions from
first principles. Typical MATLAB code for this is shown below. It produces similar results for 𝒖𝟏
to 𝒖𝟏𝟔 as those in the previous approach (as copied below from top-down from left to right).
0.4887 0.5970 0.7291 0.8904
0.9774 1.1939 1.4582 1.7808
1.4661 1.7908 2.1872 2.6711
1.9546 2.3875 2.9160 3.5614

% Typical Matlab Code for solving the PDE using Iterative scheme
%% Part 2biv)
%Specifying parameters
h=0.4; %Width of step along x axis (mesh spacing)
k=0.2; %Width of step along y axis (mesh spacing)
x_start = 0; %Start value of x

x_end = 2; %End value of x

y_start = 0; %Start value of y
y_end = 1; %End value of y

x=x_start:h:x_end; %x values for the grid points

y=y_start:k:y_end; %y values for the grid points

nx = (x_end - x_start)/h + 1; %Number of grid points along x axis

ny = (y_end - y_start)/k + 1; %Number of grid points along y axis

tolerance=1.0e-5; %Tolerance OR Stopping criteria for errors

%Initialising solutions
u=zeros(nx,ny); %Zero initial values of (n+1)th or Current
un=zeros(nx,ny); %Zero initial values of (n)th iteration or Previous

%Applying boundary conditions

u(1,:)=0; %boundary condition at x=x_start, LEFT
u(nx,:)=2*exp(y); %boundary condition at x=x_end, RIGHT
u(:,1)=x; %boundary condition at y=y_start, BOTTOM
u(:,ny)=x*exp(1); %boundary condition at y=y_end, TOP

%Iterative scheme
max_error = 1; %initializing the error value to allow the while loop to start
while( max_error > tolerance )
un=u;
for i = 2:nx-1
for j = 2:ny-1

Page | 50
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

u(i,j)=1/(2*(h^2+k^2))*(k^2*(un(i-1,j)+un(i+1,j))+h^2*(un(i,j-
1)+un(i,j+1))...
-h^2*k^2*(x(i)*exp(y(j))));
end
end
error = u-un;
max_error = max(error(:)); %to obtain maximum error of the u values
end
disp(u(2:nx-1,2:ny-1))

%% Part 2bv)
%Comparison of numerical and analytical solutions
n_x_int = length(x) - 2; %number of internal points along x
n_y_int = length(y) - 2; %number of internal points along y

%Exact Analytical solution

u_exact = @(x,y) x.*exp(y);
%Initialization of table

u_int = zeros(n_x_int*n_y_int,7);

u_int(:,1) =repmat((1:n_x_int)',n_y_int,1); % Values of i

r=repmat((1:n_y_int)',1,n_x_int)';
u_int(:,2) = r(:); % Values of j
u_int(:,3) =repmat((x(1,2:n_x_int+1))',n_y_int,1); % values of x
r_2=repmat(y(2:n_y_int+1),n_x_int,1);
u_int(:,4) = r_2(:); % values of y

r_3 = u(2:nx-1,2:ny-1);
u_int(:,5) = r_3(:); % numerical results
u_int(:,6) = feval(u_exact,u_int(:,3),u_int(:,4)); % exact analytical results
%percentage error between analytical and numerical solutions
u_int(:,7) = ((u_int(:,5)-u_int(:,6))./u_int(:,6))*100; %percentage error
disp(u_int)
%%

(v)
The numerical results compared to the analytical solutions for the PDE is as summarized in Table 1
below. The analytical and numerical results are in close agreement with the highest percentage
error being 0.041649%. The observed negligible errors are due to the error terms in the
𝜕2𝑢 𝜕2 𝑢
Taylor series expansion for 𝜕𝑥 2 and 𝜕𝑦 2 from the third and higher order derivatives of 𝑢(𝑥, 𝑦).
[2]

Page | 51
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Table 1 applies to both Alternative Approaches 1 and 2.

Table 1: Summary of Numerical and Analytical results of the PDE.

i j 𝑥𝑖 𝑦𝑗 Numerical Exact Solution % error
𝑢𝑖,𝑗 𝑢𝑖,𝑗 = 𝑥𝑖 𝑒 𝑦𝑗
1 1 0.4 0.2 0.488717 0.488561 0.031834
2 1 0.8 0.2 0.977422 0.977122 0.030725
3 1 1.2 0.2 1.46609 1.465683 0.027776
4 1 1.6 0.2 1.954638 1.954244 0.020148
1 2 0.4 0.4 0.596978 0.59673 0.041649
2 2 0.8 0.4 1.193939 1.19346 0.040172
3 2 1.2 0.4 1.790838 1.79019 0.03623
4 2 1.6 0.4 2.387541 2.38692 0.026045
1 3 0.4 0.6 0.729114 0.728848 0.036535
2 3 0.8 0.6 1.45821 1.457695 0.035312
3 3 1.2 0.6 2.187243 2.186543 0.032018
4 3 1.6 0.6 2.91607 2.91539 0.023312
1 4 0.4 0.8 0.890408 0.890216 0.021475
2 4 0.8 0.8 1.780804 1.780433 0.020849
3 4 1.2 0.8 2.67116 2.670649 0.019137
4 4 1.6 0.8 3.561379 3.560865 0.01441
(within the code) [3]

The comparison of the numerical solution with the analytical solution is as shown in Table 2
below. Table 2 applies to the Alternative Approach 3.
Table 2: Summary of Numerical and Analytical results of the PDE (Iterative Scheme).
i j 𝑥𝑖 𝑦𝑗 Numerical Exact Solution % error
𝑦𝑗
𝑢𝑖,𝑗 𝑢𝑖,𝑗 = 𝑥𝑖 𝑒
1 1 0.4 0.2 0.488702 0.488561 0.028869
2 1 0.8 0.2 0.977399 0.977122 0.028352
3 1 1.2 0.2 1.466067 1.465683 0.026175
4 1 1.6 0.2 1.954624 1.954244 0.019414
1 2 0.4 0.4 0.596955 0.59673 0.037765
2 2 0.8 0.4 1.193901 1.19346 0.036992
3 2 1.2 0.4 1.790801 1.79019 0.034133
4 2 1.6 0.4 2.387518 2.38692 0.025062
1 3 0.4 0.6 0.72909 0.728848 0.033319
2 3 0.8 0.6 1.458172 1.457695 0.032738
3 3 1.2 0.6 2.187205 2.186543 0.030281
4 3 1.6 0.6 2.916046 2.91539 0.022515

Page | 52
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

1 4 0.4 0.8 0.890393 0.890216 0.019866

2 4 0.8 0.8 1.78078 1.780433 0.019532
3 4 1.2 0.8 2.671137 2.670649 0.018269
4 4 1.6 0.8 3.561364 3.560865 0.014003

ENGF0004 Example: PDEs

1. Elliptic PDE: Solve the Laplace equation 𝑢𝑥𝑥 + 𝑢𝑦𝑦 = 0 in the rectangle 0 < 𝑥 < 𝑎, 0 < 𝑦 < 𝑏 with the
𝜕𝑢 𝜕𝑢
boundary conditions: At 𝑥 = 0, 𝜕𝑥 = 0; At 𝑥 = 𝑎, 𝜕𝑥 = 0; At 𝑦 = 0, 𝑢 = 0; At 𝑦 = 𝑏, 𝑢 = 𝑓(𝑥), where f is a

given arbitrary function. The problem is shown schematically below. The desired answer is therefore the
specification of 𝑢(𝑥, 𝑦) for all points within and on the given rectangle.

SOLUTION:
Apply the separation of variables method: 𝑢(𝑥, 𝑦) = 𝑔(𝑥)ℎ(𝑦). Then the PDE becomes
𝑔′′ (𝑥) ℎ ′′ (𝑦) 𝑔′′ (𝑥) ℎ ′′ (𝑦)
𝑔′′ (𝑥) ∙ ℎ(𝑦) + 𝑔(𝑥) ∙ ℎ′′ (𝑦) = 0, or 𝑔(𝑥)
+ ℎ(𝑦)
= 0, or 𝑔(𝑥)
=− ℎ(𝑦)
= 𝑘,constant.

We set both sides of the PDE equal to a constant k, because the only way for a function of x to equal a
different function of y (recall x and y are independent variables) is for both function to be constant.
The real constant k can be positive, zero, or negative. We consider each case in turn.

Case 1, 𝒌 > 𝟎:

Page | 53
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Then 𝑔′′ (𝑥) = 𝑘𝑔(𝑥), and 𝑔(𝑥) = 𝐴𝑒 𝑃𝑥 + 𝐵𝑒 −𝑃𝑥 , where 𝑃 = ξ𝑘 > 0.

𝜕𝑢
Apply the BC at 𝑥 = 0: 0 = 𝜕𝑥 (0, 𝑦) = 𝑔′ (0) ∙ ℎ(𝑦) = 𝐴𝑒 𝑃∙0 + 𝐵𝑒 −𝑃∙0 or 0 = 𝑃(𝐴 − 𝐵) ∙ ℎ(𝑦). By

assumption, 𝑃 ≠ 0, leading to 𝐴 = 𝐵, and 𝑔(𝑥) = 𝐴(𝑒 𝑃𝑥 + 𝑒 −𝑃𝑥 ). [Note that ℎ(𝑦) = 0 can be rejected
because it leads to the trivial solution 𝑢 ≡ 0.]

𝜕𝑢
Next, apply the BC at 𝑥 = 𝑎: 0 = 𝜕𝑥 (𝑎, 𝑦) = 𝑔′ (𝑎) ∙ ℎ(𝑦) = 𝐴𝑃(𝑒 𝑃∙𝑎 + 𝑒 −𝑃∙𝑎 )ℎ(𝑦). Rejecting ℎ(𝑦) = 0

and 𝑃 = 0 again, there remain 𝐴 = 0 or the term in parentheses equal 0. If 𝐴 = 0, 𝑔(𝑥) = 0, leading to
𝑢(𝑥, 𝑦) = 0, which can be rejected. For real P, the term in parentheses can only equal 0 for 𝑃 = 0, which
would contradict the assumption of Case 1, that 𝑃 > 0 because 𝑘 > 0.

We conclude that both BCs cannot be satisfied. Hence Case 1 does not contribute any solutions.

Case 2, 𝒌 = 𝟎:
Then 𝑔′′ (𝑥) = 0, or 𝑔(𝑥) = 𝐸𝑥 + 𝐹, where E and F are integration constants whose values have to be
determined.
𝜕𝑢
Apply the BC at 𝑥 = 0, 0 = (0, 𝑦) = 𝑔′ (0) ∙ ℎ(𝑦) = 𝐸ℎ(𝑦). As before, we reject ℎ(𝑦) = 0 and conclude
𝜕𝑥

that 𝐸 = 0 and therefore that 𝑔(𝑥) = 𝐹,.

𝜕𝑢
Apply the BC at 𝑥 = 𝑎, 0 = 𝜕𝑥 (𝑎, 𝑦) = 𝑔′ (𝑎) ∙ ℎ(𝑦) = 0 ∙ ℎ(𝑦), no further information is gained.

So Case 2 contributes the solution 𝑔(𝑥) = 𝐹, where the value of F has not yet been determined.

Case 3, 𝒌 < 𝟎:
Then 𝑔′′ (𝑥) = 𝑘𝑔(𝑥) = −𝑃2 𝑔(𝑥), where 𝑃 = ξ−𝑘 > 0. Hence 𝑔(𝑥) = 𝐿 cos(𝑃𝑥) + 𝑀 sin(𝑃𝑥), where the
integration constants L and M have to be determined.
𝜕𝑢
Apply the BC at 𝑥 = 0: 0 = (0, 𝑦) = 𝑔′ (0) ∙ ℎ(𝑦) = ℎ(𝑦){−𝐿𝑃sin(0) + 𝑀𝑃cos(0)} = 𝑀𝑃ℎ(𝑦). As before,
𝜕𝑥

we reject ℎ(𝑦) = 0 and 𝑃 = 0 to conclude that 𝑀 = 0, resulting in 𝑔(𝑥) = 𝐿 cos(𝑃𝑥).

𝜕𝑢
Apply the BC at 𝑥 = 𝑎, 0 = 𝜕𝑥 (𝑎, 𝑦) = −𝐿𝑃sin(𝑃𝑎)ℎ(𝑦). As before we reject ℎ(𝑦) = 0 and 𝐿 = 0 to
𝑛𝜋
conclude that sin(𝑃𝑎) = 0, or 𝑃𝑛 𝑎 = 𝑛𝜋, 𝑛 = 1, 2,3, …or 𝑃𝑛 = 𝑎
,𝑛 = 1, 2,3, … .

This set of solutions to 𝑃𝑛 (which are called eigenvalues) leads to the most general solution
𝑛=1 𝐿𝑛 cos (𝑃𝑛 𝑥), where the integrity of constants 𝐿𝑛 have to be determined.
𝑔(𝑥) = ∑∞

Page | 54
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

We now turn to the solution of ℎ(𝑦), noting that only 𝑘 = 0 and 𝑘 < 0 need to be considered. For 𝒌 = 𝟎,
ℎ′′ (𝑦) = 0, resulting in ℎ(𝑦) = 𝛼𝑦 + 𝛽, where 𝛼 and 𝛽 are integration constants.
Apply the BC at 𝑦 = 0: 0 = 𝑢(𝑥, 0) = 𝑔(𝑥)ℎ(0) = 𝑔(𝑥)𝛽. Rejecting 𝑔(𝑥) = 0 for reasons given earlier,
we conclude 𝛽 = 0, and ℎ(𝑦) = 𝛼𝑦.

So 𝑢(𝑥, 𝑦) = 𝑔(𝑥)ℎ(𝑏) = 𝐹𝛼𝑦, or

𝑢(𝑥, 𝑦) = 𝑇𝑦 (1).

where 𝑇 = 𝐹𝛼 is a constant. We will return to Equation (1) below.

𝑛𝜋
For the case 𝒌 < 𝟎, ℎ′′ (𝑦) = −𝑘ℎ(𝑦). Since −𝑘 > 0, the solution is 𝑅𝑛 𝑒 𝑃𝑛𝑦 + 𝑆𝑛 𝑒 −𝑃𝑛𝑦 , where 𝑃𝑛 = ξ−𝑘 = , as
2𝑎

before, and 𝑅𝑛 and 𝑆𝑛 are integration constants.

Apply the BC at 𝑦 = 0:
0 = 𝑢(𝑥, 0) = 𝑔(𝑥)ℎ(0) = 𝑔(𝑥)(𝑅𝑛 +𝑆𝑛 ). We conclude that 𝑅𝑛 = −𝑆𝑛 , and

ℎ(𝑦) = 𝑅𝑛 (𝑒 𝑃𝑛𝑦 − 𝑒 −𝑃𝑛𝑦 ) = 2𝑅𝑛 sinh (𝑃𝑛 𝑦).

So the general solution so far is
𝑢(𝑥, 𝑦) = 𝑇𝑦 + ∑∞
𝑛=1 𝐿𝑛 cos (𝑃𝑛 𝑥). 2𝑅𝑛 sinh (𝑃𝑛 𝑦) (2)

We are now ready to apply the BC at 𝑦 = 𝑏:

∞

𝑢(𝑥, 𝑏) = 𝑓(𝑥) = 𝑇𝑏 + ෍ 𝐿𝑛 cos (𝑃𝑛 𝑥). 2𝑅𝑛 sinh (𝑃𝑛 𝑏)

𝑛=1

Noting that everything but the cosine terms are independent of x, we define

𝑄0
2
= 𝑇𝑏, 𝑄𝑛 = 2𝐿𝑛 𝑅𝑛 sinh (𝑃𝑛 𝑏), (3)

and arrive at
∞
𝑄0
𝑓(𝑥) = + ෍ 𝑄𝑛 cos (𝑃𝑛 𝑥)
2
𝑛=1

But this is a classical Fourier cosine series, and we know that the coefficients Q are defined by:

2 𝑎 𝑛𝜋𝑥
𝑄𝑛 = 𝑎 ∫0 𝑓(𝑥) cos ቀ 𝑎
ቁ 𝑑𝑥 , 𝑛 = 0,1,2, … (4)

Note the n=0 term (the ‘constant’ term with respect to x), which is characteristic of Fourier cosine series.
We conclude that

Page | 55
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

∞ 𝑛𝜋𝑦
𝑄0 𝑦 𝑛𝜋𝑥 sinh ቀ 𝑎 ቁ
𝑢(𝑥, 𝑦) = + ෍ 𝑄𝑛 cos ቀ ቁ (5)
2𝑏 𝑎 𝑛𝜋𝑏
𝑛=1 sinh ቀ 𝑎 ቁ

with the coefficients 𝑄𝑛 being specified by Equation (4).

2. Parabolic PDE: A long rod (which can therefore be considered one-dimensional) is subjected to an
initial temperature distribution along its axis. One end of the rod is maintained at a constant
temperature, while the other end is kept insulated, i.e., there is no heat transfer from the rod into the
surroundings at that end point. The heat equation describes this situation: 𝑢𝑡 = 𝑘𝑢𝑥𝑥 , for 0 < 𝑥 <
𝐿, 0 < 𝑡 < ∞, where u(x,t) is the temperature at a point x in the rod at time t. The relevant initial and
boundary conditions are:
At t = 0, 𝑢(𝑥, 0) = 𝑓(𝑥), where f(x) is a given but arbitrary function;
At 𝑥 = 0, 𝑢(0, 𝑡) = 0;
𝜕𝑢
At 𝑥 = 𝐿, 𝜕𝑥
(𝐿, 𝑡) = 0. As we’ve discussed in class, remember that this means that the derivative is
taken first, and the limit second.
Solve for the temperature u(x,t), which gives the temperature at any point x at any time t.

SOLUTION:
Because one variable (x) has a finite domain of definition, try the separation of variables. As usual, assume
𝑢(𝑥, 𝑡) = 𝑔(𝑥)ℎ(𝑡). Substituting this into the PDE gives 𝑔(𝑥). ℎ′ (𝑡) = 𝑘. 𝑔′′(𝑥). ℎ(𝑡), or
𝑔′′ (𝑥) ℎ ′ (𝑡)
𝑔(𝑥)
= 𝑘.ℎ(𝑡) = −𝑃2 , where −𝑃2 is a constant because the function of x on the left-hand side can equal an
independent function of t on the right only by both being equal to a constant. The choice of the negative
constant is because the first-order ODE in time will produce an exponential solution, as we will see shortly,
and the temperature cannot increase exponentially with time, because it will become unbounded. So a
negative exponential solution is needed to be physically reasonable, leading to the choice of −𝑃2 . (The
squared form is only for convenience, as will become clear in what follows.)
The solution for g(x) must satisfy 𝑔′′ (𝑥) = −𝑃2 . 𝑔(𝑥), which implies 𝑔(𝑥) = 𝐿 cos(𝑃𝑥) + 𝑀 sin(𝑃𝑥), where the
integration constants L and M have to be determined.

Applying the first BC at x = 0 leads to 0 = 𝐿 + 0, or 𝐿 = 0. So 𝑔(𝑥) = 𝑀 sin(𝑃𝑥).

Applying the second BC at x = L leads to 0 = 𝑀𝑃 cos(𝑃𝐿). Since M = 0 and P = 0 lead to the trivial solution u =
0, we conclude that cos(𝑃𝐿) = 0, i.e., there are an infinity of solutions for

2𝑛−1 𝜋
𝑃𝑛 = 2
.𝐿 ,𝑛 = 1,2,3, …

of the form 𝑔𝑛 (𝑥) = 𝑅𝑛 sin(Pn x), where the 𝑅𝑛 are integration constants.

Page | 56
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

2𝑛−1 𝜋 2
For each n, the solution of h(t) satisfies ℎ′ (𝑡) = −𝑘 ቀ 2
. 𝐿 ቁ ℎ(𝑡). This linear first-order ODE can be solved
2𝑛−1 𝜋 2
immediately to give ℎ𝑛 (𝑡) = 𝑆𝑛 . exp {−𝑘 ቀ 2 . 𝐿 ቁ 𝑡}, where the 𝑆𝑛 are integration constants.
The general solution is therefore
∞ 2n − 1 π 2𝑛 − 1 𝜋 2
𝑢(𝑥, 𝑡) = ෍ 𝑄𝑛 . sin (. x). exp ቊ−𝑘 ( . ) 𝑡ቋ (1)
𝑛=1 2 L 2 𝐿
where 𝑄𝑛 = 𝑅𝑛 . 𝑆𝑛 has combined the integration constants. These constants are evaluated by applying the
2n−1 π
initial condition: 𝑓(𝑥) = 𝑢(𝑥, 0) = ∑∞
𝑛=1 𝑄𝑛 . sin ቀ 2
. L xቁ. This is clearly a Fourier series, with the
coefficients given by:
2 𝐿 2n − 1 π
𝑄𝑛 = න 𝑓(𝑥). sin ( . x) 𝑑𝑥 (2)
𝐿 0 2 L

So, the full solution is given by Equation (1), with the coefficients 𝑄𝑛 being determined by (2).

3. Hyperbolic PDE: Consider the wave equation 𝑦𝑡𝑡 = 𝑎2 𝑦xx , for 0 < 𝑥 < 𝐿, 0 < 𝑡 < ∞. This equation
describes the elastic oscillations of a bar of length L, where y(x,t) represents the displacement of the point
x in the bar as a function of time t. The relevant initial and boundary conditions are:
At t = 0, 𝑦(𝑥, 0) = 𝑦𝑡 (𝑥, 0) = 0;
At 𝑥 = 0, 𝑦(0, t) = 0;
𝐹
At 𝑥 = 𝐿, 𝑦𝑥 (𝐿, 𝑡) = 𝐸;
where 𝑎, F and E are positive constants. Solve for the displacement of the end of the bar, y(L,t), as a
function of time. [Note that y(L,t) is different from y(x,t)!]

SOLUTION:
Use the Laplace transform on t to give ℒ{𝑦(𝑥, 𝑡)} = 𝑌(𝑥, 𝑠).

Then the PDE becomes 𝑠 2 𝑌(𝑥, 𝑠) − 0 − 0 = 𝑎2 𝑌𝑥𝑥 (𝑥, 𝑠)

𝑠2
∴ 𝑌𝑥𝑥 (𝑥, 𝑠) = 𝑌(𝑥, 𝑠)
𝑎2
The solution to this linear homogeneous ODE is
𝑠 𝑠
𝑌(𝑥, s) = A cosh ቀ 𝑥ቁ + 𝐵 sinh ቀ 𝑥ቁ
𝑎 𝑎
Apply the BC at 𝑥 = 0 in the Laplace-transformed space:
𝑌(0, s) = 0 = A
𝑠
∴ 𝑌(𝑥, 𝑠) = B sinh ቀ 𝑥ቁ
𝑎
𝐹
The BC at 𝑥 = L becomes 𝑌𝑥 (L, s) = 𝐸𝑠

Page | 57
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

𝐹 𝑠 𝑠
∴ = 𝐵 ∙ ∙ cosh ቀ 𝐿ቁ
𝐸𝑠 𝑎 𝑎
𝐹𝑎 1 1
∴B= ∙ ∙ ,
𝐸 𝑠2 coshቀ 𝑠 𝐿ቁ
𝑎

resulting in the solution in the transformed space:

𝑠
𝐹𝑎 1 sinhቀ𝑎𝑥ቁ
𝑌(𝑥, s) = ∙ ∙ (1)
𝐸 𝑠2 coshቀ 𝑠 𝐿ቁ
𝑎

It follows that the displacement of the end of the bar is

𝑠
𝐹𝑎 1 sinhቀ𝑎𝐿ቁ 𝐹𝑎 1 𝑠
𝑌(L, t) = ∙ ∙
𝐸 𝑠2 coshቀ 𝑠 𝐿ቁ
= ∙
𝐸 𝑠2
∙ tanh ቀ𝑎 𝐿ቁ (2)
𝑎

Referring to Laplace transform tables, the triangular wave H(t, b) shown in the figure below has the Laplace
transform:
1 𝑏𝑠
ℒ{H(t, b)} = 𝑠2 tanh ቀ 2 ቁ (3)

Comparing Equation (2) and (3), we conclude that

𝐹𝑎 2𝐿
𝑦(L, t) = ∙ 𝐻 (𝑡, )
𝐸 𝑎

The end of the bar therefore also exhibits a triangular wave pattern in time:

Page | 58
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Note: do not confuse this triangular H-function with the Heaviside step function, which is also sometimes
denoted by H and sometimes by u.

Page | 59
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Appendix: Fourier series

Lecture 1
∞
𝑎0
(1) 𝑓(𝑥) = + ෍ 𝑎𝑟 cos(𝑟𝑥) + 𝑏𝑟 sin(𝑟𝑥)
2
𝑟=1

𝑓(𝑥) = 𝑓(𝑥 + 2𝑘𝜋), where 𝑘 is an integer or function is periodic of principal period 2𝜋.

Consequently, when we define such function, we touch at intervals of 2𝜋, typically −𝜋 < 𝑥 ≤ 𝜋.

We will show that

1 𝜋
𝑎0 = න 𝑓(𝑥)𝑑𝑥
𝜋 −𝜋
1 𝜋
𝑎𝑟 = න 𝑓(𝑥) cos(𝑟𝑥) 𝑑𝑥 𝑟 = 1,2,3 …
𝜋 −𝜋
1 𝜋
𝑏𝑟 = න 𝑓(𝑥) sin(𝑟𝑥) 𝑑𝑥
𝜋 −𝜋

lim 𝑓(𝑥) = 1
𝑥→𝑎 − 1

lim 𝑓(𝑥) = 1⁄2

𝑥→𝑎 +
1/2
(The left hand limit ≠ right hand limit)

∞ 𝑎 2𝑎 3𝑎
𝑎0
𝑓(𝑥) = + ෍ 𝑎𝑟 cos(𝑟𝑥) + 𝑏𝑟 sin(𝑟𝑥)
2
𝑟=1

Dirichlet's theorem in the next lecture

Fourier coefficients we will justify the form of the Fourier coefficients provided we are aware of the
following results:
𝜋 0, 𝑟≠𝑠
(1) න cos(𝑟𝑥)sin(𝑠𝑥)𝑑𝑥 = { 2𝜋, 𝑟=𝑠=0
−𝜋 𝜋, 𝑟=𝑠>0

𝜋 0, 𝑟≠𝑠
(2) න sin(𝑟𝑥)sin(𝑠𝑥)𝑑𝑥 = { 0, 𝑟=𝑠=0
−𝜋 𝜋, 𝑟=𝑠>0

Page | 60
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

𝜋
(3) න sin(𝑟𝑥)cos(𝑠𝑥)𝑑𝑥 = {0, ∀𝑟&𝑠
−𝜋

𝜋
0, 𝑟>0
(4) න cos(𝑟𝑥)𝑑𝑥 = {
−𝜋
2𝜋, 𝑟=0

𝜋
(5) න sin(𝑟𝑥)𝑑𝑥 = 0, ∀𝑟
−𝜋

𝜋 𝜋
cos(𝑟𝑥)
(5) න sin(𝑟𝑥)𝑑𝑥 = − [ ]
−𝜋 𝑟 −𝜋
𝑥
1
= − [cos(𝜋) − cos(𝜋)] = 0 −𝑥
r

cos(𝑥) is an even function. A function 𝑓 (𝑥) is even iff

cos(𝑥) = cos(−𝑥)
𝑓(𝑥) = 𝑓(−𝑥)

Similarly, a function 𝑓 − (𝑥) is odd iff

𝑓(𝑥) = −𝑓(−𝑥)
ℎ
𝜋 𝑦
(4) න cos(𝑟𝑥)𝑑𝑥 if 𝑟 = 0, cos(0) = 1
−𝜋
𝑥

𝜋
න d𝑥 = [x]𝜋−𝜋 = 𝜋 + 𝜋 = 2𝜋 −𝑥
−𝜋
−𝑦
𝜋 𝜋
sin(𝑟𝑥)
න cos(𝑟𝑥)𝑑𝑥 = [ ]
−𝜋 𝑟 −𝜋
1 sin (𝑥) = −sin(−𝑥)
= [sin(𝑟𝜋) − sin(−𝑟𝜋)] 𝑦
𝑟 sin(−𝑥) = − ⁄ℎ = −sin (𝑥)
1
= ∙ 2sin(𝑟𝜋) = 0
𝑟
𝜋
න cos(𝑟𝑥)𝑑𝑥 = 0, ∀𝑟
−𝜋

𝜋
(1) 𝐼 = න cos(𝑟𝑥)cos(𝑠𝑥)𝑑𝑥 if 𝑟 = 𝑠 = 0
−𝜋

𝜋
𝐼 = න 𝑑𝑥 = [𝑥]𝜋−𝜋 = 2𝜋
−𝜋

cos(𝐴 + 𝐵) = cos(𝐴)cos(𝐵) − sin(𝐴)sin(𝐵)

Page | 61
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

cos(𝐴 − 𝐵) = cos(𝐴)cos(𝐵) + sin(𝐴)sin(𝐵)

cos(𝐴 + 𝐵) + cos(𝐴 − 𝐵) = 2 × cos(𝐴) × cos(𝐵) let 𝐴 = 𝑟𝑋

𝐵 = 𝑆𝑋
1
[cos(𝑟𝑥 + 𝑠𝑥) + cos(𝑟𝑥 − 𝑠𝑥)] = 2cos(𝑟𝑥)cos(𝑠𝑥)
2

1 𝜋
I = න (cos[(𝑟 + 𝑠)𝑥] + cos[(𝑟 − 𝑥)𝑥])𝑑𝑥
2 −𝜋
𝜋
1 sin[(𝑟 + 𝑠)𝑥] sin[(𝑟 − 𝑠)𝑥]
= [ + ]
2 𝑟+𝑠 𝑟−𝑠 −𝜋

1 sin(𝑛𝜋) = 0
= [(𝑟 − 𝑠)sin(𝑟 + 𝑠)𝑥 + (𝑟 + 𝑠)sin(𝑟 − 𝑠)𝑥]𝜋−𝜋
2(𝑟 2 − 𝑠2) ∀𝑛

1
= [(𝑟 − 𝑠)sin(𝑟 + 𝑠)𝜋 + (𝑟 + 𝑠)sin(𝑟 − 𝑠)𝜋]
2(𝑟 2 − 𝑠 2 )
−(−(𝑟 − 𝑠)𝑠𝑖𝑛(𝑟 + 𝑠)𝜋 − (𝑟 + 𝑠)𝑠𝑖𝑛(𝑟 − 𝑠)𝜋)] = 0, ∀𝑟 ≠ 𝑠

𝑟=𝑠>0 sin2(𝑟𝑥) = 1 − cos2 (𝑟𝑥)

1 𝜋 cos 2(𝑟𝑥) − sin2 (𝑟𝑥)) = cos (2𝑟𝑥)
𝐼 = න cos 2(𝑟𝑥)𝑑𝑥
2 −𝜋
2cos 2 (𝑟𝑥) − 1 = cos(2𝑟𝑥)
1
cos 2 (𝑟𝑥) = [1 + cos(2𝑟𝑥)]
2
1 𝜋
∴ 𝐼 = න 1 + cos(2𝑟𝑥)𝑑𝑥
2 −𝜋
1
= 𝜋 =𝜋
sin(2𝑟𝑥)
2 ൤𝑥 + 2𝑟 ൨ −𝜋

𝜋
(3) 𝐼 = න sin(rx) cos(𝑟𝑥)𝑑𝑥 = 0 ∀ r&s
−𝜋

sin(𝐴 + 𝐵) + sin(𝐴 − 𝐵) = 2sin(𝐴)cos(𝐵)

sin(𝐴 + 𝐵) = sin(𝐴)cos(𝐵) + cos(𝐴)sin(𝐵)

sin(𝐴 − 𝐵) = sin(𝐴)cos(𝐵) − cos(𝐴)sin(𝐵) for 𝐴 = 𝑟𝑥, 𝐵 = 𝑠𝑥,

1
[sin(𝑟𝑥 + 𝑠𝑥) + sin(𝑟𝑥 − 𝑠𝑥)] = sin(𝐴)cos(𝐵)
2

Page | 62
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

1 𝜋
𝐼 = න sin[(𝑟 + 𝑠)𝑥] + sin[(𝑟 − 𝑠)𝑥]𝑑𝑥
2 −𝜋
𝜋
1 cos[(𝑟 + 𝑠)𝑥] cos[(𝑟 − 𝑠)𝑥]
=− [ + ]
2 𝑟+𝑠 𝑟−𝑠 −𝜋
−1
= [(𝑟 − 𝑠)cos(𝑟 + 𝑠)𝑥 + (𝑟 + 𝑠)cos(𝑟 − 𝑠)𝑥]𝜋−𝜋
2(𝑟 2 − 𝑠 2 )
−1
= [(𝑟 − 𝑠)cos(𝑟𝜋 + 𝑠𝜋) + (𝑟 + 𝑠)cos(𝑟𝜋 − 𝑠𝜋)]𝜋−𝜋 = 0 (𝑄. 𝐸. 𝐷. )
2(𝑟 − 𝑠 2 )
2

Now, let’s create the Fourier series,

∞
𝑎0
𝑓(𝑥) = + ෍ 𝑎𝑟 cos(𝑟𝑥) + 𝑏𝑟 sin(𝑟𝑥)
2
𝑟=1

multiply the series by cos(𝑥) and integrate from 𝜋 to −𝜋

𝜋 ∞
𝑎0 𝜋 𝜋
න 𝑓(𝑥)cos(𝑥)𝑑𝑥 = න cos(𝑠𝑥)𝑑𝑥 + න cos(𝑠𝑥) {෍ 𝑎𝑟 cos(𝑟𝑥) + 𝑏𝑟 sin(𝑟𝑥)} 𝑑𝑥
−𝜋 2 −𝜋 −𝜋 𝑟=1

"uniform convergence"
𝜋 ∞ ∞
𝑎0 𝜋 𝜋 𝜋
න 𝑓(𝑥)cos(𝑥)𝑑𝑥 = න cos(𝑠𝑥)𝑑𝑥 + ෍ න 𝑎𝑟 cos(𝑟𝑥)sin(𝑟𝑥)𝑑𝑥 + ෍ න 𝑏𝑟 cos(𝑠𝑥)sin(𝑟𝑥)𝑑𝑥
−𝜋 2 −𝜋 −𝜋 −𝜋
𝑟=1 𝑟=1

when 𝑠 = 0
𝜋
𝑎0 1 𝜋
න 𝑓(𝑥)𝑑𝑥 = · 2𝜋 ⇒ 𝑎0 = න 𝑓(𝑥)𝑑𝑥
−𝜋 2 𝜋 −𝜋

Page | 63
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

Lecture 2

We had
𝜋 ∞ ∞
𝑎0 𝜋 𝜋 𝜋
න 𝑓(𝑥)cos(𝑥)𝑑𝑥 = න cos(𝑠𝑥)𝑑𝑥 + ෍ න 𝑎𝑟 cos(𝑟𝑥)sin(𝑟𝑥)𝑑𝑥 + ෍ න 𝑏𝑟 cos(𝑠𝑥)sin(𝑟𝑥)𝑑𝑥
−𝜋 2 −𝜋 −𝜋 −𝜋
𝑟=1 𝑟=1

when 𝑠 is non-zero, only non-zero term on the R. H. S. occurs in the first when 𝑠 = 𝑟

𝜋 𝜋
න 𝑓(𝑥)cos(𝑥)𝑑𝑥 = 𝑎𝑟 න cos(𝑟𝑥) cos(𝑟𝑥) 𝑑𝑥 = 𝜋𝑎𝑟
−𝜋 −𝜋
𝜋
1
𝑎𝑟 = න 𝑓(𝑥)cos(𝑟𝑥)𝑑𝑥
𝜋 −𝜋

Similarly, to obtain the Fourier coefficients 𝑏𝑟 , multiply

∞
𝑎0
𝑓(𝑥) = + ෍ 𝑎𝑟 cos(𝑟𝑥) + 𝑏𝑟 sin(𝑟𝑥)
2
𝑟=1

by 𝑠𝑖𝑛(𝑠𝑥) and integrate between −𝜋 and 𝜋. We easily find

𝜋 𝜋
න 𝑓(𝑥)sin(𝑟𝑥)𝑑𝑥 = 𝑏𝑟 න 𝑠𝑖𝑛(𝑟𝑥)sin(𝑟𝑥)𝑑𝑥 = 𝜋𝑏𝑟
−𝜋 −𝜋

1 𝜋
𝑏𝑟 = න 𝑓(𝑥)sin(𝑟𝑥)𝑑𝑥
𝜋 −𝜋

In conclusion, provided one function is periodic of period 2𝜋, we represent the function as
∞
𝑎0
𝑓(𝑥) = + ෍ 𝑎𝑟 cos(𝑟𝑥) + 𝑏𝑟 sin(𝑟𝑥)
2
𝑟=1

1 𝜋
𝑎0 = න 𝑓(𝑥)𝑑𝑥
𝜋 −𝜋

Page | 64
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

1 𝜋 𝑟 = 1,2,3 …
𝑎𝑟 = න 𝑓(𝑥)cos(𝑟𝑥)𝑑𝑥
𝜋 −𝜋
1 𝜋 𝑟 = 1,2,3 …
𝑏𝑟 = න 𝑓(𝑥)sin(𝑟𝑥)𝑑𝑥
𝜋 −𝜋

Lecture 3

Find the Fourier series representation of the non-periodic function 𝑓(𝑥) = 𝑥 2 on −𝜋 ≤ 𝑥 ≤ 𝜋.

First extend function 𝑓(𝑥) so 𝑓(𝑥 + 2𝑛𝜋) = 𝑓(𝑥)

We know

1 𝜋 𝑓(𝑥) = 𝑓(−𝑥) ⇒ it is even ⇒ 𝑏𝑟 = 0

𝑎0 = න 𝑓(𝑥)𝑑𝑥
𝜋 −𝜋
𝐿
1 𝜋 𝐿
= න x 2 𝑑𝑥 න 𝑓(𝑥)𝑑𝑥 = 2 න 𝑓(𝑥)𝑑𝑥
𝜋 0
−𝐿
1 3𝜋 0
= [𝑥 ]−𝜋
3𝜋 𝐿 𝐿
2𝜋 2 න 𝑓(𝑥)𝑑𝑥 = 2 න 𝑓(𝑥)𝑑𝑥
= ⏟−𝐿 0
3
if 𝑓(𝑥)is even

Page | 65
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

1 𝜋
𝑎𝑟 = න 𝑓(𝑥) cos(𝑟𝑥) 𝑑𝑥
𝜋 −𝜋 𝑟 = 1,2,3 …

𝑑𝑣 𝑑𝑢
න𝑢 𝑑𝑡 = 𝑢𝑣 − න 𝑣 𝑑𝑡
𝑑𝑡 𝑑𝑡

12 𝜋 2
𝑎𝑟 = න 𝑥 cos(𝑟𝑥) 𝑑𝑥 𝑟 = 1,2,3 …
𝜋 0

Let 𝑢 = 𝑥 2 ,
𝑑𝑢 𝑑𝑣 sin(𝑟𝑥)
= 2𝑥, = cos(𝑟𝑥), 𝑣=
𝑑𝑥 𝑑𝑥 𝑟
𝜋
2 𝑥 2 · sin(𝑟𝑥) 2 𝜋
𝑎𝑟 = ቊ[ ] − න 𝑥 · sin(𝑟𝑥)𝑑𝑥 ቋ
𝜋 𝑟 0
𝑟 0
−4 𝜋
= න 𝑥 · sin(𝑟𝑥)𝑑𝑥
𝑟𝜋 0

𝜋
2 𝑥 2 · sin(𝑟𝑥) 2 𝜋
𝑎𝑟 = ቊ[ ] − න 𝑥 · sin(𝑟𝑥)𝑑𝑥 ቋ
𝜋 𝑟 0
𝑟 0
Let 𝑢 = 𝑥,
−4 𝜋
= න 𝑥 · sin(𝑟𝑥)𝑑𝑥 𝑑𝑢 𝑑𝑣 −cos(𝑟𝑥)
𝑟𝜋 0 = 1, = sin(𝑟𝑥), 𝑣 =
𝜋 𝑑𝑥 𝑑𝑥 𝑟
𝜋
−4 −𝑥 · cos(𝑟𝑥) cos(𝑟𝑥)
= ቊ[ ] +න 𝑑𝑥 ቋ
𝑟𝜋 𝑟 0 0 𝑟
𝜋
−4 −𝜋 sin(𝑟𝑥)
= [ cos(𝑟𝜋)] + [ ]
𝑟𝜋 𝑟 𝑟2 0
cos(𝑟𝜋) = (−1)𝑟
4 2 × 𝜋2
𝑎𝑟 = 2 (−1)𝑟 𝑎0 =
𝑟 3
4
𝑏𝑟 = 0, ∀𝑟 𝑎1 = 2 (−1)1 = −4
1
4 1
𝑎2 = 2 (−1)2 = 4 ( 2 )
2 2

𝜋2 cos(𝑥) cos(2𝑥) cos(3𝑥) cos(4𝑥)

𝑓(𝑥) = − 4ቆ − + − +⋯ቇ
2 1 22 32 42

∞
𝜋2 (−1)𝑟
𝑓(𝑥) = +4෍ cos(𝑟𝑥)
3 𝑟2
𝑟=1

First partial sum is

𝜋2
𝑆1 =
3

Page | 66
Partial Differential UCL Engineering
Teaching material
Equations Dr Vasos Pavlika

𝜋 2 4cos(𝑥)
𝑠2 = −
3 12

𝜋 2 4cos(𝑥) 4cos(2𝑥)
𝑠3 = − +
3 12 22

Page | 67