An Introduction to Numerical Methods

1
An Introduction to Numerical Methods

• For Weather Prediction
• by
• Joerg Urban
• office 012
• Based on lectures given by Mariano Hortal

2
Shallow water equations in 1 dimension
Non linear equations
h
advection
adjustement
diffusion
u velocity along x direction h absolute
height g acceleration due to gravity K
diffusion coefficient
Velocity is equal in layers vertical direction
(shallow)
3
Linearization
uU0u hH h
Const. perturb. in the x-comp. of
velocity Const. perturb. in the height of the
free surface
Small perturbations
Substitute and drop products of perturbations
4
Classification of PDEs

• Boundary value problems

D
open domain
GD
f, g known function L, B differential
operator f unknown function of x
boundary

• Initial value problems (most important to us)

f unknown function of t
t
t0
5
Classification of PDEs (II)

• Initial and boundary value problems
• Eigenvalue problems

f unknown eigenfunction ? eigenvalue of
operator L
6
Existence and uniqueness of solutions
y(t0)y0 (initial value problem)

• Does it have a solution?
• Does it have only one solution?
• Do we care?

If it has one and only one solution it is called
a well posed problem
7
Picards Theorem
Let
and
be continuous in the rectangle
then, the initial-value problem
Has a unique solution y(t) on the interval
Finding the solution (not analytical) Numerical
methods (finite dimensions)
8
Discretization

• Finite differences
• Spectral
• Finite elements

Transform the continuous differential
equation into a system of ordinary algebraic
equations where the unknowns are the numbers fj
9
The Lax-Richtmeyer theorem

• Convergence
• Consistency
• Stability
• Lax-Richtmeyer theorem

Discretized solution ---------gt continuous
solution discretization finer and finer
Discretized equation ---------gt continuous
equation
Discretized solution bounded
If a discretization scheme is consistent and
stable then it is convergent, and vice versa
10
Finite Differences - Introduction
1
j-1
j
j1
2
N
N1
lt------------------------------- L
----------------------------------gt
Taylor series expansion
11
Finite differences approximations
forward approximation
Consistent if are bounded
backward approximation
adding both
centered differences
Consistent if are bounded
12
Finite differences approximations (2)
Also
fourth order approximation to the first derivative
Using the Taylor expansion again we can easily
get the second derivative
second order approximation of the second
derivative
13
The linear advection equation
initial and boundary conditions
We start with a guess
Substituting we get
Eigenvalue problems for
With periodic B.C. ? can only have certain
(imaginary) values where k is the wave number
The general solution is a linear combination of
several wave numbers
14
The linear advection equation (2)
The analytic solution is then
Propagating with speed U0
No dispersion
For a single wave of wave number k the frequency
is ?kU0
Energy
For periodic B.C.
15
Space discretization
centered second-order approximation
Try
results in
whose solution is
with
c
The phase speed c depends on k
dispersion
U0
k?x p ---gt ?2?x gt c0
k?x
p
16
Group velocity
Continuous equation
Discretized equation
-U0 for k?xp
Approximating the space operator introduces
dispersion
17
Time discretization
In addition to our 2nd order centered approx. for
the space derivative we use a 1st order forward
approx. for the time derivative
Try
Courant-Friedrich- Levy number
Substituting we get
?aib
Also another dispersion is introduced, as we have
approximated the operator ?/ ?t
18
Three time level scheme (leapfrog)
This scheme is centered (second order accurate)
in both space and time
exponential
Try a solution of the form
If ?k gt 1 solution unstable if ?k 1
solution neutral if ?k lt 1 solution damped
Substituting
physical mode
?x---gt0 ?t ---gt0
computational mode
19
Stability analysis
Energy method
We have defined
We have discretized t and hence the discretized
analog of E(t) is En
For periodic boundary conditions
fn N1 fn1
If Enconst,
than the scheme is stable
20
Example of the energy method
upwind if U0gt0 downwind if U0lt0
a0 gt U00 no motion a1
?t ?x/U0
En1En if
En1 gt En unstable
a gt 0 gt U0 gt 0 upwind a lt 1
U0 ?t/ ?x lt 1 CFL cond.
damped
En1 lt En if
21
Von Neumann method
Consider a single wave
if ?k lt 1 the scheme is damping for this
wave number k if ?k 1 ?k the scheme
is neutral if ?k gt 1 for some value of k,
the scheme is unstable
alternatively
if Im(?) gt 0 scheme
unstable if Im(?) 0 scheme
neutral if Im(?) lt 0 scheme
damping
Vf ?/k vg??/?k
22
Stability of some schemes

• Forward in time, centered in space (FTCS) scheme
• Upwind or downwind

using Von Neumann, we find
scheme unstable
upwind if U0 gt 0 downwind if U0 lt 0
Using Von Neumann, we find
a lt 0 downwind a gt 1 CFL limit
a(a-1) gt 0 unstable
-1/4 lt a(a-1) lt 0 gt 0 a 1 stable
damped scheme
23
Stability of some schemes (cont)

• Leapfrog

Using von Neumann we find a1 as stability
condition
As a reminder a is the Courant-Friedrich-Levy
number
24
Stability of some schemes (cont)

• Lax Wendroff scheme

From a Taylor expansion in t we get
Substitution the advection equation we
Discretization
Applying Von Neumann we can find that a 1
-----gt stable
25
Stability of some schemes (cont)

• Implicit centered scheme

We replace the space derivation by the average
value of the centred space derivation at time
level n-1 and n1
using von Neumann
Always neutral, however an Expensive implicit
equation need to solved
Dispersion worse than leapfrog
26
Intuitive look at stability
If the information for the future time step
comes from inside the interval used for the
computation of the space derivative, the scheme
is stable. Otherwise it is unstable
x, x point where the information
comes from (xj-U0?t) Interval used
for the computation of ?f/?x
Downwind scheme (unstable)
U0?t
x
Upwind scheme (conditionally stable)
x if a lt 1 x if a gt 1
x x
CFL number gt fraction of ?x traveled in ?t
seconds
Leapfrog (conditionally stable)
x x
Implicit (unconditionally stable)
27
Dispersion and group velocity
Vf ?/k vg??/?k
Leapfrog
vg
vf
U0
K-N
??t
Implicit
p/2
p
28
Effect of dispersion
Initial
Leapfrog
implicit
29
Two-dimensional advection equation
Using von Neumann, assuming a solution of the form
we obtain
using
we obtain, for ? 1 the condition
where ?s ?x ?y
This is more restrictive than in one dimension by
a factor
30
Non linear advection equation Continuous form
u(x,t)
u(x,0)
Change in shape even for the continuous form
One Fourier component ukeikx no longer moving
with constant speed but interacting with other
components
Fourier decomposition valid at each individual
time but it changes amplitude with time
No analytical solution!
31
Energy conservation
Define again

Discretization in space
periodic B.C.
First attempt
Second attempt
terms joined by arrows cancel from consecutive
js
32
Aliassing
Aliasing occurs when the non-linear interactions
in the advection term produce a wave which is
too short to be represented on the grid.
Consider the product
in the interval 0x 2p
Minimum wavelength
1
2
n
N1
Maximum wave number representable with the
discretized grid
lt---------------L-------------gt
33
Aliassing (cont.)
Trigonometrical manipulations lead to
sin(kxj)-sin(2kM-k)xj
Therefore, it is not possible to distinguish wave
numbers k and (2kM-k) on the grid.
wave number k ? wave number 2kM-k
x
x
x
34
Non-linear instability
If k1k2 is misrepresented as k1 there is
positive feedback, which causes instability
k1 2kM - (k1k2) ----------gt 2k12kM-k2
2kM ?2k1 ? kM 2?x ?1 4?x
These wavelengths keep storing energy and total
energy is not conserved
We can remove energy from the smallest
wavelengths by - Fourier filtering -
Smoothing - Diffusion - Use some other
discretization (e.g. semi-Lagrangian)