Numerical methods for solving
the heat equation, the wave
equation and Laplace’s equation
(Finite difference methods)
Mona Rahmani
January 2019
Numerical methods are important tools to simulate
different physical phenomena.
Numerical simulation of a rotor
Numerical simulation of flow around a
Courtesy of NASA’s Ames Research Centre Porsche 956
Courtesy of NASA’s Ames Research Centre
Finite Difference Methods
• Idea: Use finite difference quotients to approximate the
derivatives in the PDE.
• Remember that the definition of the derivative of a function f(x) is:
• If is sufficiently small:
First derivative approximations
Now, consider the Taylor series expansion of
From which we can find
The truncation error
A forward difference approximation
With first order accuracy
First derivative approximations
How can we improve the accuracy of approximation of ?
Subtract the second line from the first line and solve for
The truncation error
A central difference approximation
With second order accuracy
Second derivative approximations
Again consider the two Taylor series expansions:
This time add the two lines together and solve for
The truncation error
A central difference approximation
With second order accuracy
Solving the 1D heat equation
Consider the initial-boundary value problem:
Boundary conditions (B. C.’s):
Initial condition (I. C.):
Step 1- Define a discretization in space and time: t
time step
k+1,
time step
k,
x0 = 0 xN = 1.0
x
Solving the 1D heat equation
Step 2 - Discretize the PDE. Use a forward difference scheme for the
time derivative and a central difference scheme for the space
derivative:
Or:
Which in the index notation is:
Solving the 1D heat equation
Step 3 - Write the discrete equations for all nodes in a matrix
format and solve the system:
Solving the 1D heat equation
Step 3 - Write the discrete equations for all nodes in a matrix
format and solve the system:
The boundary conditions
Solving the 1D heat equation
The discrete approximation of the 1D heat equation:
Numerical stability - for this scheme to be numerically stable,
you have to choose sufficiently small time steps
Numerical accuracy - the numerical accuracy of this scheme is first
order in time and second order in space, i.e. the error scales linearly
with and quadratically with .
Solving the 1D wave equation
Consider the initial-boundary value problem:
Boundary conditions (B. C.’s):
Initial conditions (I. C.’s):
Step 1- Define a discretization in space and time: t
time step
k+1,
time
step k,
time
step k-1, xN = 1.0
x0 = 0
x
Solving the 1D wave equation
Step 2 - Discretize the PDE. Use a central difference scheme
for both time and space derivatives:
Solving for gives:
The Courant numer
Solving the 1D wave equation
Step 3 - Write the discrete equations for all nodes in a matrix
format and solve the system:
A three-level scheme in time
Solving the 1D wave equation
Step 3 - Write the discrete equations for all nodes in a matrix
format and solve the system:
Solving the 1D wave equation
Step 3 - Write the discrete equations for all nodes in a matrix
format and solve the system:
The boundary conditions
Solving the 1D wave equation
A note on time advancing at t =0:
Since the numerical scheme involves three levels of time steps, to advance
to , you need to know the nodal values at and .
Use the two initial conditions to write a new numerical scheme
at :
I.C. 1:
I.C. 2:
Or:
Discrete wave equation
at :
Combine them:
Solving the 1D wave equation
The discrete approximation of the 1D wave equation:
Numerical stability - for this scheme to be numerically stable,
you have to choose sufficiently small time steps
Numerical accuracy - the numerical accuracy of this scheme is
second order in time and second order in space, i.e. the error scales
quadratically with and .
Solving Laplace’s equation
Consider the boundary value problem:
Boundary conditions (B. C.’s):
Step 1- Define a discretization in x and y:
y y
x x
0 1
The numerical mesh
The physical domain N+1 points in x direction,
M+1 point in y direction
Solving Laplace’s equation
Step 2 - Discretize the PDE. Use a central difference scheme
for space derivatives in x and y directions:
If :
The node (n,m) is linked to its 4 neighbouring nodes as illustrated in the finite difference
stencil:
• This finite difference stencil is valid for the
interior of the domain:
• The boundary values are found from the
boundary conditions.
Solving Laplace’s equation
Step 3 - Solve the system by Jacobi iteration:
Take successive neighbour averages at each iteration k+1 th:
Until there is small change in the solution (i.e. the solution has converged), as
measured by:
Matlab codes for numerical solutions of the
heat, the wave and Laplace’s equations:
• You can program the methods explained before in Matlab
(of course, there are many other options, e.g. Python, C+
+, Fortran, etc.)
• All the Matlab codes are uploaded on the course
webpage.
• For each code, you only need to change the input data
and maybe the plotting part. The solver is already there!
• Figures will normally be saved in the same directory as
where you saved the code.
The Matlab code for the 1D heat equation
PDE:
Set the diffusion coefficient here
Set the domain length here
B.C.’s:
Tell the code if the B.C.’s prescribe the value of u
(Dirichlet type ) or its derivative (Neumann type)
Set the values of the B.C.’s on each side I.C.:
Specify an initial value as a function of x
Specify the number of grid points, i.e. N+1,
remember increasing N increases the accuracy.
Set the initial and final time of the
computation. (T_final should be large
enough to get a steady state solution)
Set the number of time steps to be taken to
go from t_0 to t_final. Remember if this
number is too low, you might get a
numerical instability.
The Matlab code for the 1D heat equation
Results: PDE:
0.8 B.C.’s:
0.6 Time increasing
I.C.:
u(x,t)
0.4
0.2 The steady-state solution
0
0 0.2 0.4 0.6 0.8 1
x
The Matlab code for the 1D wave equation
PDE:
Set the wave speed here
Set the domain length here
B.C.’s:
Tell the code if the B.C.’s prescribe the value of u
(Dirichlet type ) or its derivative (Neumann type)
I.C.’s:
Set the values of the B.C.’s on each side
Specify the initial value of u and the initial
time derivative of u as a function of x
Specify the number of grid points, i.e. N+1.
Set the initial and final time of the
computation.
Set the number of time steps to be
taken to go from t_0 to t_final.
The Matlab code for the 1D wave equation
Results:
1
PDE:
0.5
u(x,t)
B.C.’s:
-0.5
-1
0 0.2 0.4 0.6 0.8 1 I.C.’s:
x
Standing waves
The Matlab code for Laplace’s equation
PDE:
Specify the domain size here
B.C.’s:
Set the types of the 4 boundary
Set the B.C.’s on each side of the
rectangle
Specify the number of grid points in x and
y directions, i.e. N+1 and M+1.
The Matlab code for Laplace’s equation
Results:
PDE:
2
1.5
1.5
1
1
0.5 B.C.’s:
0 0.5
2
3
1 2
1 0
y 0 0 x
2 1.8
1.6
1.5 1.4
1.2
1 1
y
0.8
0.5 0.6
0.4
0 0.2
0 1 2 3
x
