Notes Wave Combined

Numerical methods for Earth Sciences
Nicolas Brantut
Department of Earth Sciences
University College London, UK
2
Disclaimer These lecture notes have been put together from a number of sources, in-
cluding Press et al (2007) and Chapra (2012). They are therefore not meant to be original
and are not a substitute to real textbooks; rather, they are meant to be a dense summary
of some of the points that will be covered in class.
Do not distribute!
Contents
1 Integrals 5
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Method of rectangles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Midpoint Method
e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 Trapezoidal rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2 Ordinary Di↵erential Equations 11

2.1 Some generalities about ODEs . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.1 First order ODEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.2 Higher order ODEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.3 Initial value vs. Boundary value problems . . . . . . . . . . . . . . . 11
2.2 Euler’s method for initial value problems . . . . . . . . . . . . . . . . . . . . 12
2.2.1 Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.2 Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3 Modified Euler – Midpoint method . . . . . . . . . . . . . . . . . . . . . . . 13
2.3.1 Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3.2 Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.4 Stability of a numerical method . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.5 Runge-Kutta methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.5.1 Second order RK methods . . . . . . . . . . . . . . . . . . . . . . . . 14
2.5.2 Fourth order RK methods . . . . . . . . . . . . . . . . . . . . . . . . 15
2.6 Systems of ODEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.7 Adaptative stepping in ODE solvers . . . . . . . . . . . . . . . . . . . . . . 16
2.8 Two-points boundary value problems . . . . . . . . . . . . . . . . . . . . . . 19
2.8.1 The finite di↵erence method: some generalities . . . . . . . . . . . . 19
2.8.2 Finite di↵erence method for the heat flow problem (1) . . . . . . . . 20
2.8.3 Finite di↵erence method for the heat flow problem (2) . . . . . . . . 21
010
2.8.4 Going further . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3 Partial Di↵erential Equations 25

3.1 Explicit finite di↵erences for parabolic PDEs . . . . . . . . . . . . . . . . . 25
3.1.1 Concept . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1.2 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1.3 Heterogeneities and nonlinearities . . . . . . . . . . . . . . . . . . . . 26
3.1.4 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.2 Explicit finite di↵erence method for the wave equation . . . . . . . . . . . . 30
3.2.1 Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.2.2 Discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.2.3 Initialisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2.4 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2.5 Numerical dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3
4 CONTENTS
3.3 A note about normalisation . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.4 Benchmarking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.5 Fully implicit finite di↵erence method . . . . . . . . . . . . . . . . . . . . . 35
3.5.1 Concept . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.5.2 Implementing boundary conditions . . . . . . . . . . . . . . . . . . . 36
3.5.3 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.5.4 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.6 The Crank-Nicholson Method . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.6.1 Concept . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.6.2 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.6.3 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.7 Nonuniform gird steps in finite di↵erences methods . . . . . . . . . . . . . . 48
3.8 The road to 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.8.1 Explicit method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.8.2 Implicit methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
e
4 Root finding 53
4.1 Bisection method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.2 Newton-Raphson method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.3 Newton-Raphson in multiple dimensions . . . . . . . . . . . . . . . . . . . . 54
List of Codes 57
Bibliography 59
Chapter 1
1.1
O
Integrals
Introduction
国 1和分
_
敀分
两
Here we want to compute the value of integrals like
0
Z b
I= f (x)dx,
a
or Z x
F (x) = f (x0 )dx0 .
a
In real life problems, analytical, or exact, solutions cannot be found, so we need to
compute numerical approximations. Sometimes, we do not even know a mathematical form
for f (x), which is given to us as a sequence of numerical values (e.g., an accelerogram).
What to do then?
Here are very simple methods, all based on a simple idea that an integral is the “area
under the curve”.
All these methods require the discretisation of the interval over which the integral is
computed: we go from a mathematical continuum to a set of discrete values.
è
1.2 Method of rectangles
Let us discretize the interval [a, b] into N + 1 distinct points, labelled xn (n = 1, . . . , N ),
where x1 = a and xN +1 = b. The spacing between two points is then
b a
x= .
N
The rectangle method consists in approximating the integral I as (Figure 1.1):
I⇡ xf (a) + xf (a + x) + . . . + xf (b x).
A more generic way to formulate this approximation is:

N
X
I⇡ (xn ) ⇥
fEE x.
n=1
Using this rule, we can also estimate the integral F (x) at the points x = xn . By
definition, we have F (a) = 0. Then, we have i
F (a + x) ⇡ f (a) x,
5
公瓯
6 tw CHAPTER 1. INTEGRALS
豳
鬬
1
−1
tatty_t.it
xcosx
−2 Ax 1
Ithreshddty
−3
么 5
−4
At
0 π/5 2π/5 3π/5 4π/5 π
x
Figure 1.1: Integral by method of rectangles: we approximate the integral by the sum of
rectangles of height f (xn ) and width xn+1 xn .
F (a + 2 x) ⇡ F (a + x) + f (a + x) x,
etc. This can be rewritten as:
n
X1
F (xn ) ⇡ f (xk ) x.
k=1
Example: f (x) = x cos(x). The exact answers are known:
I= 2, F (x) = 1 + cos(x) x sin(x).
The piece of Matlab code 1.1 computes approximations of I and F using the method of
rectangles.
Note: The method given above is using points x1 to xN to compute the integral, so we
never use the value of f at point xN +1 = b. We are therefore using the sum of “left”
rectangles. We can rewrite a similar method using the “right” rectangles, i.e., using the
points x2 to xN +1 . The formulas are then:
N
X +1
I⇡ f (xn ) ⇥ x,
n=2
and
n
X
F (xn ) ⇡ f (xk ) x.
k=2
1.3 Midpoint Method 中会1tn

Here we are again using rectangles, but taking their heights not at one or the other top
corners, but in the middle of each interval. The numerical approximation becomes then:
N
X
2ˇ
I⇡ f (xn + x/2) ⇥ x,
n=1
1.4. TRAPEZOIDAL RULE 7
Code 1.1: Integration of x cos(x) using the method of rectangles.
N = 100+1;
x = linspace(0,pi,N);
h = x(2) - x(1);
y = x.*cos(x);
I = 0;
for n=1:N-1
I = I + y(n)*h;
end
disp(['Result is ' num2str(I)])
disp(['Error is ' num2str(abs(I+2))])
% Cumulative integral
F = zeros(size(y));
F(1) = 0;
for n=1:N-1
F(n+1) = F(n) + h*y(n);
end
%exact solution
Fex = -1 + x.*sin(x)+cos(x);
%plots
figure;
subplot 211
plot(x,F,'.', x,Fex,'r-');
xlabel('x');
ylabel('F(x)');
legend('Numerical','Analytical');
subplot 212
plot(x, abs(F-Fex), '.')
xlabel('x');
ylabel('Error | F-F {exact } | ')
and
n
X
F (xn + x/2) ⇡ f (xk + x/2) x.
k=1
Note here that the integral F is then approximated at the middle of each interval.
1.4 Trapezoidal rule

Here we approximate the area under the curve by a series of trapezoids. The approximation
is (Figure 1.2):
N
X f (xn ) + f (xn+1 )
I⇡ ⇥ x,
2
n=1
and
n
X1 f (xk ) + f (xk+1 )
F (xn ) ⇡ x.
2
k=1
For non-uniform spacing between points xn (i.e., x is not constant anymore), basic
8 CHAPTER 1. INTEGRALS
−1
xcosx
−2
−3
−4
0 π/5 2π/5 3π/5 4π/5 π
x
Figure 1.2: Integral by method of trapezoids: we approximate the integral by the sum of
trapezoids with tips at f (xn ) and f (xn+1 ) and width xn+1 xn .
geomtrical considerations can be used to find:

N
X f (xn ) + f (xn+1 )
I⇡ (xn+1 xn ),
2
n=1
and
n
X1 f (xk ) + f (xk+1 )
F (xn ) ⇡ (xk+1 xk ).
2
k=1
Example: Take the same function f (x) = x cos(x) as before. The program 1.2 computes
approximations of I and F using the trapzoidal rule.
1.4. TRAPEZOIDAL RULE 9
Code 1.2: Integration of x cos(x) using the trapezoidal rule.
N = 100+1;
x = linspace(0,pi,N);
h = x(2) - x(1);
y = (x).*cos(x);
I = 0;
for n=1:N-1
I = I + 0.5*(y(n)+y(n+1))*h;
end
disp(['Result is ' num2str(I)])
disp(['Error is ' num2str(abs(I+2))])
% Cumulative integral
F = zeros(size(y));
F(1) = 0;
for n=1:N-1
F(n+1) = F(n) + 0.5*(y(n)+y(n+1))*h;
end
%exact solution
Fex = -1 + x.*sin(x)+cos(x);
%plots
figure;
subplot 211
plot(x,F,'.',x,Fex,'r-');
xlabel('x');
ylabel('F(x)');
legend('Numerical','Analytical');
subplot 212
plot(x, abs(F-Fex), '.')
xlabel('x');
ylabel('Error | F-F {exact } | ')
10 CHAPTER 1. INTEGRALS
0 常钦幽
iiii
10z t.ie
Chapter 2
pnt 偏邀
Ordinary Di↵erential Equations

00
然
型
2.1 Some generalities about ODEs
2.1.1 First order ODEs
Ordinary Di↵erential Equations (ODEs) are of the form:
y 0 = f (t, y), (+BC) (2.1)
where y(t) is the dependent variable, and t is the independent variable (such as time, or
space coordinate). The dependent variable y can be a scalar, or a vector. In the latter
case, we say that Equation 2.1 is a system.
2.1.2 Higher order ODEs

It is always possible to reformulate a higher order ODE into a system of first order ODEs.
For instance, let us consider the following 2nd order ODE:
2 x on
y 00 (t) + p(t)y 0 (t) + q(t)y(t) + r(t) = 0, (2.2)
多成
where p, q, r are functions of t only. If we define
0
z(t) = y (t),
then Equation 2.2 is rewritten as
Oegtizdnieieyy
⇢
z 0 (t) = p(t)z(t) q(t)y(t) r(t)
y 0 (t) = z(t)
Note: This trick is often helpful, but it is sometimes better to solve the second order
problem directly if possible (we’ll see that in a few weeks!).
2.1.3
繼以
边军
Initial value vs. Boundary value problems
业
七
0
Abbreviated IVP and BVP. For IVPs, we are given the value of y at a given t = t0 , and
are looking for y(t) for t > t0 . For BVPs, we are given the value of y(x) at one or several
points in space (often the case in high order ODEs), such as x = x0 and x = x1 , and are
looking for the value of y(x) for x 2 [x0 , x1 ].
11
业
i
12 CHAPTER 2. ORDINARY DIFFERENTIAL EQUATIONS
Code 2.1: Solution of y 0 = x 2xy, y(0) = 1 using Euler’s method. The function f
called in the script can be written in the same script, or, better, placed in a separate file
010
called f.m.
公
% generate vector of x values
1 示5
N = 150;
x = linspace(0, 5, N); i.r
Dx = x(2) - x(1);
% initialise vector of solutions

y = zeros(1, N); 业
E.is
% initial condition
y(1) = 1;
t
yuniaeg.eu
% compute solution using Euler's method
for n=1:N-1 1
y(n+1) = y(n) + Dx*f(x(n), y(n));
end
% plot solution
plot(x, y, 'o')
以上 gatingˇ
hold on
plot(x, 0.5*(1 + exp(-x.ˆ2)))
5ㄩ
function out = f(t, z)
out = t - 2*t*z;
end
numerioalsolutin
ne tt
n
2.2 Euler’s method for initial value problems
2.2.1 Principle
Approximate derivative with slope:
上
y(t + h) y(t)
y 0 (t) ⇡
h
for small h.
Initial condition y(t = t0 ) = y0 . Algorithm is:
yi+1 = yi + hf (ti , yi )
where yi denotes the value of y at point ti .
Example: Let’s look at the following ODE:

y0 = x 2xy, y(x = 0) = 1.
The code 2.1 solves this equation using Euler’s method. Figure 2.1 illustrates the result
using a very coarse discretization (and therefore a very approximate result!).
等
2.2.2 Error
噩
Write the Taylor expansion of y around point t:
h2 00
__0ˋ
y(t + h) = y(t) + hy 0 (t) + y (t) + . . .
2
Local error of truncation is O(h2 ). Global error is O(h) (because we make the sum of all
the local errors, and they cumulate !).
2.3. MODIFIED EULER – MIDPOINT METHOD 13
1
approx. euler method
0.9
ě
0.8
y(x)
closed form
0.7
0.6
0.5
0 0.5 1 1.5 2
x
Figure 2.1: Solution of y 0 = x 2xy, y(0) = 1 using Euler’s method (dots) and closed-
form solution. We visualise that the approximate solution is computed by extrapolating
the local slope of the curve.
2.3
2.3.1
Modified Euler – Midpoint method
Principle
2 OphldngcI_tnt.gl
Similar idea as Euler’s method, but improve our estimate of slope:
thx l
凹
y(t + h) ⇡ y(t) + h ⇥ (slope)
where we take the slope as the derivative at the midpoint between t and t + h:
slope = y 0 (t + h/2) = f (t + h/2, y(t + h/2)).
But we don’t know y(t + h/2) ! The idea is to approximate it as
y(t + h/2) ⇡ y(t) + (h/2)f (t, y(t)).
The algorithm is thus:
yi+1 = yi + hf ti + h/2, yi + (h/2)f (ti , yi ) .
Example: Let’s look at the following ODE, same as before:
y0 = x 2xy, y(x = 0) = 1.
l
Code 2.2 solves this ODE using the midpoint method.
2.3.2 Error
Local error of truncation is O(h3 ). Global error is O(h2 ).
2.4 Stability of a numerical method

Stability 6= Accuracy. A method is unstable if the errors grow exponentially for a problem
for which there is a bounded solution. Stability depends on
• the type of di↵erential equation we are solving,

Code 2.2: Solution of y 0 = x o 2xy, y(0) = 1 using the midpoint method.
% generate vector of x values

N = 150;
x = linspace(0, 5, N);
Dx = x(2) - x(1);
% initialise vector of solutions

y = zeros(1, N);
% initial condition
y(1) = 1;
% compute solution using the midpoint method

for n=1:N-1
y(n+1) = y(n) + Dx*f(x(n) + 0.5*Dx, y(n) + 0.5*Dx*f(x(n), y(n)));
end
% plot solution
plot(x, y, 'o')
hold on
plot(x, 0.5*(1 + exp(-x.ˆ2)))
• the numerical method,
• the step size
Stability can be assessed using the test case
y0 = ay, y(0) = y0 .
Euler’s method is shown to be stable only is the step size is h < 2/a.
2.5 Runge-Kutta methods

杉掠以
Idea is: approximate the local derivative (slope) using a linear combination of slopes taken
at various positions along the nominal step.
General form of RK methods:
n
X
0
yi+1 = yi + h aj kj ,
j=1
上
where aj are constants, and kj are estimates of the slope. Constants are chosen to achieve
a given order of the method, by equating the expression of yi+1 toa Taylor expansion of
y(t + h).
事
2.5.1 Second order RK methods
This is when n = 2.
yi+1 = yi + h(a1 k1 + a2 k2 ).
Several choices are possible:

2.6. SYSTEMS OF ODES 15
贰贰
The midpoint method: It is in fact a type of RK2 method, where a1 = 0, a2 = 1, and
k1 = f (ti , yi ),
✓ ◆
h h
k2 = f ti + , yi + f (ti , yi ) .
2 2
Heun’s method: This is the case when a1 = a2 = 1/2, and

k1 = f (ti , yi ),
k2 = f (ti + h, yi + hf (ti , yi )) .
ǎiiiiiiiiiie
2.5.2 Fourth order RK methods
It is the the standard ODE solver. It is a fourth order method (n = 4). The method is:
e
yi+1 = yi + h(a1 k1 + a2 k2 + a3 k3 + a4 k4 ),
-
where
a1 = 1/6, a2 = 1/3, a3 = 1/3, a4 = 1/6,
and
k1 = f (ti , yi ),
业
k2 = f (ti + h/2, yi + (h/2)k1 ),
是
k3 = f (ti + h/2, yi + (h/2)k2 ),
k4 = f (ti + h, yi + hk3 ).
2.6 Systems of ODEs

Recall that ODEs are in general in the form
y 0 = f (t, y) ( + BC).
丝
In this form, y can be a vector of unknown functions. In practice, it means that we can
write:
y10 = f1 (t, {y1 , . . . , yn }),
来
y20 = f2 (t, {y1 , . . . , yn }),
..
.
上 A
yn0 = fn (t, {y1 , . . . , yn }).
Example The Lorentz system (Lorentz, 1963) consists in three unknown functions
x(t), y(t), z(t) governed by the following ODEs:
dx
= (y x),
dt
堂 io
ee dy
= x(⇢ z) y,
dt
dz
= xy z, (2.3)
dt
where , ⇢and are constant parameters. This system is nonlinear : it contains products
of the unknowns, e.g., the product xy in the third equation.
0
The Matlab code 2.3 is the function that computes the derivatives (output) as func-
tions of the values of (x, y, z) and other parameters.
Taking = 10, ⇢ = 28 and = 8/3, and starting with x(0) = 1, y(0) = 1 and z(0) = 1,
we can solve the equations using RK4. The Matlab script 2.4 does the job.
āo
16
YE CHAPTER 2. ORDINARY DIFFERENTIAL EQUATIONS
Code 2.3: Matlab function that returns derivatives of the vector system of ODE of the
Lorenz problem. 1
function [dX] = f(X, rho,sigma, beta)

0
yiy.in
% [dX] = F(X, rho,sigma, beta)
%
%This function computes the derivative of variables (x,y,z) of the Lorentz
%system.
%
%input:
% X: vector of 3 elements containing the dependent variables, X(1)=x,
% X(2)=y, X(3)=z
0g.ro y.no
列
% rho: parameter of the Lorentz system
% sigma: parameter
% beta: parameter
%
%output:
% dX: column vector of 3 elements containing the time derivatives of the
% dependent variable: dX=(dx,dy,dz).
240
x=X(1);
上
y=X(2);
z=X(3);
dx = sigma*(y-x);
dy = x*(rho-z)-y; tunatindydt tt.gl
dz = x*y-beta*z;
dX=[dx;
dy;
dz];
agaeryugkr
1
以 29
2.7
fg
Adaptative stepping in ODE solvers
巡
When solving equations numerically, how should we choose the step size h ? We need
stability and accuracy. One option is to impose a constant h through the computation,
but this is not necessarily desirable or convenient. Indeed, if the solution y(t) changes
rapidly, we might need small steps, whereas if it remains more or less constant the time
steps could be quite large.
One way to adapt the time step to the behaviour of the solution is to try doubling it
(i.e., use 2h instead of h), and see if we can retain an target accuracy. If not, then we
e 噬
would have to go back and divide the step by 2.
Overall, the idea is to compute the numerical solution twice:
5
• a solution y1 using a step of 2h (one “big” leap from tn to tn + 2h), ix
• a solution y2 using a step of h (two small leaps, from tn to tn + h and then from
tn + 2h).
Note these two solutions are approximations made at a time tn + 2h. But obviously y2
should be more accurate than y1 because the time step is smaller. Now let us give ourselves
a target accuracy ✏, and compare ✏ to the di↵erence between our two solutions,
If
= |y2 y1 |.
四善
< ✏, it means that we can safely double the time step, that is, y1 is almost as
good an approximation as y2 , within our tolerance ✏.
2.7. ADAPTATIVE STEPPING IN ODE SOLVERS 17
Code 2.4: Script that computes the solution of the Lorenz problem using RK4 method.
%parameters:
r=28;
s=10;
b=8/3;
%discretisation:
e
h=0.01;
t=0:h:100;
%initialisation of solution
N=length(t);
X=zeros(3,N); %NOTE: this is now a MATRIX !
%initial condition
X(:,1)=[1;1;1]; %we need a column here ! x,y,z = 1,1,1
%solution using RK4

for i=1:N-1
k1 = f(X(:,i), r,s,b);
k2 = f(X(:,i)+0.5*h*k1, r,s,b);
三
k3 = f(X(:,i)+0.5*h*k2, r,s,b);
k4 = f(X(:,i)+h*k3, r,s,b);
X(:,i+1) = X(:,i) + h*(k1/6 + k2/3 + k3/3 + k4/6);

end
%put solution matrix X into separate variables for simpler handling:

x = X(1,:);
y = X(2,:);
z = X(3,:);
%now we can plot stuff

plot3(x,y,z, 'color', [0 114 178]/255);
axis equal;
grid on;
box on;
xlabel('x'), ylabel('y'), zlabel('z')
view(24,18)
If > ✏, then clearly doubling the time step is clearly not good, but we also don’t
even know if our solution y2 is good enough, so we have to divide the time step by 2, and
try again.
The code 2.5 is an implementation of this technique for the ODE y 0 = x 2xy with
initial condition y(0) = 1. The one thing to notice in the code is how a while loop is used
in conjunction with a counter n instead of a simple for loop, because we do not know in
advance how many time steps we will have (they keep changing)!
Code 2.5: Implementation of adptative time-stepping in the RK4 method.
% initial step size

h = 0.01;
%target error
epsilon = 1e-3;
%initialise solution with aguess for number of points
N = 100;
y = zeros(1,N); x = zeros(1,N);
%max x value
xmax = 5;
% initial conditions
y(1) = 1; x(1) = 0;
% main computation loop

n=1;
while x(n)<=xmax
% increment x:
x(n+1) = x(n) + 2*h;
%solution using 2h
k1 = func(x(n), y(n));
k2 = func(x(n) + h, y(n) + h*k1);
k3 = func(x(n) + h, y(n) + h*k2);
k4 = func(x(n) + 2*h, y(n) + 2*h*k3);
y1 = y(n) + 2*h*(k1/6 + k2/3 + k3/3 + k4/6);
%solution using h (twice)

k2 = func(x(n) + h/2, y(n) + h*k1/2);
k3 = func(x(n) + h/2, y(n) + h*k2/2);
k4 = func(x(n) + h, y(n) + h*k3);
y interm = y(n) + h*(k1/6 + k2/3 + k3/3 + k4/6);
k1 = func(x(n) + h, y interm);
k2 = func(x(n) + h + h/2, y interm + h*k1/2);
k3 = func(x(n) + h + h/2, y interm + h*k2/2);
k4 = func(x(n) + h + h, y interm + h*k3);
y2 = y interm + h*(k1/6 + k2/3 + k3/3 + k4/6);
%take the difference

Delta = abs(y2 - y1);
dittere
%take optimal step size:
if Delta<epsilon
y(n+1) = y2;
finite
h = 2*h;
n = n+1;
else
h = h/2;
end
end
e
%cleanup
x(n:end) = []; y(n:end) = [];
%% plot solution
plot(x, y, '.', x, 0.5*(1+exp(-x.ˆ2)));
xlabel('x'); ylabel('y(x)');
%% function definition
function derivative = func(x,y)
derivative = x - 2*x*y;
照
nti.li
end
2.8. TWO-POINTS BOUNDARY VALUE PROBLEMS 19
上 in
2.8 Two-points boundary value problems
General form:
y 00 (x) + p(x)y 0 (x) + q(x)y(x) + r(x) = 0,
with boundary conditions expressed as either:
___70
1. Dirichlet conditions (values given in terms of y):
EOO.it114
y(x = x0 ) = y0 ,
y(x = x1 ) = y1 ,
2. Neumann conditions (values given in terms of y 0 ):
y 0 (x = x0 ) = q0 ,
y 0 (x = x1 ) = q1 ,
3. Or mixed conditions:
字
i
应
y(x = x0 ) = y0 ,
0
y (x = x1 ) = q1 .
Geophysical example: The steady-state heat flux in the Earth’s crust is governed by
the following second order ODE for temperature T as a function of depth x:
d2 T
142
T
h
= , (2.4)
dx2 k
where h is the volumetric heat production rate (here we will use h = 3 ⇥ 10 6 W m 3 ),
and k is the thermal conductivity of rocks (here we will use k = 4 W m 1 K 1 ).

We will approach this problem using two di↵erent boundary conditions:
1. Constant temperatures at two points:
T (x = 0) = T0 , T (x = L) = TL ,
三
where in practice we will choose T0 = 0 C, and TL = 800 C.
2. Or a constant temperature at one end and a constant flux at the other end:
dT
器上
k = q, T (x = L) = TL ,
1 dx x=0
where in practice we will choose q/k = 36e 3Km 1 and TL = 800 C.
2.8.1 The finite di↵erence method: some generalities

The key idea is to, again, use a so-called “finite di↵erence” approximation of the deriva-
tives, in order to obtain a system of algeabraic equation for the solution at given points,
or nodes. We discretise the space variable into positions xi , and define the value of the
function at those positions as T (xi ) = Ti , for i = 0, . . . , n.
When approximating a derivative with as finite di↵erence, several choices are possible:
巒
Y
20
y
以
CHAPTER 2. ORDINARY DIFFERENTIAL EQUATIONS
放 0
1. Forward finite di↵erence (first order accuracy):
dT T (x + x) T (x) Ti+1 Ti
⇡ = ,
dx x x
-
00
2. Backward finite di↵erence (first order accuracy):
dT T (x) T (x x) Ti Ti 1
⇡ = ,
dx x x
-
3. Centered finite di↵erence (second order accuracy):
扩
0 nA
dT T (x + x) T (x x) Ti+1 Ti 1
⇡ = .
dx 2 x 2 x
For higher order derivatives, the centered finite di↵erence method yields:
d2 T Ti+1 2Ti + Ti 1
2
⇡ ,
dx x2
d3 T (1/2)Ti+2 Ti+1 Ti 1 (1/2)Ti 2
⇡ ,
dx3 x3
d4 T Ti+2 4Ti+1 + 6Ti 4Ti 1 + Ti 2
4
⇡ .
dx x4
上
etc.
2.8.2 Finite di↵erence method for the heat flow problem (1)
The FD scheme for Equation 2.4 is:
0
Ti+1 2Ti + Ti 1 h
= for i = 1, . . . , n 1
x2 k
ǜiǖijiii
and
T0 = T (x = 0) = T0 ,
Tn = T (x = L) = TL ,
where the corresponding spatial nodes are x = xi (i = 0, . . . , n), where x0 = 0 and xn = L.

There are n + 1 unknowns, which are Ti (i = 0, . . . , n), and we have written a lin-
ear system of n + 1 equations relating these unknowns to the parameters and boundary
ÉIÈÉÉ
conditions of the problem. This linear system can be put in a convenient matrix form:
AX = B (2.5)
where 0 1
1 0 0 0 ··· 0 0
B 1 2 1 0 ··· 0 C 0
B C
B 0 1 2 1 ··· 0 C 0
B C
B .. .. .. .. C ..
A=B C,
A⑬
B . . . . C .
B 0 ··· 0 si1 2 1 i 0 C
B ijibi C
@ 0 ··· 0 0 1 2 1 A
0 ··· 0 0 0 0 1
中国
州
i国
2.8. TWO-POINTS BOUNDARY VALUE PROBLEMS 品 21
and the vectors X and B are

0 1 0 1
道
T0 T0
B T1 C B x2 h/k C
B C B C
B T2 C B x2 h/k C
夘
B C B C
X=B .. C, B=B .. C.
B . C B . C
B C B C
@ Tn 1 A @ x2 h/k A
Tn TL
Now, solving the problem amounts to solving the linear system (2.6). One could
imagine doing it like so:
X = A 1 B,
but in practice it is completely unnecessary to compute A 1 explicitly. In Matlab, solving
下
a linear system is done with the forward slash operator:
X = A\B;
For more information about what this does and why, you are strongly encouraged to look
into Matlab’ help file (or type doc \ at the command line prompt).
The Matlab script 2.6 computes the solution to our heat flow problem (2.4), and
compares it with the closed-form solution.
2.8.3 Finite di↵erence method for the heat flow problem (2)
Now we want to solve the problem (2.4), where the temperature T is given at position
x = L, and the temperature gradient dT /dx is given at position x = 0. We can write this
latter boundary condition using our centered finite di↵erence scheme as:
eo
dT T1 T 1
1 dx x=0
⇡
2 x
= q/k.
The problem here is that we have not defined any node at index -1 ! How do we go about
this ? The idea is quite simple: we can just make one extra “ghost” node at index i = 1,
compute everything as normal, and just discard the result for T at that node. In practice,
the FD scheme for Equation 2.4 is:
e
Ti+1 2Ti + Ti 1 h
= for i = 0, . . . , n 1
x2 k
with the boundary conditions:
T1 T 1 = 2 x ⇥ q/k,
Tn = T (x = L) = TL .
We have now a linear system of n + 2 Equations for n + 2 unknowns, which are our
T 1 , T0 , T1 , . . . , Tn . This linear system can be put in a convenient matrix form:
AX = B (2.6)
where 0 1
1 0 1 0 ··· 0 0
B 1 2 1 0 ··· 0 C 0
B C
B 0 1 2 1 ··· 0 C 0
B
B .. .. ..
尹. C
C ..
A=B . . . .. C, .
B C
B 0 ··· 0 1 2 1 0 C
B C
@ 0 ··· 0 0 1 2 1 A
0 ··· 0 0 0 0 1
Code 2.6: Solution of boundary-value problem (2.4) with fixed T boundary conditions
using finite di↵erences.
% parameters
L = 35e3;
N = 5;
Dx = L/(N-1);
x = 0:Dx:L;
h = 3e-6; %W/mˆ3
k = 4; %W/m/K
%boundary conditions
TS = 0;
TL = 800;
0N
%input matrices
A = zeros(N);
A(1,1) = 1;
A(N,N) = 1;
for m=2:N-1
A(m,m-1) = 1;
A(m,m) = -2;
A(m,m+1) = 1;
end
B = zeros(N,1);
B(1) = TS;
B(N) = TL;
for m=2:N-1
B(m) = -h*Dxˆ2/k;
end
%solve
T = A\B;
plot(x, T, 'o')
hold on
Tclosed = TS + (TL-TS)*(x/L) - 0.5*(h/k)*(x-L).*x;
plot(x, Tclosed,'r')
and the vectors X and B are

0 1 0 1
B
B
B
T 1
T0 C
C
C
B
B
B
0
2 xq/k
x2 h/kC
C
C
B T1 C B x2 h/kC
X=B .. C, B=B .. C.
B . C B . C
B C B C
@ Tn 1 A @ x h/k A
2
Tn TL
The Matlab script 2.7 solves this problem, and compares it with the closed-form solu-
tion.
2.8.4 Going further

Equipped with this method, more complicated problems can be solved, for instance a
similar heat flow problem with a nonhomogeneous distribution of heat sources:
d2 T h(x)
= ,
dx2 k
2.8. TWO-POINTS BOUNDARY VALUE PROBLEMS 23
where h(x) can be given by either a mathematical formula, or even by acutal data collected
somewhere !
Similarly, if the right hand side of Equation (2.4) depends on the temperature T itself,
we can also use a finite di↵erence method to arrive at a system of algebraic equations
relating the temperature vector at all nodes T0 , . . . , Tn to the parameters of the problem.
If the right-hand-side dependency is linear, we can get the solution in one step like we did
above. If the dependency is nonlinear, we need to use iterative methods to solve for the
temperature, but that is another story.
Code 2.7: Solution of boundary-value problem (2.4) with fixed T boundary condition at
x = L and fixed dT /dx at x = 0, using finite di↵erences.
% parameters
L = 35e3;
N = 500;
Dx = L/(N-1);
x = 0:Dx:L;
h = 3e-6; %W/mˆ3
k = 4; %W/m/K
TL = 800;
qoverk = 36e-3; %K/m
%input matrices
A = zeros(N+1);
x
A(1,1) = -1;
A(1,3) = 1;
A(N+1,N+1) = 1;
for m=2:N
A(m, m-1) = 1;
A(m, m) = -2;
是选
A(m, m+1) = 1;
end
B = zeros(N+1,1);
B(1) = 2*Dx*qoverk; 1
B(N+1) = TL;
for m=2:N
B(m) = -h*Dxˆ2/k;
end
0127
%solve
T = A\B; 1
竺 X.T
%delete first element (ghost node at index -1)
00
T(1) = [];
%plot
plot(x, T, 'o');
hold on
维
Tclosed = TL + qoverk*(x-L) - 0.5*(h/k)*(x.ˆ2 - Lˆ2);
plot(x, Tclosed, 'r')
是
州 004
为
之
00
⑦ 业
ˋ
Chapter 3 哵
Partial Di↵erential Equations
Tx
3.1 Explicit finite di↵erences for parabolic PDEs
3.1.1 Concept
✓ ◆
tī
A typical example (of geophysical interest !) is the heat equation (here in one dimension):
器上
@T @ @T
= ↵ , (3.1)
@t @x @x
where T is the temperature, t is the time, x is the spatial coordinate, and ↵ is the thermal
di↵usivity. For a constant value of ↵, Equation (3.1) is rewritten as
的 0幽
@T @2T
1
- = ↵ . (3.2)
@t @x2
We discretise the space along nodes at positions xi , and time with steps tn . The grid
spacing is denoted x and the time step is denoted t. We denote T (xi , tn ) ⇡ Tin (this
is not an equality because Tin is the approximate, numerical solution).
A classic method is to use centered finite di↵erences in space, and forward finite dif-
ferences in time (FTCS), which leads to the following numerical scheme:
⑤ 00
Tin+1 Tin n
Ti+1 2Tin + Tin 1
前的
=↵ , (3.3)
t x2
which can be rewritten as:
Tin+1 = Tin + s(Ti+1
n
2Tin + Tin 1 ), (3.4)
É 027 三
where
s = ↵ t/ x2 .
发放
This scheme is called explicit because we can compute the temperature field at time step
n + 1 based only on its value at time step n.
Boundary conditions are implemented in just the same way as explained in Lecture 3
for boundary value problems. Note that they can now depend on time !
3.1.2 Stability
显此
The stability of the explicit finite di↵erence method for the di↵usion equation depends on
the time step, grid spacing and di↵usivity. One can show that stability is ensured when
-
which corresponds to a limitation on the maximum
0
s < 1/2,
time step:
eǎi
1 x2
t< .
2 ↵
25
26 CHAPTER 3. PARTIAL DIFFERENTIAL EQUATIONS
Important: Remember that stability has nothing to do with accuracy, so a stable

scheme might still be very inaccurate if the time or space steps are taken too large.
3.1.3 Heterogeneities and nonlinearities

For nonconstant ↵, the FTCS scheme for Equation (3.1) is written:
0_
Tin+1 = Tin + ↵ (T n Tin ) ↵i 1/2 (Ti
n
Tin 1 ) ,
x2 i+1/2 i+1
where
↵i+1/2 = ↵(xi +
e x/2), ↵i 1/2 = ↵(xi
o
x/2)
Note that if ↵ is not a function of space ↵(x) but of temperature itself ↵(T ), the numerical
scheme is the same but we then use
n 1 n
↵i+1/2 ⇡ ↵(Ti+1 ) + ↵(Tin )
热
2
and
1
↵in 1/2 ⇡ ↵(Tin ) + ↵(Tin 1 ) .
2
Here, the stability condition in terms of maximum time step is:
⇢
1 x2
t < tmax = min ,
i 2 ↵i+1/2
which essentially means that we are limited by the largest value of ↵ encountered in our
problem.
3.1.4 Examples
三
Cooling of oceanic lithosphere
We solve Equation (3.2) with the following initial and boundary conditions:
ee_
T (x, t = 0) = Tm , T (x = 0, t) = Ttop , T (x = L, t) = Tm ,
with numerical values ↵ = 10 6 m2 s 1 , Tm = 1200 C, Ttop = 0 C, and L = 100 km.

We will use I space steps and N time steps. The numerical scheme is:
✓ ◆
00
n+1 n ↵ t n
Ti = Ti + (Ti+1 2Tin + Tin 1 ), i = 1, . . . , I 1.
x2
0
and we have the boundary conditions
e
T0n = Ttop , TIn = Tm , n = 0, . . . , N
and initial condition

Ti0 = Tm , i = 0, . . . , I.
If we want x = 1 km, we’ll have a total of I = L/ x + 1 = 101 nodes. We therefore
have to choose t small enough to ensure stability of the scheme:
1 x2
t< tmax = ⇡ 15800 years.
2 ↵
The Matlab code 3.1 implements this method and generates the plots.
At
me
3.1. EXPLICIT FINITE DIFFERENCES FOR PARABOLIC PDES 27
Cooling with nonuniform di↵usivity

We solve the same problem as above but we now use
↵(x) = ↵0 exp( x/h).
The numerical scheme is the following
o
t
Tin+1 = Tin + ↵ Tn Tin ↵i Tin Tin 1 , i = 1, . . . , I 1,
x2 i+1/2 i+1 1/2
where
↵i±1/2 = ↵0 exp (xi ± x/2)/h),
and we have the boundary conditions
T0n = Ttop , TIn = Tm , n = 0, . . . , N
and initial condition

Ti0 = Tm , i = 0, . . . , I.
If we want x = 1 km, we’ll have a total of I = L/ x + 1 = 101 nodes. We therefore
have to choose t small enough to ensure stability of the scheme:
1 x2
t< tmax = ⇡ 75400 years.
2 ↵0
The Matlab code 3.2 implements this method and generates the plots.
Code 3.1: Matlab implementation of explicit FD method to solve the plate cooling
problem.
%% Parameters
L = 100e3; %m
alpha = 1e-6; %mˆ2/s
Tm = 1200; %C
Ttop = 0; %C
tmax = 50e6* 365*24*3600; % in seconds !!
%% Discretisation
Dx = 1e3;
Dtmax = 0.5*Dxˆ2/alpha;
Dt = 0.5*Dtmax;
%time vector:
t = 0:Dt:tmax;
%space vector:
x = 0:Dx:L;
%number of nodes in time and space:

N = length(t);
I = length(x);
%% initialisation
T = zeros(N,I);
%initial conditions:
T(1,:) = Tm;
T(1,1) = Ttop;
%% iterations
for n=1:N-1
a
%boundary conditions:
T(n+1, 1) = Ttop;
T(n+1, I) = Tm;
%interior points:
for in=2:I-1
T(n+1, in) = T(n, in) + (alpha*Dt/Dxˆ2)*(T(n,in+1) ...
- 2*T(n,in) + T(n, in-1));
end
end
%% plots
figure;
[c, h] = contour(t/(1e6*24*3600*365), -x/1e3, T', 'k');
clabel(c, h);
xlabel('time (My)')
ylabel('depth (km)')
3.1. EXPLICIT FINITE DIFFERENCES FOR PARABOLIC PDES 29
Code 3.2: Matlab implementation of explicit FD method to solve the plate cooling
problem with heterogeneous heat di↵usivity.
%% Parameters
L = 100e3; %m
alpha0 = 2e-6; %mˆ2/s
h=50e3; %m
Tm = 1200; %C
Ttop = 0; %C
tmax = 50e6* 365*24*3600; % in seconds !!
%% Discretisation
Dx = 1e3;
Dtmax = 0.5*Dxˆ2/alpha0;
Dt = 0.5*Dtmax;
%time vector:
t = 0:Dt:tmax;
%space vector:
x = 0:Dx:L;
%number of nodes in time and space:

N = length(t);
I = length(x);
%% initialisation
T = zeros(N,I);
%initial conditions:
e
T(1,:) = Tm;
T(1,1) = Ttop;
%% iterations
for n=1:N-1
%boundary conditions:
T(n+1, 1) = Ttop;
T(n+1, I) = Tm;
%interior points:
for in=2:I-1
0
a ip = alpha0*exp(-(x(in)+0.5*Dx)/h);
a im = alpha0*exp(-(x(in)-0.5*Dx)/h);
T(n+1, in) = T(n, in) + (Dt/Dxˆ2)*(a ip *(T(n,in+1) - T(n,in)) -...
a im *(T(n,in) - T(n,in-1)));
end
end
%% plots
figure;
[c, hc] = contour(t/(1e6*24*3600*365), -x/1e3, T', 'k');
clabel(c, hc);
xlabel('time (My)')
ylabel('depth (km)')
3.2 Explicit finite di↵erence method for the wave equation

Solving the wave equation numerically is one of the main task in modern seismology. Here
避
we can only brush the surface of the problem, and the reader is referred to the excellent
book “Computational seismology” by Igel for a more complete reference and examples on
the subject.
3.2.1 Problem statement

We will consider a simple case of a scalar wave equation in 1D, which governs the evolution
of a displacement u as a function of time t and space x following
eo
@2u 2
2@ u
= c , (3.5)
@t2 @x2
where c is the wavespeed in the material. We need to complement this equation with
initial and boundary conditions. Initial conditions are
0
u(t = 0, x) = f (x), (3.6)
@u
(t = 0, x) = g(x), 以 (3.7)
@t
where f and g are two given functions of position x. Di↵erent types of boundary conditions
can be used, either by fixing the value of u (Dirichlet boundary condition) or the value
of @u/@t (Neumann boundary condition) at the boundaries of the medium. Here, we will
develop our algorithm taking the example of Dirichlet boundary conditions at the two
end-points, x = 0 and x = L:
u(t, x = L) = b(t),
0
u(t, x = 0) = a(t), (3.8)
(3.9)
where a and b are given functions of time t.
3.2.2 Discretisation
Let us choose a uniformly spaced mesh in both space and time, and define time steps
tn = n t and space steps xi = i x, where n = 0, . . . , N , i = 0, . . . , I, and t is the
(constant) time step size, and x is the (constant) space step size.
t2
y
We approximate the derivatives in Equation 3.5 with centered finite di↵erences:
a
@2u un+1 2uni + uni 1
x
y
i
E
⇡ ,
@t2 tn ,xi t2
@2u uni+1 2uni + uni 1
⇡ ,
@x2 tn ,xi x2
烾
where we used the notation uni = u(tn , xi ). Combining these approximations with the
wave equation for indices i = 1, . . . , I 1 and n = 2, . . . , N , we arrive at:
un+1 2 n 2
)uni + 2 n
uni 1
0nn.ee
i = ui+1 + 2(1 ui 1 , (3.10)
where
t
=c . (3.11)
x
Equation (3.10) is an explicit formulation, since it allows us to compute the value of u at
a given time step based on values of u at previous time steps.
3.2. EXPLICIT FINITE DIFFERENCE METHOD FOR THE WAVE EQUATION 31
Equation (3.10) is valid for i = 1, . . . , I 1 and n = 2, . . . , N (i.e., we have excluded

the boundary points as well as the first two time steps). The boundary conditions, in
discretised form, are written simply as:
un0 = an = a(tn ),
unI = bn = b(tn ),
so that the value of u at these nodes is always known. Therefore, at positions x1 and xI 1,
our algorithm (3.10) becomes
un+1
1 = 2 n
u2 + 2(1 2
)un1 + 2
an un1 1
,
un+1
I 1 =
2
bn + 2(1 2
)unI 1 + 2 n
uI 2 unI 1
1.
The unknowns of the problem form a vector

0 1
u1
B C
u = @ ... A (3.12)
uI 1
and Equation (3.10) can be rewritten in matrix form as
000
un+1 = B · un un 1
+ cn , (for n = 2, . . . , N ) (3.13)
where 0 1
2(1 2) 2 0 0
B 2 2(1 2) 2 0 C
B C
B .. C
B=B
B 0 2 2(1 2) . 0 C
C (3.14)
B .. .. C
@ 0 0 . . 2 A
0 0 0 2 2(1 2)
and 0 1
2a
n
B 0 C
B C
Bn .. C
c =B . C. (3.15)
B C
@ 0 A
2b
n
3.2.3 Initialisation
Our algorithm (3.13) uses values of u at time steps n and n 1 to compute the value of
u at time step n + 1: we need to have information on two preceeding steps in order to
march forward in time. How do we initialise the procedure?
We are given two initial conditions, (3.6) and (3.7). The initial condition (3.6) gives
us directly what the displacement is at t = 0, which in terms of our discretised variables
yields
u0i = fi = f (xi ). (3.16)
Now in order to be able to compute u2i (the third time step), we also need to know what
u1i is. But this is not given to us directly! We only know @u/@t at t = 0 (Equation 3.7).
We can use that information to produce a numerical estimate of u1i as follows. Firstly, let
us remark that, by definition:
u1i = u(0 + t, xi ).
Using a Taylor expansion, valid for small t, we get

@u 1 @2u
u(0 + t, xi ) ⇡ u(0, xi ) + t+ t2 + . . . .
@t 0,xi 2 @t2 0,xi
Now we can use our condition (3.7),

@u
(0, xi ) = gi = g(xi ),
@t
and use our discretised wave equation to substitute for the second derivative:
@2u c2
⇡ (u0 2u0i + u0i 1 ),
@t2 0,xi x2 i+1
where we remark again that u00 = a0 and u0I = b0 and u0i = fi . Combining, we obtain:
2 2
u1i = fi+1 + (1 + 2
)fi + fi 1 + gi t, (3.17)
2 2
which can be put in matrix form analogue to (3.13):
1 1
u1 = B · f + tg + c0 . (3.18)
2 2
Example: Solve the wave equation (3.5) for x 2 [0, 1] with boundary conditions
2 2
u(t, 0) = 0, u(t, 1) = 0 and initial conditions u(0, x) = e (x 1/2) /(0.02) and u̇(0, x) = 0.
The code 3.3 shows a full implementation of the explicit finite di↵erence method to find
this solution.
3.2.4 Stability
We can analyse the stability of the numerical method by choosing a trial solution of the
wave equation, and determining the conditions under which the numerical approximation
of that solution will (or will not) grow exponentially with time. Here, we choose a solution
of the form u(t, x) = est e jkx , where s is a growth rate, k is a wavenumber, and j 2 = 1.
Using our discretised time tn = n t and space xi = i x, the trial solution at time step n
and space step i is
uni = esn t e jki x ,
and substituting this form in Equation (3.10) (ignoring boundary conditions for the time
being), we obtain
e2s t
2es t ( 2 ik x
e + 2(1 2
)) + 2
e ik x
1 = 0.
The growth of the solution between two successive time step is given by a = es t , and we
know that the solution will be unstable is |a| > 1. In terms of a, and after some algebraic
manipulations, the above equation is rewritten as
a2 2a 1 2 2
sin2 (k x/2) + 1 = 0,
which can be solved to find
p
2
a± = ± 2 1, =1 2 sin2 (k x/2).
We observe that |a| > 1 if > 1, and |a| = 1 is  1. The condition  1 is equivalent
to c t x < 1. We have therefore arrived at our stability condition:
t
Stability requires c  1.
x
This condition is also known as the Courant-Friedrichs-Lewy (CFL) criterion.
3.2. EXPLICIT FINITE DIFFERENCE METHOD FOR THE WAVE EQUATION 33
3.2.5 Numerical dispersion

The discretisation in space and time, and subsequent approximations of the derivatives,
lead to potential problems and inconsistencies with the expected physical properties of the
solution. One important such problem is the fact that the phase velocity of the numerical
solution depends on the discretisation, so that di↵erent wavelengths do not propagate at
the same velocity. This phenomenon is called numerical dispersion. Let us now analyse
why this is.
Assuming we have chosen a stable numerical scheme, i.e., c t x < 1, then we know
that |a| = 1, so that s must be purely imaginary. This means that our trial solution is a
plane wave, and we can write a as
a = ei! t ,
where ! is the pulsation of the plane wave solution. From the solution a+ given above1 ,
we can express the pulsation as a function of wavenumber and the other parameters:
2
cos(! t) = =1 2 sin2 (k x/2),
so that
1 ⇥ ⇤
!=acos 1 2 2 sin2 (k x/2) .
t
Now, we can ask ourselves what is the phase velocity of our numerical solution, and
compare it with the (true) phase velocity c. By definition, the phase velocity is the ratio
of pulsation to wavenumber; let us call c⇤ the phase velocity of the numerical wave:
! 1 ⇥ ⇤
c⇤ = = acos 1 2 2
sin2 (k x/2) .
k ! t
We directly observe that c⇤ 6= c! Our numerical wave has a phase velocity that depends
on the discretisation, and on the wavenumber k! It means there is dispersion: here the
disperstion is purely numerical, since our initial wave equation (the true physical model)
assumed no such dispersion.
We can finally ask ourselves under which conditions we can approach c⇤ ! c (this
would be desirable to have a more faithful numerical solution). For small x, we can
approximate
sin(k x/2) ⇠ k x/2,
so that
1 ⇥ ⇤
c⇤ ⇠ acos 1 c2 k 2 t2 /2 .
k t
Now assuming also that t approaches 0, we use the Taylor expansion of acos near 1 to
show that
1 p
c⇤ ⇠ 2 ⇥ (c2 k 2 t2 /2),
k t
so that we have shown that
lim c⇤ = c.
t!0
x!0
What this means is that we need to have both small time and space steps to minimise the
phenomenon of numerical dispersion.
Example: When running the code 3.3, we observe that for long computational times
the shape of the propagating pulse becomes distorted, with “trailing” wavelets behind the
main pulse. This distortion is directly due to numerical dispersion, with high wavenumbers
propagating at lower phase velocity (hence the “trail” of high frequency wavelets).
1
Using a would lead to negative !, which is unphysical.
Code 3.3: Solution of the wave equation using finite di↵erences.
% solving the wave equation in 1d using the CTCS method
% parameters
v = 1;
L = 1;
% boundary conditions
a = @(t) 0;
b = @(t) 0;
% initial conditions
initial u = @(x) exp(-(x-0.5).ˆ2/(0.02ˆ2));
initial udot = @(x) 0;
% discretisation
I = 301;
x = linspace(0,L,I+2); %add boundary points
Dx = x(2)-x(1);
xint = x(2:end-1); %without boundary points
% CFL criterion
Dtmax = Dx/v;
Dt = 0.5*Dtmax;
N = 4000;
t = Dt*(0:N-1);
% Matrix B
B = zeros(I);
s2 = (v*Dt/Dx)ˆ2;
B(1,1) = 2*(1-s2);
B(1,2) = s2;
for n=2:I-1
B(n,n-1) = s2;
B(n,n) = 2*(1-s2);
B(n,n+1) = s2;
end
B(I,I-1) = s2;
B(I,I) = 2*(1-s2);
% initialisation
u = zeros(I,N);
u(:,1) = initial u(xint);
c = zeros(I,1);
c(1) = s2*a(t(1))/2;
c(I) = s2*b(t(1))/2;
g = initial udot(xint);
u(:,2) = 0.5*B*u(:,1) + Dt*g + c;
% main loop
for n=2:N-1
c(1) = s2*a(t(n));
c(I) = s2*b(t(n));
u(:,n+1) = B*u(:,n) - u(:,n-1) + c;
end
%plot solution
figure;
for k=1:10:N
plot(xint, u(:,k));
ylim([-1,1]);
drawnow;
end
3.3. A NOTE ABOUT NORMALISATION 35
3.3 A note about normalisation

[This is physics ! But essential for numerical computing]
Equation (3.2) is written with parameters and variables having physical dimensions
(by default, we assume SI units: T is K, t in seconds, etc). This is often not very general,
and physical behaviour typically depends on nondimensional combinations of parameters,
not individual ones. In terms of numerical computing, there can also be problems if our
code combines both very large and very small numbers (due to numerical accuracy as well
as the inherent approximations when dealing with floats).
It is always best to normalise the governing equations, based on natural characteristic
scales of our problem. For instance, let’s use Equation (3.2) with boundary conditions
T (x = 0, t) = T0 , T (x = L, t) = TL ,
a natural choice of characteristic scales is TL T0 for the temperature, L for space, and
L2 /↵ for time (this comes from dimensional analysis!). We can define the following nondi-
mensional quantities:
T T0 x t
T̃ = , x̃ = , t̃ = 2 ,
TL T0 L L /↵
and then rewrite Equation (3.2) in a completely equivalent (but more general) way:
@ T̃ @ 2 T̃
= ,
@ t̃ @ x̃2
with boundary conditions
T̃ (x̃ = 0, t̃) = 0, T̃ (x̃ = 1, t̃) = 1.
This is clearly easier to handle ! We have much less parameters to deal with and to fix
in our codes. Plus, most variables are now of the order of 1, which helps in terms of
numerical computing !
3.4 Benchmarking
This is a key concept in numerical computing. Because we use numerical methods to
obtain solutions when closed-form solutions cannot be found, we always have to ask our-
selves whether we can trust our solutions or not. Even if the numerical algorithm can
be mathematically proven to be a good approximation, our own implementation may be
bugged. So what to do ?
The idea is to benchmark our code by using it in situations where a closed-form solution
does exist, and check its validity and accuracy.
In the case of our heat equation, most (if not all) closed-form solutions for a range of
geometries are given by Carslaw and Jaeger (1959). If your code does not converge to the
given closed-form solution in a simple case, then there is no chance you can trust it to do
anything else !
Note that benchmarking is never a proof that a code works in a general case. But
it is an essential step when developing any numerical code, otherwise people will be very
suspicious !
3.5 Fully implicit finite di↵erence method

-
3.5.1 Concept
Here again we use the typical example (of geophysical interest!) of the heat equation in
one dimension, equation (3.2).
装
2 oe
One way to solve this problem numerically is to use centered finite di↵erences in space,
and backward finite di↵erences in time (BTCS), which leads to the following numerical
scheme:
Tin+1 Tin T n+1 2Tin+1 + Tin+1
= ↵ i+1 1
, (3.19)
t x2
which can be rewritten as:
n+1
sTi+1 + (1 + 2s)Tin+1 + sTin+1 n
1 = Ti , (3.20)
a0
where
s = ↵ t/ x2 .
This scheme is called implicit because we compute the temperature field at time step n + 1
and node i based on its value at other nodes i + 1 and i 1 at the same time step n + 1.
The BTCS scheme is first order in time, and second order in space.
Equation (3.20) is essentially a linear system, which can be put in matrix form as
AX = B, (3.21)
where A is a tridiagonal matrix, X is the vector formed by {T0n+1 , . . . , TIn+1 }, and B is

typically {T0n , . . . , TIn }, with the first and last elements potentially modified due to the
boundary conditions. In general B also includes the source terms.
3.5.2 Implementing boundary conditions

Dirichlet (imposing T ): For instance, this could be:
T (x = 0, t) = Tb (t),
X
which in terms of indices would translate into
T0n = Tb (tn ).
i ii
é
The matrix A is then of the form:
0 1
1 0 0 0 0 ···
B s (1 + 2s) s 0 0 ··· C
B C
A=B 0 s (1 + 2s) s 0 ··· C ,
@ A
.. .. .. .. .. ..
. . . . . .
and 0 1
T0n+1
B T1n+1 C
B C
X=B T2n+1 C,
@ A
..
.
and 0 1
Tb (tn+1 )
B T1n C
B C
B=B T2n C.
@ A
..
.
Note here that I do not show the last lines of A, X and B since they also change
depending on the boundary conditions !
3.5. FULLY IMPLICIT FINITE DIFFERENCE METHOD
Neumann (imposing @T /@x): For instance, this could be:
@T
0
5经 37
= q0 (t),
@x x=0
which in terms of indices would translate into
T1n+1 T n+1
1
= q0 (tn+1 ).
2 x
The matrix A is then of the form:
0 1
1 0 1 0 0 ···
B s (1 + 2s) s 0 0 ··· C
B C
A=B 0 s (1 + 2s) s 0 ··· C ,
@ A
.. .. .. .. .. ..
. . . . . .
and 0 1
T n+11
B T0n+1 C
B C
昌
B T1n+1 C
X=B C,
B T2n+1 C
@ A
..
.
and 0 1
2 xq0 (tn+1 )
B T0n C
B C
B T1n C
B=B C.
B T2n C
@ A
..
.
3.5.3 Stability
Again here we assume a solution of the form
Tin = e n t jki x
e ,
and substitute it into Equation (3.20). Defining the amplification factor as
Tin+1 t
a= =e ,
Tin
after some algebra we find that
1
a= ,
1 + 4s sin2 (k x/2)
e___
which ensures that |a| < 1 whatever the time and space steps. The method is then always
stable. But not necessarily accurate ! Again, accuracy is not the same as stability.
3.5.4 Example
We solve the heat equation
o
@T @2T
=↵ 2,
@t @x
with boundary conditions:
@T q(t)
T (x = L) = Tb , = ,
@x x=0 2↵⇢c
with q(t) = q0 exp( t/t0 ). The parameter values are ↵ = 10 6 m2 s 1 , ⇢c = 2.6 ⇥

106 Pa K 1 , q0 = 30 ⇥ 106 Pa m s 1 , t0 = 0.05 s, and L = 0.01 m.
Let us write the full numerical scheme, using the fully implicit finite di↵erence method.
We can split the description into three parts:
• For i = 0, . . . , I 1 (i.e., for interior points, not directly a↵ected by the boundary
conditions):
n+1
sTi+1 + (1 + 2s)Tin+1 sTin+1 n
1 = Ti .
• Boundary condition at x = 0: We introduce the extra “ghost” node at i = 1, and

obtain
__
q(tn+1 )
T1n+1 T n+1
1 = x .
↵⇢c
• Boundary condition at x = L (i.e., last element i = I):
TIn+1 = Tb .
In terms of matrices, the problems is now
AX = B,
where 0 1
1 0 1 0 0 ··· 0
B s (1 + 2s) s 0 0 ··· 0 C
B C
B 0 s (1 + 2s) s 0 ··· 0 C
B C
B .. .. .. .. .. .. C
A=B . . . . . . C,
B C
B 0 ··· 0 s (1 + 2s) s 0 C
B C
@ 0 ··· 0 0 s (1 + 2s) s A
0 ··· 0 0 0 0 1
and 0 1
T n+11
B
B T0n+1 C
C
B T1n+1 C
B C
B T2n+1 C
X=B C,
B .. C
B . C
B C
@ T n+1 A
I 1
TIn+1
and 0 1
1
xq(tn+1 )/(↵⇢c)
B T0n C
B C
B T1n C
B C
B T2n C
B=B C.
B .. C
B . C
B C
@ TIn+1 A
1
Tb
The code 3.4 is a Matlab implementation of that method.
3.5. FULLY IMPLICIT FINITE DIFFERENCE METHOD 39
Benchmarking: At this point, how can we make sure that our code gives a reasonable
solution to the problem? There is a priori no analytical solution available2 . But a solution
can be found in a slightly simpler case, when q(t) = q0 (constant flux at x = 0). It reads
(Carslaw and Jaeger, 1959):
"r ✓ ◆#
q0 ↵t x2 /4↵t x x
T (x, t) = e erfc p .
↵⇢c ⇡ 2 2 ↵t
The code 3.5 computes the numerical and closed-form solutions, and compares them. This
looks fine, and it gives us confidence that our code can be used in a more general case.
But be careful: it is by no means a proof that the code works! In general, more than one
benchmark solutions are used, to encompass all the possible cases and check all potential
problems.
2
In general there isn’t: otherwise we would not be solving it numerically!
Code 3.4: Matlab implementation of implicit FD method to solve the problem of example
3.5.4.

rhoc = 2.6e6; %Pa/K
q0 = 30e6; %Pa m /s
t0 = 0.05; %s
Tb = 0; %to simplify
%domain size:
L= 0.01; %m
%computed time:
tmax = 1; %s
I = 200+1;
x = linspace(0,L,I);
Dx = (x(2)-x(1));
N = 200;
t = linspace(0,tmax,N);
Dt = (t(2)- t(1));
%initialise T:
T = zeros(I+1,N);
s = alpha*Dt/Dxˆ2; %this does not change !
%get the matrix A: also constant over time !

A = zeros(I+1);
%fill A:
A(1,1) = -1;
A(1,3) = 1;
for i=2:I
A(i,i-1) = -s;
A(i,i) = (1+2*s);
A(i,i+1) = -s;
end
A(I+1,I+1) = 1;
%compute for all times:

for n=1:N-1
%assign vector B:
B = T(:,n);
%the first element is different:
B(1) = -Dx*q0/(alpha*rhoc)*exp(-t(n+1)/t0);
%and the last as well:
B(I+1) = Tb;
%solve the system AXn = -Xn+1 + Bn+1
X = A\B;
%store the result into our array T:
T(:,n+1) = X;
end
%extract only T0 to TI (not T-1):

T(1,:) = []; %this deletes the first line
%plots
figure;
pcolor(x,t,T');
xlabel('x (m)');
ylabel('t (s)');
colorbar;
3.5. FULLY IMPLICIT FINITE DIFFERENCE METHOD 41
Code 3.5: Matlab implementation of implicit FD method and corresponding analytical

benchmark for the problem of example 3.5.4.

rhoc = 2.6e6; %Pa/K
q0 = 3e6; %Pa m /s
t0 = 0.05; %s
%domain size:
L= 0.01; %m
%computed time:
tmax = 1; %s
I = 200+1;
Dx = (x(2)-x(1));
N = 200;
Dt = (t(2)- t(1));
%initialise T:
T = zeros(I+1,N);

%get the matrix A: also constant over time !
A = zeros(I+1);
%fill A:
A(1,1) = -1;
A(1,3) = 1;
for i=2:I
A(i,i-1) = -s;
A(i,i) = (1+2*s);
A(i,i+1) = -s;
end
A(I+1,I+1) = 1;
%compute for all times:

for n=1:N-1
%assign vector B:
B = T(:,n);
%the first element is different:
B(1) = -Dx*q0/(alpha*rhoc);
%and the last as well:
B(I+1) = Tb;
X = A\B;
T(:,n+1) = X;
end

%closed form:
Tclosed = zeros(I,N);
for i=1:I
Tclosed(i,:) = q0/alpha/rhoc *(sqrt(alpha*t/pi).*exp(-x(i)ˆ2./(4*alpha*t)) ...
-x(i)/2 *erfc(x(i)./(2*sqrt(alpha*t))) );
end
figure;
plot(t,T(1,:),'o', t, Tclosed(1,:),'-')
3.6 The Crank-Nicholson Method

3.6.1 Concept
Here again, we solve the heat equation (3.2)
逗层
The idea here is to use the average of implicit and explicit finite di↵erence schemes.
This amounts to using a centered time, centered space scheme, which should be second
order accurate in both space and time. The Crank-Nicholson scheme reads:
的 4
1
Tin+1 Tin
n_i
↵
- = T n+1 2Tin+1 + Tin+1 n
1 + Ti+1 2Tin + Tin 1 ,
t 2 x2 i+1
- -
which is in fact centered in time at tn+1/2 . This can be more conveniently rewritten as:
at
s n+1 s n+1 s n s
II
Ti+1 + (1 + s)Tin+1 s)Tin + Tin 1 ,
n
Ti 1 = Ti+1 + (1 (3.22)
2 2 2 2
which also correspond to a linear system of the form
PX = B,
where
B = QY + C.
With these notations, we have

0 1 0 1
T0n+1 T0n
B C B C
X = @ ... A , Y = @ ... A ,
TIn+1 TIn
and P and Q being tridiagonal matrices. The additional matrix C corresponds to addi-
tional terms associated with boundary conditions and source terms. The best to see what
they look like to go through our example again.
3.6.2 Stability
If we assume a solution of the form
Tin = e n t jki x
e ,
and substitute it into Equation 3.22, we can work out the amplification factor defined as
a = Tin+1 /Tin = e t
.
After some algebra, we find that
1 2s sin2 (k x/2)
a= ,
1 + 2s sin2 (k x/2)
which ensures that |a| < 1 regardless of the time and space steps. The method is therefore
unconditionnally stable. But again not necessarily accurate!
3.6. THE CRANK-NICHOLSON METHOD 43
3.6.3 Example
We solve the heat equation
@T @2T
=↵ 2,
@t @x
with boundary conditions:
@T q(t)
T (x = L) = Tb , = ,
@x x=0 2↵⇢c
with q(t) = q0 exp( t/t0 ). The parameter values are ↵ = 10 6 m2 s 1 , ⇢c = 2.6 ⇥
106 Pa K 1 , q0 = 30 ⇥ 106 Pa m s 1 , t0 = 0.05 s, and L = 0.01 m.
Let us write the full numerical scheme, using the Crank-Nicholson method. We can
split the description into three parts:
• For i = 0, . . . , I 1 (i.e., for interior points, not directly a↵ected by the boundary
conditions):
s n+1 s n+1 s n s
T + (1 + s)Tin+1 T = Ti+1 + (1 s)Tin + Tin 1 .
2 i+1 2 i 1 2 2
• Boundary condition at x = 0: We introduce the extra “ghost” node at i = 1.
Remember that the finite di↵erence approximation for the spatial derivative should
be averaged between time steps n and n + 1:
!
1 T1n+1 T n+1 1 T1n T n1 q(tn+1 ) + q(tn )
+ = ,
2 2 x 2 x 2↵⇢c
Tip
which is more conveniently rewritten as:
x q(tn+1 ) + q(tn )
T1n+1 T n+1 n
1 =T 1 T1n .
0 ↵⇢c
• Boundary condition at x = L (i.e., last element i = I):

TIn+1 = Tb .
PED
ei 1
In terms of matrices, the problems is now 化
PX = QY + C,
where 0 1
1 0 1 0 0 ··· 0
B s/2 (1 + s) s/2 0 0 ··· 0 C
B C
B 0 s/2 (1 + s) s/2 0 ··· 0 C
B C
B .. .. .. .. .. .. C
P=B . . . . . . C,
B C
B 0 ··· 0 s/2 (1 + s) s/2 0 C
B C
@ 0 ··· 0 0 s/2 (1 + s) s/2 A
i 01
0 ··· 0 0 0 0 1
Qi
and 0
1 0 1 0 0 ··· 0
B s/2 (1 s) s/2 0 0 ··· 0 C
B C
B 0 s/2 (1 s) s/2 0 ··· 0 C
B C
B C
Q=B
B
..
.
..
.
..
.
..
.
..
.
..
. C,
C n
B 0 ··· 0 s/2 (1 s) s/2 0 C
B C
@ 0 ··· 0 0 s/2 (1 s) s/2 A
0 ··· 0 0 0 0 0
品品
44
l 上品
CHAPTER 3. PARTIAL DIFFERENTIAL EQUATIONS
and 0 1 0 1
T n+1 1 T n1
B
B T0n+1 C
C
B
B T0n C
C
B T1n+1 C B T1n C
B C B C
B
X=B T02n+1 C
C, 为 B
Y =B T2n C
C,
B .. C B .. C
B . C B . C
B C B C
@ T n+1 A @ Tn A
I 1 I 1
TIn+1 TIn
and 0 1
x(q(tn+1 ) + q(tn ))/(↵⇢c)
B 0 C
B C
B
C=B EIE ... si C
C.
B C
@ 0 A
Tb
n
The codes 3.6 and 3.7 give a Matlab implementation of that method.
Benchmarking: The benchmarking can be done in the same way as for the other meth-
ods, by using a known solution to a simpler problem. This is left for homework.
A short script to make some good quality figures The code 3.8 can be run after
you have computed the temperature as a function of time and space in the example. This
is an example of how to plot a complex solution and extract only the key aspects to show
the readers.
Code 3.6: Matlab implementation of Crank-Nicholson method for the heat flow problem
of example 3.6.3. This part of the code only shows how the fill in the matrices.

rhoc = 2.6e6; %Pa/K
q0 = 30e6; %Pa m /s
t0 = 0.05; %s
%domain size:
L= 0.01; %m
%computed time:
tmax = 1; %s
%generate vector of positions:

I = 255+1;
Dx = (x(2)-x(1));
%generate vector of times:

N = 256;
Dt = (t(2)- t(1));
%CFL number:
%initialise T:
T = zeros(I+1,N);
%initial condition at t=0:
T(I+1,1) = Tb;
% fill the matrix P: also constant over time !

P = zeros(I+1);
P(1,1) = -1;
P(1,3) = 1;
for i=2:I
P(i,i-1) = -0.5*s;
P(i,i) = (1+s);
P(i,i+1) = -0.5*s;
end
P(I+1,I+1) = 1;
% fill the matrix Q: also constant over time !

Q = zeros(I+1);
Q(1,1) = 1;
Q(1,3) = -1;
for i=2:I
Q(i,i-1) = 0.5*s;
Q(i,i) = (1-s);
Q(i,i+1) = 0.5*s;
end
Code 3.7: Matlab implementation of Crank-Nicholson method for the heat flow problem
of example 3.6.3. This part of the code shows the computation of the solution and plots
it.
% COMPUTATION
for n=1:N-1
%intialise rhs vector B:
Y = T(:,n);
B = Q*Y;
%update first element of B:
B(1) = B(1) - 2*Dx*q0/(alpha*rhoc)*0.5*(exp(-t(n)/t0)+exp(-t(n+1)/t0));
%and the last one:
B(I+1) = Tb;

X = P\B;

T(:,n+1) = X;
end

%plots
figure;
pcolor(x,t,T');
xlabel('x (m)');
ylabel('t (s)');
colorbar;
Code 3.8: Matlab script to produce a good quality figure from computations performed
in codes 3.6 and 3.7.
figure;
subplot 211
plot(t, T(1,:));
hold all;
plot(t, T(10,:));
plot(t, T(50,:));
plot(t, T(100,:));
xlabel('{\itt} (s)');
ylabel('{\itT} (K)');
legend('{\itx} = 0',...
['{\itx} = ' num2str(x(10)*1e3) ' mm'],...
['{\itx} = ' num2str(x(50)*1e3) ' mm'],...
['{\itx} = ' num2str(x(100)*1e3) ' mm']);
subplot 212
plot(x*1e3, T(:,4));
hold all;
plot(x*1e3, T(:,16));
plot(x*1e3, T(:,32));
plot(x*1e3, T(:,64));
xlim([0 5])
xlabel('{\itx} (mm)');
ylabel('{\itT} (K)');
legend(['{\itt} = ' num2str(t(4)) ' s'],...
['{\itt} = ' num2str(t(16)) ' s'],...
['{\itt} = ' num2str(t(32)) ' s'],...
['{\itt} = ' num2str(t(64)) ' s']);
3.7 Nonuniform gird steps in finite di↵erences methods

Moving forward to solving PDEs using finite di↵erence methods, we are now going to
explore more general numerical scheme accounting for nonuniform grid spacing.
Let us use the example of the heat equation (here in one dimension), equation (3.2).
spacing is denoted xi and the time step is denoted t. We denote T (xi , tn ) ⇡ Tin (this
Here, we did not make the choice of having a constant value of x, and therefore we
need to track the value of the spacing between each node xi . Hence, we denote
xi = xi+1 xi .
(Note that if we have I nodes, with i ranging from 0 to I 1, we have I 1 spacings)

Using nonconstant step sizes, the finite di↵erence approximation of the second deriva-
tive in x is given by:
✓ ◆
@2T 2 Ti+1 Ti Ti Ti 1
⇡ ,
@x2 xi + xi 1 xi xi 1
so that the explicit finite di↵erence method is written as

✓ n ◆
Tin+1 Tin 2↵ Ti+1 Tin Tin Tin 1
= ,
t xi + xi 1 xi xi 1
and the implicit finite di↵erence method is then

!
n+1
Tin+1 Tin 2↵ Ti+1 Tin+1 Tin+1 Tin+1
1
= .
t xi + xi 1 xi xi 1
Example The code 3.9 is an implementation of the explicit finite di↵erence method
solving the heat equation with the following boundary conditions:
T (t, 0) = T0 + A sin(2⇡t/⌧ ), T (t, L) = T0 ,
i.e., a surface temperature that oscillates with a period ⌧ and amplitude A around a fixed
value T0 , and a temperature at depth L that is fixed at T0 .
Furthermore, if we also have nonuniform di↵usivity, then the approximation becomes:
x
✓ ◆ ✓ ◆
e
@ @T 2 Ti+1 Ti Ti Ti 1
↵(x) ⇡ ↵i+1/2 ↵i 1/2 ,
@x @x xi + xi 1 xi xi 1
where we have as before
↵i+1/2 = ↵(xi + xi /2),

↵i 1/2 = ↵(xi xi 1 /2).
Then, the explicit and implicit finite di↵erence scheme can be found as above, by taking
the value of T at time step n (explicit) or n + 1 (implicit).
Note that in the implicit method with nonuniform grid, the construction of the matrix
problem is more subtle since the elements of the matrix depend on the position (index i).
3.7. NONUNIFORM GIRD STEPS IN FINITE DIFFERENCES METHODS 49
Code 3.9: Implementation of nonconstant grid spacing to solve the heat equation with
oscillating imposed T at the top and constant T at some position x = L.
%% Parameters
L = 6; %in m
alpha = 1e-6; %in mˆ2/s
day = 24 * 3600; %in s
tfinal = 30 * day;
tau = 1*day;
A = 15;
T0 = 10;
%% Discretise
Dx = logspace(-2, log10(L/4),25);
x = cumsum(Dx);
Dtmax = 0.5*min(Dx)ˆ2/alpha;
Dt = 0.5*Dtmax;
t = 0:Dt:tfinal;
%number of nodes:
J = length(x);
N = length(t);
%% initialise
T = zeros(N, J);
% initial condition
T(1, :) = T0;
T(1, 1) = Tsurf(t(1), T0, A, tau);
%% solve
for n=1:N-1
T(n+1, 1) = Tsurf(t(n+1), T0, A, tau);
T(n+1, J) = T0;
%interior points:
for i=2:J-1
T(n+1, i) = T(n,i) + (alpha*Dt*2/(Dx(i) + Dx(i-1)))*(...
(T(n, i+1) - T(n,i))/Dx(i) - (T(n, i) - T(n,i-1))/Dx(i-1));
end
end
%% plots
figure;
pcolor(t,-x,T')
xlabel('time (s)')
ylabel('depth (m)')
colorbar;
在孔器器
50 器匙 CHAPTER 3. PARTIAL DIFFERENTIAL EQUATIONS
3.8 The road to 2D

The heat equation in 2D is written:
✓ ◆
@T @2T @2T
=↵ + .
@t @x2 @y 2
The idea is to discretise the space into I points in the x direction, and J points in the y
direction:
x = xi , i = 0, . . . , I 1
y = yj , j = 0, . . . , J 1.
There is then a total 趟

of I ⇥ J nodes.
以
3.8.1 Explicit method
n
If we use an explicit finite di↵erence method, this is quite easy. Just define an array Ti,j
(i.e., a 3D array), and the scheme is written as:
n+1 ✓ ◆
ㄨ 000
Ti,j n
Ti,j n n + Tn n n + Tn
Ti+1,j 2Ti,j i 1,j Ti,j+1 2Ti,j i,j 1
=↵ + .
t x2 y2
This can be solved almost directly, for instance in Matlab we will have to define:
_b e
% initialise array of T:
T = zeros(N,I,J);
and then iterate like
tnn.int
E
T(n+1, i, j) = (a*Dt)*(T(n,i+1,j)-2*T(n,i,j)+T(n,i-1,j))/Dxˆ2 ...
+ (a*Dt)*(T(n,i,j+1)-2*T(n,i,j)+T(n,i,j+1))/Dyˆ2;
e
Again, we have to be careful about the stability, and we need:
t< tmax =
1 min{ x,
2 ↵
y}2
.
3.8.2 Implicit methods

This is not conceptually more difficult, but requires lots of care when implementing. The
main di↵erence and complication of implicit methods is that we want to solve a linear
system of equations, so we need to have the unknowns as a single vector X. If you observe
what we did above for the explicit methods, you will realise that the temperature T at
time step n + 1 is no longer a vector, but a 2D array ! So we cannot use that approach
directly.
Idea: Define a unique index k for each node (i, j), i.e., define a map to convert any pair
(i, j) into a single index k. One such map is:
k = i + I ⇥ j, (3.23)
and is illustrated in Figure 3.1.

3.8. THE ROAD TO 2D 51
i\j 0 1 2
0
0 3 6
1
1 4 7
2
2 5 8
Figure 3.1: Illustration of the map (3.23). Indices i are shwon in orange, j in blue, and
the corresponding k are in green.
In terms of index pairs (i, j), the fully implicit finite di↵erence method reads:
n+1 n+1 n+1
!
Ti,j n
Ti,j Ti+1,j 2Ti,j + Tin+1
1,j
n+1
Ti,j+1 n+1
2Ti,j n+1
+ Ti,j 1
=↵ + .
t x2 y2
Our map (3.23) implies that:
(i, j) ! k,
(i + 1, j) ! k + 1,
(i 1, j) ! k 1,
(i, j + 1) ! k + I,
(i, j 1) ! k I,
so that the scheme, in terms of indices k, is

↵ t n+1 ↵ t n+1
(T 2Tkn+1 + Tkn+1
1) (T 2Tkn+1 + Tkn+1 n+1
I ) + Tk = Tkn .
x2 k+1 y 2 k+I
This can be more conveniently rewritten as
n+1 n+1
sy Tk+I sx Tk+1 + (1 + 2sx + 2sy )Tkn+1 sx Tkn+1
1 sy Tkn+1 n
I = Tk ,
Which is a nice linear system of the form AX = B ! Of course, here A is more complicated
than in the 1D case, but it is in fact only a few extra terms that are non zero.
In practice, this approach requires a lot of attention to details, especially when num-
bering your nodes and handling boundary conditions.
As usual, the advice here is the write everything very carefully on paper before doing
any computing.
1
Chapter 4 平 iii.i
0_to
Root finding
iitinniin.ae
Sometimes we want to find the zero(s) of a given function f (x), i.e., solve
0s
f (x) = 0.
If f is not simple it is in general not possible to find an expression for x is terms of common
mathematical functions. There are many numerical methods to achieve approximations
of x.
了
The first step is always to plot the function and visualise what we should expect.
Elementary analysis shoudl also be performed to check is we expect multiple roots, if we 三
t
need complex numbers, etc.
4.1 Bisection method

0Y 分纪
This method always works (but cna be slow). It requires an initial bracketing of the
solution. Say we know two numbers a, b such that f (a)  0  f (b). Then
㕺
2
三3nifiii
• compute c = (a + b)/2, iii.si
• if f (c) < 0, assign a c, else b c,
• repeat until f (c) is equal to zero within some tolerance.
Code 4.1 is a matlab implementation of the method to find the zero of f (x) = e x2 x.
4.2 Newton-Raphson method 㚇

The idea here is to start searching for a root with an initial guess x0 , and successively
determine improved approximations. In principle, we want ot find an increment such
that
9
f (x0 + ) = 0.
If we are not too far from the root, should be small and we use the approximation:
f (x0 + ) = f (x0 ) + f 0 (x) + . . .
Equating this to zero, we find
___ˊ0_ =
f (x0 )
f 0 (x0 )
.
ㄨ1
ˋ_ˊ
不 1
业
We then assign our new approximation x1 = x0 + and keep iterating until a given
tolerance is reached.
2
Code 4.2 is a matlab implementation of the method to find the zero of f (x) = e x x.
兴
53
灶少
__110
700 ㄨ2
ˋ_众们
州
54 CHAPTER 4. ROOT FINDING
上
Code 4.1: Root of f (x) = e x2 x using bisection method.
器
主
% function definition
f = @(x) exp(-xˆ2)-x;
%initial bracket
co
a=1;
b=0;
c=0.5*(a+b);
%bisection using 1e-8 tolerance 时 ftii

while abs(f(c))>1e-8
if f(c)<0
a=c;
三
else
b=c;
end
end
c=0.5*(a+b);
5
disp(['Root is x=' num2str(c)])
f 51
Code 4.2: Root of f (x) = e x2 x using Newton-Raphson method.
% function and its derivative

f = @(x) exp(-xˆ2)-x;
fp = @(x) -2*x*exp(-xˆ2) -1;
%initial guess
x0=0;
while abs(f(x0))>1e-8
x0 = x0 - f(x0)/fp(x0);
end
disp(['Root is x=' num2str(c)])
4.3 Newton-Raphson in multiple dimensions
If we are looking for the zero(s) of a vector function of multiple variable, e.g., find (x, y)
such that
莁燕
f (x, y) = 0,
g(x, y) = 0,
then the situation is a lot more complicated. Newton-Raphson is the only simple method
to get there, if our initial guess is not too bad.
装步
How does that work? Say we have an approximation for our root (x, y), and we want
to improve on this by finding the increment ( x , y ) such that
f (x + x, y + y) = 0,
g(x + x , y + y ) = 0.
4.3. NEWTON-RAPHSON IN MULTIPLE DIMENSIONS
问 55
Writing the Taylor expansion gives

@f @f
f (x + x, y + y) = f (x, y) + x + y + ...,
@x @y
g(x + x, y + y) = g(x, y) +
@g
x+
@g
y + .... 2ㄨ 1
@x @y
The increments can be found by solving the system
蹒
@f @f
f (x, y) = x + y,
@x @y
@g @g
g(x, y) = x+ y,
@x @y
which can be rewritten in matrix form as
J = F,
where 0 1
@f @f
i _i
B @x @y C ,
J = @ @g
Jocab matrix
@g A
@x @y
and = ( x , y ) and F = (f (x, y), g(x, y)). The solution of the linear system is obtained
using Matlab’s solver with the \ operator. Then we can increment the solution and repeat
the process.
Example We want to find the (local) minimum of function
p(x, y) = x4 + 2y 2 + xy 6y + 5x.
We expect that the total derivative is zero at the minimum, so we look to find the coor-
dinate (x, y) such that
@p
@x = f (x, y) = 4x3 + y + 5 = 0
@p
@y = g(x, y) = x + 4y 6 = 0
Code 4.3 is a matlab implementation of the Newton-Raphson method to find the zero of
this system.
B
A
KI 4
吧近
定
n1
56 CHAPTER 4. ROOT FINDING
eii7
Code 4.3: Finding the local minmum of x4 + 2y 2 + xy 6y + 5x using Newton-Raphson

method.
p = @(x,y) xˆ4 + 2*yˆ2 + x*y - 6*y + 5*x;

F = @(x,y) [4*xˆ3 + y + 5;
4*y + x - 6];
J = @(x,y) [12*xˆ2, 1;
1, 4];
xx = linspace(-3,3);
yy = linspace(-3,3);
ff = zeros(100);
for i=1:length(xx)
for j=1:length(yy)
pp(i,j) = p(xx(i),yy(j));
end
end
figure;
surf(xx,yy,pp);
hold on
x0=-2;
y0=2;
X=[x0;y0];
F0 = F(x0,y0);
while sum(abs(F0))>1e-8
D = J(x0,y0)\(-F0);
X = X + D;
F0 = F(X(1),X(2));
end
plot3(X(2), X(1), p(X(1),X(2)),'ko');

List of Codes
1.1 Integration of x cos(x) using the method of rectangles. . . . . . . . . . . . . 7

1.2 Integration of x cos(x) using the trapezoidal rule. . . . . . . . . . . . . . . . 9
2.1 Solution of y 0 = x 2xy, y(0) = 1 using Euler’s method. The function f

called in the script can be written in the same script, or, better, placed in
a separate file called f.m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2 Solution of y 0 = x 2xy, y(0) = 1 using the midpoint method. . . . . . . 14
2.3 Matlab function that returns derivatives of the vector system of ODE of
the Lorenz problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4 Script that computes the solution of the Lorenz problem using RK4 method. 17
2.5 Implementation of adptative time-stepping in the RK4 method. . . . . . . . 18
2.6 Solution of boundary-value problem (2.4) with fixed T boundary conditions
using finite di↵erences. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.7 Solution of boundary-value problem (2.4) with fixed T boundary condition
at x = L and fixed dT /dx at x = 0, using finite di↵erences. . . . . . . . . . 24
3.1 Matlab implementation of explicit FD method to solve the plate cooling

problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2 Matlab implementation of explicit FD method to solve the plate cooling
problem with heterogeneous heat di↵usivity. . . . . . . . . . . . . . . . . . . 29
3.3 Solution of the wave equation using finite di↵erences. . . . . . . . . . . . . . 34
3.4 Matlab implementation of implicit FD method to solve the problem of ex-
ample 3.5.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.5 Matlab implementation of implicit FD method and corresponding analytical
benchmark for the problem of example 3.5.4. . . . . . . . . . . . . . . . . . 41
3.6 Matlab implementation of Crank-Nicholson method for the heat flow prob-
lem of example 3.6.3. This part of the code only shows how the fill in the
matrices. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.7 Matlab implementation of Crank-Nicholson method for the heat flow prob-
lem of example 3.6.3. This part of the code shows the computation of the
solution and plots it. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.8 Matlab script to produce a good quality figure from computations per-
formed in codes 3.6 and 3.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.9 Implementation of nonconstant grid spacing to solve the heat equation with
oscillating imposed T at the top and constant T at some position x = L. . . 49
2
4.1 Root of f (x) = e x x using bisection method. . . . . . . . . . . . . . . . 54
2
4.2 Root of f (x) = e x x using Newton-Raphson method. . . . . . . . . . . 54
4.3 Finding the local minmum of x4 +2y 2 +xy 6y +5x using Newton-Raphson
method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
57
58 LIST OF CODES
Bibliography
Carslaw, H. S. and J. C. Jaeger (1959), Conduction of heat in solids, 2nd edition,

Oxford University Press, New York.
Chapra, S. C. (2012), Applied Numerical Methods with Matlab for Engineers and Scien-
tists, 3rd edition, McGraw Hill, New York, USA.
Igel, H. (2017), Computational Seismology, Oxford University Press, Oxford, UK.
Lorentz, E. N. (1963), Deterministic nonperiodic flow, J. Atm. Sci., 20, 130–141.
Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery (2007),

Numerical Recipes, 3rd edition, Cambridge University Press, Cambirdge, UK.
59
Numerical methods for Earth Sciences
Appendix: Julia codes
Nicolas Brantut
Department of Earth Sciences
University College London, UK
2
List of Codes
1 Integration of x cos(x) using the method of rectangles. . . . . . . . . . . . . 4

2 Integration of x cos(x) using the trapezoidal rule. . . . . . . . . . . . . . . . 5
3 Solution of y 0 = x 2xy, y(0) = 1 using Euler’s method. . . . . . . . . . . 6
4 Solution of y 0 = x 2xy, y(0) = 1 using the midpoint method. . . . . . . 6
5 Script that computes the solution of the Lorenz problem using RK4 method. 7
6 Implementation of adptative time-stepping in the RK4 method. . . . . . . . 8
7 Solution of boundary-value problem with fixed T boundary conditions using
finite differences. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
8 Solution of boundary-value problem with fixed T boundary condition at
x = L and fixed dT /dx at x = 0, using finite differences. . . . . . . . . . . . 10
9 Matlab implementation of explicit FD method to solve the plate cooling
problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
10 Matlab implementation of explicit FD method to solve the plate cooling
problem with heterogeneous heat diffusivity. . . . . . . . . . . . . . . . . . . 12
11 Solution of the wave equation using finite differences. . . . . . . . . . . . . . 13
12 Matlab implementation of implicit FD method. . . . . . . . . . . . . . . . . 14
13 Matlab implementation of Crank-Nicholson method for the heat flow problem. 15
14 Matlab script to produce a good quality figure from computations performed
in code 13. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
15 Implementation of nonconstant grid spacing to solve the heat equation with
oscillating imposed T at the top and constant T at some position x = L. . . 17
2
16 Root of f (x) = e x x using bisection method. . . . . . . . . . . . . . . . 18
2
17 Root of f (x) = e x x using Newton-Raphson method. . . . . . . . . . . . 18
18 Finding the local minmum of x4 + 2y 2 + xy 6y + 5x using Newton-Raphson
method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3
4 LIST OF CODES
Code 1: Integration of x cos(x) using the method of rectangles.
⌥ ⌅
function integ1(func,a,b,N)
x = range(a, stop=b, length=N)
I = 0.0
for n in 1:N-1
I += func(x[n])*(x[n+1] - x[n])
end
return I
end
function integ1c(func,a,b,N)
F = zeros(size(x))
for n in 1:N-1
F[n+1] = F[n] + func(x[n])*(x[n+1] - x[n])
end
return F
end
# our test function

f(x) = x*cos(x)
fi(x) = -1 + x*sin(x) + cos(x) #integral
# compute integral from 0 to ⇡ with N points

N = 100+1
I = integ1(f,0,⇡,N)
println("Result is $(I)")
println("Error is $(abs(I+2))")
# compute cumulative integral

F = integ1c(f,0,⇡,N)
# plots
using PyPlot
figure()
subplot(2,1,1)
plot(x, F, "k.", x, fi.(x), "r-")
xlabel(L"x")
ylabel(L"F(x)")
legend(["Numerical", "Analytical"])
subplot(2,1,2)
plot(x, abs.(F - fi.(x)), "k.")
xlabel(L"x")
⌃ ⇧
ylabel("Error \$F-F_{exact}\$")
LIST OF CODES 5
Code 2: Integration of x cos(x) using the trapezoidal rule.
⌥ ⌅
function integ2(func,a,b,N)
I = 0.0
for n in 1:N-1
I += 0.5*(func(x[n]) + func(x[n+1]))*(x[n+1] - x[n])
end
return I
end
function integ2c(func,a,b,N)
F = zeros(size(x))
for n in 1:N-1
F[n+1] = F[n] + 0.5*(func(x[n]) + func(x[n+1]))*(x[n+1] - x[n])
end
return F
end
# our test function

f(x) = x*cos(x)
fi(x) = -1 + x*sin(x) + cos(x) #integral
# compute integral from 0 to ⇡ with N points

N = 100+1
I = integ2(f,0,⇡,N)
println("Result is $(I)")
println("Error is $(abs(I+2))")
# compute cumulative integral

F = integ2c(f,0,⇡,N)
# plots
using PyPlot
figure()
subplot(2,1,1)
plot(x, F, "k.", x, fi.(x), "r-")
xlabel(L"x")
ylabel(L"F(x)")
legend(["Numerical", "Analytical"])
subplot(2,1,2)
plot(x, abs.(F - fi.(x)), "k.")
xlabel(L"x")
⌃ ⇧
ylabel("Error \$F-F_{exact}\$")
6 LIST OF CODES
Code 3: Solution of y 0 = x 2xy, y(0) = 1 using Euler’s method.
⌥ ⌅
function euler1(f, x, f0)
y = zeros(size(x))
y[1] = y0
N = length(x)
for n=1:N-1
y[n+1] = y[n] + (x[n+1]-x[n])*f(x[n],y[n])
end
return y
end
function func(t, z)
return t - 2*t*z
end
# generate vector of x values

N = 150
x = range(0,stop=5,length=N)
y0 = 1
# compute solution using Euler's method

y = euler1(func, x, y0)
# plot solution
e
using PyPlot
figure()
plot(x, y, "ko")
⌃ ⇧
plot(x, 0.5*(1 .+ exp.(-x.ˆ2)), "r")
Code 4: Solution of y 0 = x 2xy, y(0) = 1 using the midpoint method.
⌥ ⌅
function midpoint(f, x, f0)
y = zeros(size(x))
y[1] = y0
N = length(x)
for n=1:N-1
x = x[n+1]-x[n]
y[n+1] = y[n] + x*f(x[n] + 0.5* x,
y[n] + 0.5* x*f(x[n],y[n]))
end
return y
end
function func(t, z)
return t - 2*t*z
end
# generate vector of x values

N = 150
x = range(0,stop=5,length=N)
y0 = 1
# compute solution using Euler's method

y = midpoint(func, x, y0)
# plot solution
using PyPlot
figure()
plot(x, y, "ko")
⌃ ⇧
plot(x, 0.5*(1 .+ exp.(-x.ˆ2)), "r")
LIST OF CODES 7
Code 5: Script that computes the solution of the Lorenz problem using RK4 method.
⌥ ⌅
function florenz(X, ⇢, , )
x=X[1]
y=X[2]
z=X[3]
dx = *(y-x)
dy = x*(⇢-z)-y
dz = x*y- *z
return [dx;
dy;
dz]
end
function rk4(f, x, y0)

N = size(x,1)
h = x[2]-x[1]
y = Array{typeof(y0), 1}(undef,N)
y[1] = y0
for n=1:N-1
k1 = f(x[n], y[n])
k2 = f(x[n]+0.5*h, y[n]+0.5*h*k1)
k3 = f(x[n]+0.5*h, y[n]+0.5*h*k2)
k4 = f(x[n]+h, y[n]+h*k3)
y[n+1] = y[n] + h*(k1/6+k2/3+k3/3+k4/6)

end
return y
end
#parameters
r=28
s=10
b=8/3
N=10000 + 1
t=range(0,stop=100,length=N)
Y0=[1.;1.;1.]
#solution
Y = rk4((x,y)->florenz(y,r,s,b), t, Y0)
#extract x,y,z from solution

x = [Z[1] for Z in Y]
y = [Z[2] for Z in Y]
z = [Z[3] for Z in Y]
#plot
using PyPlot
figure()
plot3D(x,y,z)
⌃ ⇧
axis("scaled")
8 LIST OF CODES
Code 6: Implementation of adptative time-stepping in the RK4 method.
⌥ ⌅
function rk4step(f, x, y, h)
k1 = f(x, y)
k2 = f(x+h/2, y+h*k1/2)
k3 = f(x+h/2, y+h*k2/2)
k4 = f(x+h, y+h*k3)
return h*(k1/6+k2/3+k3/3+k4/6)
end
function rk4adapt(f, a, b, y0, h0, ✏)

# initial step size
h = h0
# initialise solution
y = Array{typeof(y0), 1}(undef,1)
y[1] = y0
# initialise x array
x = Array{typeof(h0), 1}(undef,1)
x[1] = a
# main loop
n = 1
while x[n]<=b
# solution using step 2h

y1 = y[n] + rk4step(f, x[n], y[n], 2h)
# solution using step h

y_interm = y[n] + rk4step(f, x[n], y[n], h)
y2 = y_interm + rk4step(f,
x[n],
y_interm,
h)
# difference
= abs(y2 - y1)
# adapt step size

if <✏
push!(y, y2)
push!(x, x[n]+2h)
h *= 2
n += 1
else
h /= 2
end
end
return x,y
end
f(x,y) = x-2x*y
x,y = rk4adapt(f, 0, 5, 1.0, 0.01, 0.01)
#plot solution
using PyPlot
figure()
plot(x, y, "k.", x, 0.5*(1 .+exp.(-x.ˆ2)), "r")
xlabel(L"x")
⌃ ⇧
ylabel(L"y(x)")
LIST OF CODES 9
Code 7: Solution of boundary-value problem with fixed T boundary conditions using

finite differences.
⌥ ⌅
# parameters
L = 35e3
N = 50
x = range(0, stop=L, length=N)
Dx = x[2]-x[1]
h = 3e-6
k = 4.0
# BCs
TS = 0.0
TL = 800.0
# matrices
A = zeros(N,N)
A[1,1] = 1
A[N,N] = 1
for m=2:N-1
A[m,m-1] = 1
A[m,m] = -2
A[m,m+1] = 1
end
B = zeros(N)
B[1] = TS
B[N] = TL
for m=2:N-1
B[m] = -h*Dxˆ2/k
end
# solve
T = A\B
# plot
using PyPlot
figure()
plot(x, T, "k.")
Tclosed = TS .+ (TL-TS)*(x/L) .- 0.5*(h/k)*(x.-L).*x
⌃ ⇧
plot(x, Tclosed, "r")
10 LIST OF CODES
Code 8: Solution of boundary-value problem with fixed T boundary condition at x = L

and fixed dT /dx at x = 0, using finite differences.
⌥ ⌅
# parameters
L = 35e3
N = 50
x = range(0, stop=L, length=N)
Dx = x[2]-x[1]
h = 3e-6
k = 4.0
# BCs
TL = 800.0
qoverk = 36e-3
# matrices
A = zeros(N+1,N+1)
A[1,1] = -1
A[1,3] = 1
A[N+1,N+1] = 1
for m=2:N
A[m,m-1] = 1
A[m,m] = -2
A[m,m+1] = 1
end
B = zeros(N+1)
B[1] = 2*Dx*qoverk
B[N+1] = TL
for m=2:N
B[m] = -h*Dxˆ2/k
end
# solve
T = A\B
# delete first element

deleteat!(T,1)
# plot
using PyPlot
figure()
plot(x/1e3, T, "k.")
Tclosed = TL .+ qoverk*(x.-L) .- 0.5*(h/k)*(x.ˆ2 .- Lˆ2)
⌃ ⇧
plot(x/1e3, Tclosed, "r")
LIST OF CODES 11
Code 9: Matlab implementation of explicit FD method to solve the plate cooling problem.
⌥ ⌅
# parameters
L = 100e3
↵ = 1e-6
Tm = 1200
Tt = 0
tmax = 50e6*365*24*3600
# discretisation
x = 1e3
tmax = 0.5* xˆ2/↵
t = tmax/2
t = 0: t:tmax
x = 0: x:L
N = size(t,1)
I = size(x,1)
# solution
T = zeros(N,I)
T[1,:] .= Tm
T[1,1] = Tt
for n=1:N-1
T[n+1, 1] = Tt
T[n+1, I] = Tm
for i=2:I-1
T[n+1,i] = T[n,i] + (↵* t/ xˆ2)*(T[n,i+1] - 2T[n,i] + T[n,i-1])
end
end
#plot
using PyPlot
figure()
c=contour(t/(1e6*24*3600*365), -x/1e3, permutedims(T), levels=0:200:1200)
c["clabel"](fmt="%1.0f")
xlabel("time (My)")
⌃ ⇧
ylabel("depth (km)")
12 LIST OF CODES
Code 10: Matlab implementation of explicit FD method to solve the plate cooling problem
with heterogeneous heat diffusivity.
⌥ ⌅
# parameters
L = 100e3
↵0 = 2e-6
h = 50e3
Tm = 1200
Tt = 0
tmax = 50e6*365*24*3600
diffth(a,x,x0) =a*exp(-x/x0)
# discretisation
x = 1e3
tmax = 0.5* xˆ2/↵0
t = tmax/2
t = 0: t:tmax
x = 0: x:L
N = size(t,1)
I = size(x,1)
# solution
T = zeros(N,I)
T[1,:] .= Tm
T[1,1] = Tt
for n=1:N-1
T[n+1, 1] = Tt
T[n+1, I] = Tm
for i=2:I-1
↵plus = diffth(↵0 , x[i]+0.5* x, h)
↵minus = diffth(↵0 , x[i]-0.5* x, h)
T[n+1,i] = T[n,i] + ( t/ xˆ2)*(↵plus*(T[n,i+1] - T[n,i]) -
↵minus*(T[n,i] - T[n,i-1]))
end
end
#plot
using PyPlot
figure()
c=contour(t/(1e6*24*3600*365), -x/1e3, permutedims(T), levels=0:200:1200)
c["clabel"](fmt="%1.0f")
xlabel("time (My)")
⌃ ⇧
ylabel("depth (km)")
LIST OF CODES 13
Code 11: Solution of the wave equation using finite differences.
⌥ ⌅
using PyPlot
function matrixB(s2)
B = zeros(I,I)
B[1,1] = 2*(1-s2)
B[1,2] = s2
for n=2:I-1
B[n,n-1] = s2
B[n,n] = 2*(1-s2)
B[n,n+1] = s2
end
B[I,I-1] = s2
B[I,I] = 2*(1-s2)
return B
end
function timeloop!(u, t, B, s2, a, b, I, N)

c = zeros(I)
for n=2:N-1
c[1] = s2*a(t[n])
c[I] = s2*b(t[n])
u[n+1] = B*u[n] - u[n-1] + c
end
return nothing
end
# parameters
v = 1.0
L = 1.0
# boundary conditions
a(t) = 0.0
b(t) = 0.0
# initial conditions
initial_u(x, x0 , w) = exp(-(x-x0 )ˆ2/(wˆ2))
initial_udot(x) = 0.0
# discretisation
I = 301
x = range(0,L,length=I+2)
x = x[2]-x[1]
xint = x[2:end-1]
# CFL criterion
tmax = x/v
t = 0.5* tmax
N = 4000
t = t*(0:N-1)
# Matrix B
s2 = (v* t/ x)ˆ2
B = matrixB(s2)
# initialisation
u = Array{Array{Float64,1}}(undef,N)
u[1] = initial_u.(xint, 0.5, 0.02)
c = zeros(I)
c[1] = s2*a(t[1])/2
c[I] = s2*b(t[1])/2
g = initial_udot.(xint)
u[2] = 0.5*B*u[1] + t*g + c
# compute solution:
timeloop!(u, t, B, s2, a, b, I, N)
#plot solution
figure();
for k=1:10:N
cla()
plot(xint, u[k]);
ylim(-1,1)
yield()
⌃ ⇧
end
14 LIST OF CODES
Code 12: Matlab implementation of implicit FD method.
⌥ ⌅
# parameters
↵ = 1e-6
⇢c = 2.6e6
q0 = 30e6
t0 = 0.05
Tb = 0
L = 0.01
tmax = 1
# flux
q(t,t0,q0) = q0*exp(-t/t0)
# discretisation
I = 200+1
x = range(0,stop=L,length=I)
x = x[2]-x[1]
N = 200
t = range(0,stop=tmax,length=N)
t = t[2]-t[1]
s = ↵* t/ xˆ2
# initialise T
T = zeros(I+1, N)
# make matrix A
A = zeros(I+1, I+1)
A[1,1] = -1
A[1,3] = 1
for i=2:I
A[i,i-1] = -s
A[i,i] = 1+2s
A[i,i+1] = -s
end
A[I+1,I+1] = 1
# compute
for n=1:N-1
# assign vector B
B = T[:,n]
B[1] = - x/(↵*⇢c) * q(t[n+1],t0 ,q0 )
B[I+1] = Tb
#solve
T[:,n+1] = A\B
end
# delete solution at ghost node

T = T[2:I+1,:]
# closed form solution using constant flux q0

using SpecialFunctions
Tclosed(x, t, q0, alpha, rhoc) = q0/alpha/rhoc*(

sqrt(alpha*t/⇡)*exp(-xˆ2/(4*alpha*t))
-(x/2) *erfc(x/(2*sqrt(alpha*t))) )
# plots
using PyPlot
figure()
pcolor(x,t,permutedims(T))
xlabel("\$x\$ (m)")
ylabel("\$t\$ (s)")
colorbar()
#figure()
⌃ ⇧
#plot(t,T[1,:],"k.", t, Tclosed.(0,t,q0 ,↵,⇢c),"r")
LIST OF CODES 15
Code 13: Matlab implementation of Crank-Nicholson method for the heat flow problem.
⌥ ⌅
# parameters
↵ = 1e-6
⇢c = 2.6e6
q0 = 30e6
t0 = 0.05
Tb = 0
L = 0.01
tmax = 1
# flux
q(t,t0,q0) = q0*exp(-t/t0)
# discretisation
I = 200+1
x = range(0,stop=L,length=I)
x = x[2]-x[1]
N = 200
t = range(0,stop=tmax,length=N)
t = t[2]-t[1]
s = ↵* t/ xˆ2
# initialise T
T = zeros(I+1, N)
T[I+1,1] = Tb
# make matrix P
P = zeros(I+1, I+1)
P[1,1] = -1
P[1,3] = 1
for i=2:I
P[i,i-1] = -s/2
P[i,i] = 1+s
P[i,i+1] = -s/2
end
P[I+1,I+1] = 1
# compute
for n=1:N-1
# assign vector B
B = zeros(I+1)
B[1] = T[1,n] - T[3,n] - x/(↵*⇢c)*(
q(t[n+1],t0 ,q0 ) + q(t[n],t0 ,q0 ))
for i=2:I
B[i] = T[i-1,n]*s/2 + T[i,n]*(1-s) + T[i+1,n]*s/2
end
B[I+1] = Tb
#solve
T[:,n+1] = P\B
end
# delete solution at ghost node

T = T[2:I+1,:]
# plots
using PyPlot
figure()
pcolor(x,t,permutedims(T))
xlabel("\$x\$ (m)")
ylabel("\$t\$ (s)")
colorbar()
# closed form solution using constant flux q0

using SpecialFunctions
Tclosed(x, t, q0, alpha, rhoc) = q0/alpha/rhoc*(

sqrt(alpha*t/⇡)*exp(-xˆ2/(4*alpha*t))
-(x/2) *erfc(x/(2*sqrt(alpha*t))) )
#figure()
⌃ ⇧
#plot(t,T[1,:],"k.", t, Tclosed.(0,t,q0 ,↵,⇢c),"r")
16 LIST OF CODES
Code 14: Matlab script to produce a good quality figure from computations performed
in code 13.
⌥ ⌅
using PyPlot
figure()
subplot(2,1,1)
plot(t, T[1,:], "-", linewidth=1)
xlabel("\$t\$ (s)")
ylabel("\$T\$ (K)")
legend(["\$x=0\$";
"\$x=$(round(x[10]*1e3,digits=0))\$ mm";
"\$x=$(round(x[50]*1e3,digits=0))\$ mm";
"\$x=$(round(x[100]*1e3,digits=0))\$ mm"])
subplot(2,1,2)
plot(x*1e3, T[:,4], "-", linewidth=1)
xlim(0,5)
xlabel("\$x\$ (mm)")
ylabel("\$T\$ (K)")
legend(["\$t=$(round(t[4]*1e3,digits=0))\$ s";
"\$t=$(round(t[16]*1e3,digits=0))\$ s";
"\$t=$(round(t[32]*1e3,digits=0))\$ s";
⌃ ⇧
"\$t=$(round(t[64]*1e3,digits=0))\$ s"])
LIST OF CODES 17
Code 15: Implementation of nonconstant grid spacing to solve the heat equation with
oscillating imposed T at the top and constant T at some position x = L.
⌥ ⌅
# parameters
L = 6
↵ = 1e-6
day = 24*3600
tfinal = 30day
⌧ = 1day
A = 15
T0 = 10
Tsurf(t, ampl, t0, T0) = T0 + ampl*sin(2⇡*t/t0)
# discretise
x = exp10.(range(-2, stop=log10(L/4), length=25))
x = cumsum( x)
tmax = 0.5*minimum( x)ˆ2/↵
t = 0.5* tmax
t = 0: t:tfinal
J = size(x,1)
N = size(t,1)
T = zeros(N,J)
T[1,:] .= T0
T[1,1] = Tsurf(t[1], A, ⌧, T0 )
for n=1:N-1
T[n+1, 1] = Tsurf(t[n+1], A, ⌧, T0 )
T[n+1, J] = T0
for i=2:J-1
T[n+1,i] = T[n,i] + (2↵* t)/( x[i]+ x[i-1])*(
(T[n,i+1]-T[n,i])/ x[i] - (T[n,i]-T[n,i-1])/ x[i-1]
)
end
end
#plots
using PyPlot
pcolormesh(t, -x, permutedims(T))
xlabel("time (s)")
ylabel("depth (m)")
⌃ ⇧
colorbar()
18 LIST OF CODES
Code 16: Root of f (x) = e x2 x using bisection method.
⌥ ⌅
function bisection(func, a, b, tol)
c = (a+b)/2
while abs(func(c))>tol
if func(c)<0
a=c
else
b=c
end
c = (a+b)/2
end
return c
end
f(x) = exp(-xˆ2)-x
⌃ ⇧
x = bisection(f, 1, 0, 1e-8)
Code 17: Root of f (x) = e x2 x using Newton-Raphson method.
⌥ ⌅
function newtonraphson1d(func, funcp, x0, tol)
x = x0
while abs(func(x))>tol
x -= func(x)/funcp(x)
end
return x
end
f(x) = exp(-xˆ2) - x
fp(x) = -2*x*exp(-xˆ2) - 1
⌃ ⇧
x = newtonraphson1d(f, fp, 0, 1e-8)
Code 18: Finding the local minmum of x4 + 2y 2 + xy 6y + 5x using Newton-Raphson

method.
⌥ ⌅
function newtonraphson(func, jac, x0, tol)
x = copy(x0)
ff = func(x)
while maximum(abs.(ff))>tol
jj = jac(x)
x += jj\(-ff)
ff = func(x)
end
return x
end
p(x,y) = xˆ4 + 2*yˆ2 + x*y - 6*y + 5*x
F(x) = [4*x[1]ˆ3 + x[2] + 5

4*x[2] + x[1] - 6]
J(x) = [12*x[1]ˆ2 1
1 4]
⌃ ⇧
x = newtonraphson(F, J, [-2;2], 1e-8)

Notes Wave Combined

Uploaded by

Copyright:

Available Formats

Notes Wave Combined

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Notes Wave Combined

Uploaded by

Copyright:

Available Formats

Numerical methods for Earth Sciences

2 Ordinary Di↵erential Equations 11

3 Partial Di↵erential Equations 25

3.3 A note about normalisation . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

A more generic way to formulate this approximation is:

Example: f (x) = x cos(x). The exact answers are known:

I= 2, F (x) = 1 + cos(x) x sin(x).

1.3 Midpoint Method 中会1tn

Code 1.1: Integration of x cos(x) using the method of rectangles.

1.4 Trapezoidal rule

geomtrical considerations can be used to find:

Code 1.2: Integration of x cos(x) using the trapezoidal rule.

Ordinary Di↵erential Equations

y 0 = f (t, y), (+BC) (2.1)

2.1.2 Higher order ODEs

then Equation 2.2 is rewritten as

% initialise vector of solutions

Example: Let’s look at the following ODE:

slope = y 0 (t + h/2) = f (t + h/2, y(t + h/2)).

But we don’t know y(t + h/2) ! The idea is to approximate it as

y(t + h/2) ⇡ y(t) + (h/2)f (t, y(t)).

The algorithm is thus:

yi+1 = yi + hf ti + h/2, yi + (h/2)f (ti , yi ) .

Example: Let’s look at the following ODE, same as before:

2.4 Stability of a numerical method

• the type of di↵erential equation we are solving,

Code 2.2: Solution of y 0 = x o 2xy, y(0) = 1 using the midpoint method.

% generate vector of x values

% initialise vector of solutions

% compute solution using the midpoint method

• the numerical method,

• the step size

Stability can be assessed using the test case

2.5 Runge-Kutta methods

Several choices are possible:

Heun’s method: This is the case when a1 = a2 = 1/2, and

2.6 Systems of ODEs

function [dX] = f(X, rho,sigma, beta)

%solution using RK4

X(:,i+1) = X(:,i) + h*(k1/6 + k2/3 + k3/3 + k4/6);

%put solution matrix X into separate variables for simpler handling:

%now we can plot stuff

Code 2.5: Implementation of adptative time-stepping in the RK4 method.

% initial step size

% main computation loop

%solution using h (twice)

%take the difference

2. Neumann conditions (values given in terms of y 0 ):

and k is the thermal conductivity of rocks (here we will use k = 4 W m 1 K 1 ).

1. Constant temperatures at two points:

where in practice we will choose q/k = 36e 3Km 1 and TL = 800 C.

2.8.1 The finite di↵erence method: some generalities

where the corresponding spatial nodes are x = xi (i = 0, . . . , n), where x0 = 0 and xn = L.

and the vectors X and B are

and the vectors X and B are

2.8.4 Going further

Important: Remember that stability has nothing to do with accuracy, so a stable

3.1.3 Heterogeneities and nonlinearities

with numerical values ↵ = 10 6 m2 s 1 , Tm = 1200 C, Ttop = 0 C, and L = 100 km.

and initial condition

tmax = 50e6* 365243600; % in seconds !!

tmax = 50e6* 365243600; % in seconds !!

u(:,2) = 0.5Bu(:,1) + Dt*g + c;