The Finite Difference Method

Heiner Igel

Department of Earth and Environmental Sciences

Ludwig-Maximilians-University Munich

Heiner Igel Computational Seismology 1 / 32

Outline

1 Introduction
Motivation
History
Finite Differences in a Nutshell

2 Finite Differences and Taylor Series

Finite Difference Definition
Higher Derivatives
High-Order Operators

3 Finite-Difference Approximation of Wave Equations

Heiner Igel Computational Seismology 2 / 32

 Introduction Motivation

Motivation

Simple concept
Robust
Easy to parallelize
Regular grids
Explicit method

Heiner Igel Computational Seismology 3 / 32

 Introduction History

History

Several Pioneers of solving PDEs with finite-difference method

(Lewis Fry Richardson, Richard Southwell, Richard Courant, Kurt
Friedrichs, Hans Lewy, Peter Lax and John von Neumann)
First application to elastic wave propagation (Alterman and Karal,
1968)
Simulating Love waves and was the frst showing snapshots of
seimsic wave fields (Boore, 1970)
Concept of staggered-grids by solving the problem of rupture
propagation (Madariaga, 1976 and Virieux and Madariaga, 1982)

Heiner Igel Computational Seismology 4 / 32

 Introduction History

History

Several Pioneers of solving PDEs with finite-difference method

(Lewis Fry Richardson, Richard Southwell, Richard Courant, Kurt
Friedrichs, Hans Lewy, Peter Lax and John von Neumann)
First application to elastic wave propagation (Alterman and Karal,
1968)
Simulating Love waves and was the frst showing snapshots of
seimsic wave fields (Boore, 1970)
Concept of staggered-grids by solving the problem of rupture
propagation (Madariaga, 1976 and Virieux and Madariaga, 1982)

Heiner Igel Computational Seismology 4 / 32

 Introduction History

History

Several Pioneers of solving PDEs with finite-difference method

(Lewis Fry Richardson, Richard Southwell, Richard Courant, Kurt
Friedrichs, Hans Lewy, Peter Lax and John von Neumann)
First application to elastic wave propagation (Alterman and Karal,
1968)
Simulating Love waves and was the frst showing snapshots of
seimsic wave fields (Boore, 1970)
Concept of staggered-grids by solving the problem of rupture
propagation (Madariaga, 1976 and Virieux and Madariaga, 1982)

Heiner Igel Computational Seismology 4 / 32

 Introduction History

History

Several Pioneers of solving PDEs with finite-difference method

(Lewis Fry Richardson, Richard Southwell, Richard Courant, Kurt
Friedrichs, Hans Lewy, Peter Lax and John von Neumann)
First application to elastic wave propagation (Alterman and Karal,
1968)
Simulating Love waves and was the frst showing snapshots of
seimsic wave fields (Boore, 1970)
Concept of staggered-grids by solving the problem of rupture
propagation (Madariaga, 1976 and Virieux and Madariaga, 1982)

Heiner Igel Computational Seismology 4 / 32

 Introduction History

History

Extension to 3D because of parallel computations (Frankel and

Vidale, 1992; Olsen and Archuleta, 1996; etc.)
Application to spherical geometry by Igel and Weber, 1995;
Chaljub and Tarantola, 1997 and 3D spherical sections by Igel et
al., 2002
Incorporation in the first full waveform inversion schemes initially
in 2D, e.g. (Crase et al., 1990) and later in 3D (Chen et al., 2007)

Heiner Igel Computational Seismology 5 / 32

 Introduction History

History

Extension to 3D because of parallel computations (Frankel and

Vidale, 1992; Olsen and Archuleta, 1996; etc.)
Application to spherical geometry by Igel and Weber, 1995;
Chaljub and Tarantola, 1997 and 3D spherical sections by Igel et
al., 2002
Incorporation in the first full waveform inversion schemes initially
in 2D, e.g. (Crase et al., 1990) and later in 3D (Chen et al., 2007)

Heiner Igel Computational Seismology 5 / 32

 Introduction History

History

Extension to 3D because of parallel computations (Frankel and

Vidale, 1992; Olsen and Archuleta, 1996; etc.)
Application to spherical geometry by Igel and Weber, 1995;
Chaljub and Tarantola, 1997 and 3D spherical sections by Igel et
al., 2002
Incorporation in the first full waveform inversion schemes initially
in 2D, e.g. (Crase et al., 1990) and later in 3D (Chen et al., 2007)

Heiner Igel Computational Seismology 5 / 32

 Introduction Finite Differences in a Nutshell

Finite Differences in a Nutshell

Snapshot in space of the

pressure field p
Zoom into the wave field
with grid points indicated
by +
Exact interpolate using
Taylor series

Heiner Igel Computational Seismology 6 / 32

 Introduction Finite Differences in a Nutshell

1D acoustic wave equation

p̈(x, t) = c(x)2 ∂x2 p(x, t) + s(x, t)

p pressure
c acoustic velocity
s source term

Approximation with a difference formula

p(x, t + dt) − 2p(x, t) + p(x, t − dt)
p̈(x, t) ≈
dt 2

and equivalently for the space derivative

Heiner Igel Computational Seismology 7 / 32

 Finite Differences and Taylor Series Finite Difference Definition

Finite Differences and Taylor Series

Heiner Igel Computational Seismology 8 / 32

 Finite Differences and Taylor Series Finite Difference Definition

Finite Differences

Forward derivative
f (x + dx) − f (x)
dx f (x) = lim
dx → 0 dx

Centered derivative
f (x + dx) − f (x − dx)
dx f (x) = lim
dx → 0 2dx

Backward derivative
f (x) − f (x − dx)
dx f (x) = lim
dx → 0 dx

Heiner Igel Computational Seismology 9 / 32

 Finite Differences and Taylor Series Finite Difference Definition

Finite Differences

Forward derivative
f (x + dx) − f (x)
dx f + ≈
dx

Centered derivative
f (x + dx) − f (x − dx)
dx f c ≈
2dx

Backward derivative
f (x) − f (x − dx)
dx f − ≈
dx

Heiner Igel Computational Seismology 10 / 32

 Finite Differences and Taylor Series Finite Difference Definition

Finite Differences and Taylor Series

The approximate sign is important here as the derivatives at point

x are not exact. Understanding the accuracy by looking at the
definition of Taylor Series:

1
f (x + dx) = f (x) + f 0 (x) dx + 2! f 00 (x) dx 2 + O(dx 3 )

Subtraction with f(x) and division by dx leads to the definition of

the forward derivative:

f (x+dx)−f (x) 1
dx = f 0 (x) + 2! f 00 (x) dx + O(dx 2 )

Heiner Igel Computational Seismology 11 / 32

 Finite Differences and Taylor Series Finite Difference Definition

Finite Differences and Taylor Series

Using the same approach - adding the Taylor Series for f (x + dx)
and f (x − dx) and dividing by 2dx leads to:

f (x+dx)−f (x−dx)
2dx = f 0 (x) + O(dx 2 )

This implies a centered finite-difference scheme more rapidly

converges to the correct derivative on a regular grid
=⇒ It matters which of the approximate formula one chooses
=⇒ It does not imply that one or the other finite-difference
approximation is always the better one

Heiner Igel Computational Seismology 12 / 32

 Finite Differences and Taylor Series Higher Derivatives

Higher Derivatives

The partial differential equations have often 2nd (seldom higher)

derivatives
Developing from first derivatives by mixing a forward and a
backward definition yields to

f (x+dx)−f (x)
− f (x)−fdx(x−dx) f (x + dx) − 2f (x) + f (x − dx)
∂x2 f ≈ dx
=
dx dx 2

Heiner Igel Computational Seismology 13 / 32

 Finite Differences and Taylor Series Higher Derivatives

Higher Derivatives

Determining the weights with which the function values have to be

multiplied to obtain derivative approximations
The Taylor series for two grid points at x ± dx, include the function
at the central point as well and multiply each by a real number

1
af (x + dx) = a f (x) + f 0 (x) dx + f 00 (x) dx 2 + . . .
 
2!
bf (x) = b [f (x)]
1
cf (x − dx) = c f (x) − f 0 (x) dx + f 00 (x) dx 2 − . . .
 
2!

Heiner Igel Computational Seismology 14 / 32

 Finite Differences and Taylor Series Higher Derivatives

Higher Derivatives

Summing up the right-hand sides of the previous equations and

comparing coefficients leads to the following system of equations:
1
a=
a+b+c =0 dx 2

a−c =0 ⇒ b=− 2
2
dx
2! 1
a+c = c=
dx 2 dx 2

Heiner Igel Computational Seismology 15 / 32

 Finite Differences and Taylor Series High-Order Operators

High-Order Operators

What happens if we extend the domain of influence for the

derivative(s) of our function f(x)?
Let us search for a 5-point operator for the second derivative

f 00 ≈ af (x + 2dx) + bf (x + dx) + cf (x) + df (x − dx) + ef (x − 2dx)

a + b + c +d +e = 0
2a + b − d − 2e = 0
1
4a + b + d + 4e =
2dx 2
8a + b − d − 8e = 0
16a + b + d + 16e = 0

Heiner Igel Computational Seismology 16 / 32

 Finite Differences and Taylor Series High-Order Operators

High-Order Operators

Using matrix inversion we obtain a unique solution

1
a = −
12dx 2
4
b =
3dx 2
5
c = −
2dx 2
4
d =
3dx 2
1
e = − .
12dx 2

with a leading error term for the 2nd derivative is O(dx 4 )

=⇒ Accuracy improvement

Heiner Igel Computational Seismology 17 / 32

 Finite Differences and Taylor Series High-Order Operators

High-Order Operators

Graphical illustration of the

Taylor Operators for the
first derivative for higher
orders
The weights rapidly
become small for
increasing distance to
central point of evaluation

Heiner Igel Computational Seismology 18 / 32

 Finite-Difference Approximation of Wave Equations

Finite-Difference Approximation of
Wave Equations

Heiner Igel Computational Seismology 19 / 32

 Finite-Difference Approximation of Wave Equations

Acoustic waves in 1D

To solve the wave equation, we start with the

simplemost wave equation:
The constant density acoustic wave equation
in 1D

p̈ = c 2 ∂x2 p + s

impossing pressure-free conditions at the two

boundaries as

p(x) |x=0,L = 0

Heiner Igel Computational Seismology 20 / 32

 Finite-Difference Approximation of Wave Equations

Acoustic waves in 1D

The following dependencies apply:

p → p(x, t) pressure
c → c(x) P-velocity
s → s(x, t) source term

As a first step we need to discretize space and time and we do that

with a constant increment that we denote dx and dt.
xj = jdx, j = 0, jmax
tn = ndt, n = 0, nmax

Heiner Igel Computational Seismology 21 / 32

 Finite-Difference Approximation of Wave Equations

Acoustic waves in 1D
Starting from the continuous description of the partial differential
equation to a discrete description. The upper index will correspond to
the time discretization, the lower index will correspond to the spatial
discretization

pjn+1 → p(xj , tn + dt)

pjn → p(xj , tn )
pjn−1 → p(xj , tn − dt)
n
pj+1 → p(xj + dx, tn )
pjn → p(xj , tn )
n
pj−1 → p(xj − dx, tn )
.

Heiner Igel Computational Seismology 22 / 32

 Finite-Difference Approximation of Wave Equations

Acoustic waves in 1D

pjn+1 − 2pjn + pjn−1

" n − 2p n + p n
#
pj+1 j j−1
= cj2 + sjn .
dt 2 dx 2

the r.h.s. is defined at same n+1

time level n i
the l.h.s. requires information n−1
from three different time levels

j −1 j j +1

Heiner Igel Computational Seismology 23 / 32

 Finite-Difference Approximation of Wave Equations

Acoustic waves in 1D

Assuming that information at time level n (the presence) and n − 1 (the

past) is known, we can solve for the unknown field pjn+1 :

dt 2 h n i
pjn+1 = cj2 pj+1 − 2pj
n
+ p n
j−1 + 2pjn − pjn−1 + dt 2 sjn
dx 2
The initial condition of our wave simulation problem is such that
everything is at rest at time t = 0:

p(x, t)|t=0 = 0 , ṗ(x, t)|t=0 = 0.

Heiner Igel Computational Seismology 24 / 32

 Finite-Difference Approximation of Wave Equations

Acoustic waves in 1D

Waves begin to radiate as soon as the source term s(x, t) starts to

act
For simplicity: the source acts directly at a grid point with index js
Temporal behaviour of the source can be calculated by Green’s
function

s(x, t) = δ(x − xs ) δ(t − ts )

where xs and ts are source location and source time and δ()
corresponds to the delta function

Heiner Igel Computational Seismology 25 / 32

 Finite-Difference Approximation of Wave Equations

Acoustic waves in 1D

A delta function contains all frequencies and we cannot expect that our
numerical algorithm is capable of providing accurate solutions
Operating with a band-limited source-time function:

s(x, t) = δ(x − xs ) f (t)

where the temporal behaviour f(t) is chosen according to our specific

physical problem

Heiner Igel Computational Seismology 26 / 32

 Finite-Difference Approximation of Wave Equations

Example

Simulating acoustic wave propagation in a 10km column (e.g. the

atmosphere) and assume an air sound speed of c = 0.343m/s. We
would like to hear the sound wave so it would need a dominant
frequency of at least 20 Hz. For the purpose of this exercise we
initialize the source time function f (t) using the first derivative of a
Gauss function.
1
− (t−t0 )2
(4f0 )2
f (t) = −8 f0 (t − t0 ) e

where t0 corresponds to the time of the zero-crossing, f0 is the

dominant frequency

Heiner Igel Computational Seismology 27 / 32

 Finite-Difference Approximation of Wave Equations

Example

What is the minimum spatial

wavelength that propagates inside the
medium?
What is the maximum velocity inside
the medium?
What is the propagation distance of
the wavefield (e.g., in dominant
wavelengths)?

Heiner Igel Computational Seismology 28 / 32

 Finite-Difference Approximation of Wave Equations

Example

Sufficient to look at the relation between frequency and wavenumber:

ω λ
c = = = λf
k T
where c is velocity, T is period, λ is wavelength, f is frequency, and
ω = 2πf is angular frequency
dominant wavelength of f0 = 20Hz
substantial amount of energy in the wavelet is at frequencies
above 20 Hz
=⇒ λ = 17m and λ = 7m for frequencies 20Hz and 50Hz,
respectively

Heiner Igel Computational Seismology 29 / 32

 Finite-Difference Approximation of Wave Equations

Matlab code fragment

% Time extrapolation
for it = 1 : nt,
(...)
% Space derivatives
for j=2:nx-1
d2p(j)=(p(j+1)-2*p(j)+p(j-1))/dxˆ 2;
end
% Extrapolation
pnew=2*p-pold+cˆ 2*d2p*dtˆ 2;
% Source input
pnew(isrc)=pnew(isrc)+src(it)*dtˆ 2;
% Time levels
pold=p;
p=pnew;
(...)

Heiner Igel Computational Seismology 30 / 32

 Finite-Difference Approximation of Wave Equations

Result

Choosing a grid increment of

dx = 0.5m −→ about 24 points per
spatial wavelength for the dominant
frequency
Setting time increment dt = 0.0012 −→
around 40 points per dominant period

Heiner Igel Computational Seismology 31 / 32

 Finite-Difference Approximation of Wave Equations

Summary

Replacing the partial derivatives by finite differences allows partial

differential equations such as the wave equation to be solved directly for
(in principle) arbitrarily heterogeneous media
The accuracy of finite-difference operators can be improved by using
information from more grid points (i.e., longer operators). The weights
for the grid points can be obtained using Taylor series

Heiner Igel Computational Seismology 32 / 32