Partial Differential Equations

Heat diffusion equation

1
Partial Differential Equations

A partial differential equation (PDE) is an

equation that involves an unknown function
and its partial derivatives.

Example :
∂ 2 u ( x, t ) ∂ u ( x, t )
=
∂x 2 ∂t
PDE involves two or more independen t variables
(in the example x and t are independen t variables)

2
 Notation

∂ 2 u ( x, t )
u xx =
∂x 2
∂ 2 u ( x, t )
u xt =
∂x ∂t
Order of the PDE = order of the highest order derivative.

3
Linear PDE
Classification
A PDE is linear if it is linear in the unknown
function and its derivatives
Example of linear PDE :
2 u xx + 1 u xt + 3 utt + 4 u x + cos(2t ) = 0
2 u xx − 3 ut + 4 u x = 0
Examples of Nonlinear PDE
2 u xx + (u xt )2 + 3 utt = 0
u xx + 2 u xt + 3 ut = 0
2 u xx + 2 u xt ut + 3 ut = 0
4
 Representing the Solution of a PDE
(Two Independent Variables)
• Three main ways to represent the solution

T ( x1 , t1 ) T=5.2

t1 T=3.5

x1
Different curves are Three dimensional The axis represent
used for different plot of the function the independent
values of one of the T(x,t) variables. The value
independent
of the function is
variable
displayed at grid
points
5
 The Heat Equation

Incorporation of the constitutive equation into the energy

equation above yields:

Dividing both sides by rCp and introducing the thermal

diffusivity of the material given by
For constant thermal properties and no heat generation.

This is often called the heat equation.

General conduction equation based on Cartesian
Coordinates
 ∂T
ρC p = ∇.k∇T + g ( x, t )
∂t
For an isotropic and homogeneous material:

∂T
ρC p 2
= k∇ T + g ( x, t )
∂t
∂T ∂ T ∂ T ∂ T 2 2 2
ρC p = k  2 + 2 + 2  + g ( x, y , z : t )
∂t  ∂x ∂y ∂z 
 General conduction equation based on Polar
Cylindrical Coordinates

∂T ∂ T 1 ∂ T ∂ T 
2 2 2
ρC p = k 2 + + 2  + g ( x, y , z : t )
∂t  ∂r r ∂θ 2
∂z 
Thermal Conductivity of Brick Masonry Walls
 Thermally Heterogeneous Materials

k = k ( x, y , z )
∂T
ρC p = ∇.k∇T + g ( x, t )
∂t

 ∂T  ∂ k ∂T   ∂T 
∂ k    ∂ k 
∂T  ∂x   ∂y   ∂z 
ρC p = + + + g ( x, y , z , t )
∂t ∂x ∂y ∂z
 ∂T ∂k ∂T ∂ 2T ∂k ∂T ∂ 2T ∂k ∂T ∂ 2T
ρC p = +k 2 + +k 2 + + k 2 + g ( x, y , z , t )
∂t ∂x ∂x ∂x ∂y ∂y ∂y ∂z ∂z ∂z
Satellite Imaging : Remote Sensing
Ultra-sound Imaging of Brain
 Steady-State One-Dimensional Conduction

• For one-dimensional steady state

conduction with no energy generation,
the heat equation reduces to :

∂T ∂k ∂T ∂ 2T
ρC p = + k 2 + g ( x, y , z , t )
∂t ∂x ∂x ∂x

• Assume a homogeneous medium with invariant thermal conductivity ( k =

constant) :
∂T ∂ 2T
ρC p = k 2 + g ( x, y , z , t )
∂t ∂x
One dimensional Transient conduction with heat generation.
Steady-State One-Dimensional Conduction
• For one-dimensional heat conduction
in a variable area geometry.
• We can devise a basic description of
the process.
• The first law in control volume form
(steady flow energy equation) with no
shaft work and no mass flow reduces
to the statement that ΣQ = 0 for all
surfaces.
• From Fourier law of conduction, the
heat transfer rate in at the left (at x) is:
 Taylor’s Theory of Continuum

Q ( x + dx ) = Q ( x ) +
dQ (x )
dx +
d 2 
Q ( x ) (dx )2
+
d 3 
Q ( x ) (dx )3
+ ......
2 3
dx dx 2! dx 3!

n d Q ( x ) (dx )
n  n n
Q( x + dx ) = Q( x ) + ∑ (− 1)
 
n
i =1 dx n!
• For a function converging & well behaving…

d  (x )
Q
Q ( x + dx ) = Q ( x ) + dx
dx x
• For a pure steady state conduction:

Q ( x + dx ) − Q ( x ) = 0
 d  (x )
Q
Q ( x ) + dx − Q ( x ) = 0
dx x

dQ ( x )
=0
dx x
Substitute Fourier’s law of conduction:

 dT 
d  − kA 
 dx 
=0
dx
x
  dT 
d  kA 
 dx 
=0
dx
If k is constant (i.e. if the material is homogeneous and properties of the
medium are independent of temperature), this reduces to

d 2T  dA  dT 
A 2 +   =0
dx  dx  dx 

Pure radial conduction through

A Sphere.

d 2T  dA  dT 
A 2 +   =0
dr  dr  dr 
Surface area of a sphere at r

dA
A = 4πr 2
& = 8πr
dr

d 2T  1  dT 
2
+   =0
dr  2r  dr 
Heat transfer through a plane slab

d 2T dT
A 2 =0⇒ = C1 ⇒ T = C1 x + C2
dx dx
Isothermal Wall Surfaces
 Wall Surfaces with Convection

d 2T dT
A 2 =0⇒ = C1 ⇒ T = C1 x + C2
dx dx

Boundary conditions:

dT
−k = h1 (T (0) − T∞1 )
dx x =0

dT
−k = h2 (T ( L) − T∞ 2 )
dx x=L
 Wall with isothermal Surface and Convection Wall

d 2T dT
A 2 =0⇒ = C1 ⇒ T = C1 x + C2
dx dx

Boundary conditions:

T ( x = 0) = T1

dT
−k = h2 (T ( L) − T∞ 2 )
dx x=L
 Electrical Circuit Theory of Heat Transfer
• Thermal Resistance
• A resistance can be defined as the ratio of a driving potential
to a corresponding transfer rate.

∆V
R=
i
Analogy:
Electrical resistance is to conduction of electricity as thermal
resistance is to conduction of heat.
The analog of Q is current, and the analog of the temperature
difference, T1 - T2, is voltage difference.
From this perspective the slab is a pure resistance to heat transfer
and we can define
 ∆T
Q=
Rth
 The composite Wall
• The concept of a thermal
resistance circuit allows ready
analysis of problems such as a
composite slab (composite planar
heat transfer surface).
• In the composite slab, the heat
flux is constant with x.
• The resistances are in series and
sum to R = R1 + R2.
• If TL is the temperature at the left,
and TR is the temperature at the
right, the heat transfer rate is
given by
 Wall Surfaces with Convection

d 2T dT
A 2 =0⇒ = C1 ⇒ T = C1 x + C2
dx dx

Boundary conditions:

dT
−k = h1 (T (0) − T∞1 ) T∞1
dx x =0
T∞2

dT
−k = h2 (T ( L) − T∞ 2 ) Rconv,1 Rcond Rconv,2
dx x=L
Heat Equation Different curve is
used for each value
of t

ice ice Temperature Temperature at

x different x at t=0
Thin metal rod insulated
everywhere except at the
edges. At t =0 the rod is
placed in ice

∂ 2 T ( x, t ) ∂ T ( x, t ) Position x
2
− =0
∂x ∂t Temperature at
T (0, t ) = T (1, t ) = 0 different x at t=h

T ( x,0) = sin(π x)

34
Examples of PDEs

PDEs are used to model many systems in many different fields of

science and engineering.

Important Examples:
• Laplace Equation
• Heat Equation
• Wave Equation

35
Laplace Equation

∂ 2u ( x , y , z ) ∂ 2u ( x , y , z ) ∂ 2u ( x , y , z )
2
+ 2
+ 2
=0
∂x ∂y ∂z
Used to describe the steady state distribution of
heat in a body.
Also used to describe the steady state
distribution of electrical charge in a body.

36
Heat Equation

∂ u ( x, y , z, t )  ∂ 2u ∂ 2u ∂ 2u 
= α  2 + 2 + 2 
∂t 
 ∂x ∂y ∂z 

The function u(x,y,z,t) is used to represent

the temperature at time t in a physical body
at a point with coordinates (x,y,z)
α is the thermal diffusivity. It is sufficient to
consider the case α = 1.

37
Simpler Heat Equation

2
∂ T ( x, t ) ∂ T ( x, t )
= x
∂t ∂x 2

T(x,t) is used to represent the temperature

at time t at the point x of the thin rod.

38
Wave Equation

∂ 2u ( x , y , z , t ) 2  ∂ 2
u ∂ 2
u ∂ 2
u 
= c  2 + 2 + 2 
∂t 2 
 ∂x ∂y ∂z 

The function u(x,y,z,t) is used to represent the

displacement at time t of a particle whose
position at rest is (x,y,z) .
The constant c represents the propagation
speed of the wave.

39
Classification of PDEs

Linear Second order PDEs are important sets of equations that are
used to model many systems in many different fields of science and
engineering.

Classification is important because:

• Each category relates to specific engineering problems.
• Different approaches are used to solve these categories.

40
Linear Second Order PDEs
Classification
A second order linear PDE (2 - independent variables)
A u xx + B u xy + C u yy + D = 0,
A, B, and C are functions of x and y
D is a function of x, y , u, u x , and u y
is classified based on (B 2 − 4 AC) as follows :
B 2 − 4 AC < 0 Elliptic
B 2 − 4 AC = 0 Parabolic
B 2 − 4 AC > 0 Hyperbolic

41
Linear Second Order PDE
Examples (Classification)
∂ 2u ( x , y ) ∂ 2u ( x , y )
Laplace Equation 2
+ 2
=0
∂x ∂y
A = 1, B = 0, C = 1 ⇒ B 2 − 4 AC < 0
⇒ Laplace Equation is Elliptic
One possible solution : u( x, y ) = e x sin y
u x = e x sin y , u xx = e x sin y
u y = e x cos y , u yy = −e x sin y
u xx + u yy = 0

42
Linear Second Order PDE
Examples (Classification)
∂ 2u ( x , t ) ∂ u ( x , t )
Heat Equation α 2
− =0
∂x ∂t
A = α , B = 0, C = 0 ⇒ B 2 − 4 AC = 0
⇒ Heat Equation is Parabolic
______________________________________
2 2
2 ∂ u ( x , t ) ∂ u ( x, t )
Wave Equation c 2
− 2
=0
∂x ∂t
A = c 2 > 0, B = 0, C = −1 ⇒ B 2 − 4 AC > 0
⇒ Wave Equation is Hyperbolic
43
 Boundary Conditions for PDEs
• To uniquely specify a solution to the PDE, a set of boundary conditions
are needed.
• Both regular and irregular boundaries are possible.

t
∂ 2u( x, t ) ∂u( x, t )
Heat Equation : α 2
− =0
∂x ∂t region of
interest
u(0, t ) = 0
u(1, t ) = 0
x
u( x,0) = sin(π x ) 1

CIS301_Topic9
44
The Solution Methods for PDEs

• Analytic solutions are possible for simple and special (idealized) cases
only.

• To make use of the nature of the equations, different methods are

used to solve different classes of PDEs.

• The methods discussed here are based on the finite difference

technique.

45
 Parabolic Equations
 Parabolic Equations
 Heat Conduction Equation
 Explicit Method
 Implicit Method
 Cranks Nicolson Method

46
 Parabolic Equations
A second order linear PDE (2 - independent variables x , y )
A u xx + B u xy + C u yy + D = 0,
A, B, and C are functions of x and y
D is a function of x, y, u, u x , and u y

is parabolic if B 2 − 4 AC = 0

47
 Parabolic Problems
∂ T ( x, t ) ∂ 2 T ( x, t )
Heat Equation : =
∂t ∂x 2
T (0, t ) = T (1, t ) = 0
T ( x,0) = sin(π x ) ice ice

* Parabolic problem ( B 2 − 4 AC = 0)
* Boundary conditions are needed to uniquely specify a solution.

48
 Finite Difference Methods
• Divide the interval x into sub-intervals, each of width h
• Divide the interval t into sub-intervals, each of width k
• A grid of points is used for
the finite difference solution
• Ti,j represents T(xi, tj) t
• Replace the derivates by
finite-difference formulas

49
Finite Difference Methods
Replace the derivatives by finite difference formulas
∂ 2T
Central Difference Formula for 2 :
∂x
∂ 2T ( x, t ) Ti −1, j − 2Ti , j + Ti +1, j Ti −1, j − 2Ti , j + Ti +1, j
2
≈ 2
=
∂x ( ∆x ) h2
∂T
Forward Difference Formula for :
∂t
∂T ( x, t ) Ti , j +1 − Ti , j Ti , j +1 − Ti , j
≈ =
∂t ∆t k

50
Solution of the Heat Equation

• Two solutions to the Parabolic Equation

(Heat Equation) will be presented:
1. Explicit Method:
Simple, Stability Problems.
2. Crank-Nicolson Method:
Involves the solution of a Tridiagonal system
of equations, Stable.

51
Explicit Method
∂T ( x, t ) ∂ 2T ( x, t )
=
∂t ∂x 2
T ( x, t + k ) − T ( x, t ) T ( x − h, t ) − 2T ( x, t ) + T ( x + h, t )
=
k h2
k
T ( x, t + k ) − T ( x, t ) = 2 (T ( x − h, t ) − 2T ( x, t ) + T ( x + h, t ) )
h
k
Define λ = 2
h
T ( x, t + k ) = λ T ( x − h, t ) + (1 − 2 λ ) T ( x, t ) + λ T ( x + h, t )

52
 Explicit Method
How Do We Compute?
T ( x, t + k ) = λ T ( x − h, t ) + (1 − 2 λ ) T ( x, t ) + λ T ( x + h, t )
means

T(x,t+k)

T(x-h,t) T(x,t) T(x+h,t)

53
 Convergence and Stability
T ( x, t + k ) can be computed directly using :
T ( x, t + k ) = λ T ( x − h, t ) + (1 − 2 λ ) T ( x, t ) + λ T ( x + h, t )

Can be unstable (errors are magnified )

1 h2
To guarantee stability, (1 − 2 λ ) ≥ 0 ⇒ λ ≤ ⇒ k ≤
2 2
This means that k is much smaller than h
This makes it slow.

54
 Convergence and Stability of the Solution
• Convergence
The solutions converge means that the solution obtained using the finite
difference method approaches the true solution as the steps
approach zero.

• Stability: ∆x and ∆t
An algorithm is stable if the errors at each stage of the computation are
not magnified as the computation progresses.

55
Example 1: Heat Equation
Solve the PDE :
∂ 2u(x,t) ∂u(x,t)
2
− =0
∂ x ∂t
u(0, t ) = u(1, t ) = 0
ice ice
u( x,0) = sin(π x )
x

Use h = 0.25, k = 0.25 to find u( x, t ) for x ∈ [0,1], t ∈ [0,1]

k
λ = 2 =4
h

56
 Example 1
∂ 2 u( x, t ) ∂ u ( x, t )
2
− =0
∂x ∂t
u ( x − h, t ) − 2u ( x, t ) + u ( x + h, t ) u ( x, t + k ) − u ( x, t )
2
− =0
h k
16(u ( x − h, t ) − 2u ( x, t ) + u ( x + h, t ) ) − 4(u ( x, t + k ) − u( x, t ) ) = 0

u ( x , t + k ) = 4 u ( x − h, t ) − 7 u ( x , t ) + 4 u ( x + h, t )

57
Example 1
u ( x , t + k ) = 4 u ( x − h, t ) − 7 u ( x , t ) + 4 u ( x + h, t )

t=1.0 0 0

t=0.75 0 0

t=0.5 0 0
t=0.25 0 0

t=0 0 0
Sin(0.25π) Sin(0. 5π) Sin(0.75π)

x=0.0 x=0.25 x=0.5 x=0.75 x=1.0

58
Example 1
u(0.25,0.25) = 4 u(0,0) − 7 u(0.25,0) + 4 u(0.5,0)
= 0 − 7 sin(π / 4) + 4 sin(π / 2) = −0.9497

t=1.0 0 0

t=0.75 0 0

t=0.5 0 0
t=0.25 0 0

t=0 0 0
Sin(0.25π) Sin(0. 5π) Sin(0.75π)

x=0.0 x=0.25 x=0.5 x=0.75 x=1.0

59
Example 1
u(0.5,0.25) = 4 u(0.25,0) − 7 u(0.5,0) + 4 u(0.75,0)
= 4 sin(π / 4) − 7 sin(π / 2) + 4 sin(3π / 4) = −0.1716

t=1.0 0 0

t=0.75 0 0

t=0.5 0 0
t=0.25 0 0

t=0 0 0
Sin(0.25π) Sin(0. 5π) Sin(0.75π)

x=0.0 x=0.25 x=0.5 x=0.75 x=1.0

60
Remarks on Example 1
The obtained results are probably not accurate
because : 1 − 2λ = −7

For accurate results : 1 − 2λ ≥ 0

h 2 (0.25) 2
One needs to select k ≤ = = 0.03125
2 2
k
For example, choose k = 0.025, then λ = 2 = 0.4
h

61
Example 1 – cont’d
u( x, t + k ) = 0.4 u( x − h, t ) + 0.2 u( x, t ) + 0.4 u( x + h, t )

t=0.10 0 0

t=0.075 0 0

t=0.05 0 0
t=0.025 0 0

t=0 0 0
Sin(0.25π) Sin(0. 5π) Sin(0.75π)

x=0.0 x=0.25 x=0.5 x=0.75 x=1.0

62
Example 1 – cont’d
u (0.25,0.025) = 0.4 u (0,0) + 0.2 u (0.25,0) + 0.4 u(0.5,0)
= 0 + 0.2 sin(π / 4) + 0.4 sin(π / 2) = 0.5414

t=0.10 0 0

t=0.075 0 0

t=0.05 0 0
t=0.025 0 0

t=0 0 0
Sin(0.25π) Sin(0. 5π) Sin(0.75π)

x=0.0 x=0.25 x=0.5 x=0.75 x=1.0

63
Example 1 – cont’d
u (0.5,0.025) = 0.4 u (0.25,0) + 0.2 u (0.5,0) + 0.4 u(0.75,0)
= 0.4 sin(π / 4) + 0.2 sin(π / 2) + 0.4 sin(3π / 4) = 0.7657

t=0.10 0 0

t=0.075 0 0

t=0.05 0 0
t=0.025 0 0

t=0 0 0
Sin(0.25π) Sin(0. 5π) Sin(0.75π)

x=0.0 x=0.25 x=0.5 x=0.75 x=1.0

64
 Crank-Nicolson Method
The method involves solving a Tridiagonal system of linear equations.
The method is stable (No magnification of error).
→ We can use larger h, k (compared to the Explicit Method).

65
Crank-Nicolson Method
Based on the finite difference method
1. Divide the interval x into subintervals of width h
2. Divide the interval t into subintervals of width k
3. Replace the first and second partial derivatives with their
backward and central difference formulas respectively :
∂u ( x, t ) u ( x, t ) − u ( x, t − k )
≈
∂t k
∂ 2 u ( x , t ) u ( x − h , t ) − 2u ( x , t ) + u ( x + h , t )
2
=
∂x h2
66
Crank-Nicolson Method
∂ 2 u ( x, t ) ∂ u ( x, t )
Heat Equation : 2
= becomes
∂x ∂t

u ( x − h , t ) − 2u ( x , t ) + u ( x + h , t ) u ( x , t ) − u ( x , t − k )
2
=
h k
k
2
(u( x − h, t ) − 2u( x, t ) + u( x + h, t ) ) = u( x, t ) − u( x, t − k )
h
k k k
− 2 u ( x − h, t ) + (1 + 2 2 ) u( x, t ) − 2 u( x + h, t ) = u( x, t − k )
h h h

67
 Crank-Nicolson Method
k
Define λ = 2 then Heat equation becomes :
h
− λ u( x − h, t ) + (1 + 2λ ) u( x, t ) − λ u( x + h, t ) = u( x, t − k )

u(x-h,t) u(x,t) u(x+h,t)

u(x,t - k)

68
Crank-Nicolson Method
The equation :
− λ u( x − h, t ) + (1 + 2λ ) u( x, t ) − λ u( x + h, t ) = u( x, t − k )
can be rewritten as :
− λ ui −1, j + (1 + 2λ ) ui , j − λ ui +1, j = ui , j −1
and can be expanded as a system of equations (fix j = 1) :
− λ u0,1 + (1 + 2λ ) u1,1 − λ u2,1 = u1,0
− λ u1,1 + (1 + 2λ ) u2,1 − λ u3,1 = u2,0
− λ u2,1 + (1 + 2λ ) u3,1 − λ u4,1 = u3,0
− λ u3,1 + (1 + 2λ ) u4,1 − λ u5,1 = u4,0

69
Crank-Nicolson Method
− λ u( x − h, t ) + (1 + 2λ ) u( x, t ) − λ u( x + h, t ) = u( x, t − k )
can be expressed as a Tridiagonal system of equations :
1 + 2λ −λ   u1,1   u1,0 + λ u0,1 
 − λ 1 + 2λ  u   u 
−λ 
   2,1  =  2,0
 − λ 1 + 2λ − λ  u3,1   u3,0 
  u   u + λ u 
 − λ 1 + 2λ   4,1   4,0 5,1 

where u1,0 , u2,0 , u3,0 , and u4,0 are the initial temperature values
at x = x0 + h, x0 + 2h, x0 + 3h, and x0 + 4h
u0,1 and u5,1 are the boundary values at x = x0 and x0 + 5h

70
Crank-Nicolson Method

The solution of the tridiagonal system produces :

The temperature values u1,1, u2,1, u3,1, and u4,1 at t = t0 + k
To compute the temperature values at t = t0 + 2k
Solve a second tridiagonal system of equations ( j = 2)
1 + 2λ −λ   u1,2   u1,1 + λ u0,2 
 − λ 1 + 2λ  u   u2,1 
−λ  2, 2  =  
 
 − λ 1 + 2λ −λ  u3,2   u3,1 
     
 − λ 1 + 2λ  u4,2  u4,1 + λ u5,2 
To compute u1,2 , u2,2 , u3,2 , and u4,2
Repeat the above step to compute temperature values at t0 + 3k , etc.

71
 Example 2
Solve the PDE :
∂ 2u( x, t ) ∂u( x, t )
2
− =0
∂ x ∂t
u(0, t ) = u(1, t ) = 0
u( x,0) = sin(π x )

Solve using Crank - Nicolson method

Use h = 0.25, k = 0.25 to find u( x, t ) for x ∈ [0,1], t ∈ [0,1]

72
Example 2
Crank-Nicolson Method
∂ 2 u ( x, t ) ∂ u ( x, t )
2
− =0
∂x ∂t
u ( x − h , t ) − 2u ( x , t ) + u ( x + h , t ) u ( x , t ) − u ( x , t − k )
2
=
h k
16(u( x − h, t ) − 2u( x, t ) + u( x + h, t ) ) − 4(u( x, t ) − u( x, t − k ) ) = 0
k
Define λ = 2 = 4
h
− 4 u ( x − h, t ) + 9 u ( x , t ) − 4 u ( x + h, t ) = u ( x , t − k )
− 4 ui −1, j + 9 ui , j − 4 ui +1, j = ui , j −1

73
Example 2
− 4u0,1 + 9u1,1 − 4u2,1 = u1,0 ⇒ 9u1,1 − 4u2,1 = sin(π / 4)
− 4u1,1 + 9u2,1 − 4u3,1 = u2,0 ⇒ −4u1,1 + 9u2,1 − 4u3,1 = sin(π / 2)
− 4u2,1 + 9u3,1 − 4u4,1 = u3,0 ⇒ − 4u2,1 + 9u3,1 = sin(3π / 4)
u1,4 u2,4 u3,4
t4=1.0 0 0
u1,3 u2,3 u3,3
t3=0.75 0 0
u1,2 u2,2 u3,2
t2=0.5 0 0
u1,1 u2,1 u3,1
t1=0.25 0 0

t0=0 0 0
Sin(0.25π) Sin(0. 5π) Sin(0.75π)
x0=0.0 x1=0.25 x2=0.5 x3=0.75 x4=1.0
74
Example 2
Solution of Row 1 at t1=0.25 sec
The Solution of the PDE at t1 = 0.25 sec is the solution
of the following tridiagonal system of equations :
 9 −4   u1,1  sin(0.25π )
 − 4 9 − 4 u  =  sin(0.5π ) 
   2,1   
 − 4 9  u3,1  sin(0.75π )
 u1,1   0.21151
  
⇒ u2,1  = 0.29912
 
u3,1   0.21151

75
Example 2:
Second Row at t2=0.5 sec
− 4u0,2 + 9u1,2 − 4u2,2 = u1,1 ⇒ 9u1,2 − 4u2,2 = 0.21151
− 4u1,2 + 9u2,2 − 4u3,2 = u2,1 ⇒ −4u1,2 + 9u2,2 − 4u3,2 = 0.29912
− 4u2,2 + 9u3,2 − 4u4,2 = u3,1 ⇒ − 4u2,2 + 9u3,2 = 0.21151
u1,4 u2,4 u3,4
t4=1.0 0 0
u1,3 u2,3 u3,3
t3=0.75 0 0
u1,2 u2,2 u3,2
t2=0.5 0 0
u1,1 u2,1 u3,1
t1=0.25 0 0

t0=0 0 0
Sin(0.25π) Sin(0. 5π) Sin(0.75π)
x0=0.0 x1=0.25 x2=0.5 x3=0.75 x4=1.0
76
Example 2
Solution of Row 2 at t2=0.5 sec
The Solution of the PDE at t2 = 0.5 sec is the solution
of the following tridiagonal system of equations :
 9 −4   u1,2   u1,1   0.21151
− 4 9 − 4 u  = u  = 0.29912
   2,2   2,1   
 − 4 9  u3,2  u3,1   0.21151
 u1,2  0.063267
  
⇒ u2,2  = 0.089473
 
u3,2  0.063267

77
Example 2
Solution of Row 3 at t3=0.75 sec
The Solution of the PDE at t3 = 0.75 sec is the solution
of the following tridiagonal system of equations :
 9 −4   u1,3   u1,2  0.063267
 − 4 9 − 4 u  = u  =  0.089473
   2,3   2, 2   
 − 4 9  u3,3  u3,2  0.063267
 u1,3  0.018924
 
⇒ u2,3  =  0.026763
 
u3,3  0.018924

78
Example 2
Solution of Row 4 at t4=1 sec
The Solution of the PDE at t4 = 1 sec is the solution
of the following tridiagonal system of equations :
 9 −4   u1,4   u1,3  0.018924
 − 4 9 − 4 u  = u  =  0.026763
   2, 4   2,3   
 − 4 9  u3,4  u3,3  0.018924
 u1,4  0.0056606
 
⇒ u2,4  = 0.0080053
 
u3,4  0.0056606

79
 Remarks
The Explicit Method:
• One needs to select small k to ensure stability.
• Computation per point is very simple but many
points are needed.

Cranks Nicolson:
• Requires the solution of a Tridiagonal system.
• Stable (Larger k can be used).

80
 Elliptic Equations
 Elliptic Equations
 Laplace Equation
 Solution

81
 Elliptic Equations

A second order linear PDE (2 - independent variables x , y )

A u xx + B u xy + C u yy + D = 0,
A, B, and C are functions of x and y
D is a function of x, y , u, u x , and u y

is Elliptic if B 2 − 4 AC < 0

CISE301_Topic9
82
 Laplace Equation
Laplace equation appears in several engineering problems such as:
• Studying the steady state distribution of heat in a body.
• Studying the steady state distribution of electrical charge in a body.

2 2
∂ T ( x, y ) ∂ T ( x, y )
2
+ 2
= f ( x, y )
∂x ∂y
T : steady state temperature at point (x, y)
f ( x, y ) : heat source (or heat sink)
83
 Laplace Equation

∂ 2 T ( x, y ) ∂ 2 T ( x, y )
2
+ 2
= f ( x, y )
∂x ∂y
A = 1, B = 0, C = 1
B 2 − 4 AC = −4 < 0 Elliptic

• Temperature is a function of the position (x and y)

• When no heat source is available f(x,y)=0

84
 Solution Technique

• A grid 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.

∂ 2 T ( x, y ) Ti +1, j − 2Ti , j + Ti −1, j ∂ 2 T ( x, y ) Ti , j +1 − 2Ti , j + Ti , j −1

≈ , ≈
∂x 2
(∆x ) 2
∂y 2
(∆y )2
85
Solution Technique
∂ 2 T ( x, y ) Ti +1, j − 2Ti , j + Ti −1, j
= ,
∂x 2
(∆x ) 2

∂ 2 T ( x, y ) Ti , j +1 − 2Ti , j + Ti , j −1
=
∂y 2
(∆y ) 2

∂ 2 T ( x, y ) ∂ 2 T ( x, y )
⇒ 2
+ 2
=0
∂x ∂y
is approximated by :
Ti +1, j − 2Ti , j + Ti −1, j Ti , j +1 − 2Ti , j + Ti , j −1
+ =0
(∆x ) 2
(∆y ) 2

86
 Solution Technique

Ti +1, j − 2Ti , j + Ti −1, j Ti , j +1 − 2Ti , j + Ti , j −1

+ =0
(∆x ) 2
(∆y ) 2

( Laplacian Difference Equation)

Assume : ∆x = ∆y = h
⇒ Ti +1, j + Ti −1, j + Ti , j +1 + Ti , j −1 − 4Ti , j = 0

87
Solution Technique

Ti , j +1

Ti −1, j Ti , j Ti +1, j

Ti , j −1

Ti +1, j + Ti −1, j + Ti , j +1 + Ti , j −1 − 4Ti , j = 0

88
Example
It is required to determine the steady state temperature at all points of a
heated sheet of metal. The edges of the sheet are kept at a constant
temperature: 100, 50, 0, and 75 degrees.

100

75 50

The sheet is divided

to 5X5 grids.

89
 Known

Example To be determined

T1, 4 = 100 T2, 4 = 100 T3, 4 = 100

T1,3 T2,3
T0,3 = 75 T3,3
T4,3 = 50
T1, 2 T2, 2 T3, 2
T0, 2 = 75 T4, 2 = 50

T1,1 T2,1 T3,1

T0,1 = 75 T4,1 = 50

T1, 0 = 0 T2, 0 = 0 T3, 0 = 0

90
 Known

First Equation To be determined

T1, 4 = 100 T2, 4 = 100

T1,3 T2,3
T0,3 = 75

T1, 2 T2, 2
T0, 2 = 75

T0,3 + T1, 4 + T1, 2 + T2,3 − 4T1,3 = 0

75 + 100 + T1, 2 + T2,3 − 4T1,3 = 0

91
 Known

Another Equation To be determined

T1, 4 = 100 T2, 4 = 100 T3, 4 = 100

T1,3 T2,3 T3,3

T1, 2 T2, 2 T3, 2

T1,3 + T2, 4 + T3,3 + T2, 2 − 4T2,3 = 0

T1,3 + 100 + T3,3 + T2, 2 − 4T2,3 = 0

92
Solution
The Rest of the Equations
 4 −1 0 −1   T1,1   75 
    
−1 4 −1 0 −1   T2,1   0 
 0 −1 4 0 0 −1   T   50 
   3,1   
−1 0 0 4 −1 0 −1   T1,2   75 
 −1 0 −1 4 −1 0 −1  T  =  0 
   2, 2   
 −1 0 −1 4 0 0 − 1  T3,2   50 
 − 1 0 0 4 − 1 0   T   175 
   1,3   
 − 1 0 − 1 4 − 1  T2,3  100 
    
 − 1 0 − 1 4 T
  3,3   150 

93
Thank You

94