Transcript Numerical Methods for Partial Differential Equations

Numerical Methods for Partial
Differential Equations
CAAM 452
Spring 2005
Lecture 7
Instructor: Tim Warburton
Summary of Small Theta Analysis
• The dominant remainder term in this analysis
relates to a commonly used, physically motivated
description of the shortfall of the method:
1 2 t
2 m
3 3
i
ct   mt  O m  
du
Dissipative
L
 c u  uˆm  uˆm  0  e
e dx dx
dt
1 2 t
2 m
3
c

t

O

i
ct

 Unstable
m
m

du
2!dx
dx
L
 c u  uˆm  uˆm  0  e
e
dt
 i 3 it
2 m
c
 mt  O m5 
Dispersive
i
ct

du
3!dx
dx
L
 c 0u  uˆm  uˆm  0  e
e
dt
 i 4 5 it
2 m
c
 mt  O m7 
i
c
t

du
Dispersive
 c 4u  uˆm  uˆm  0  e L e dx 5! dx
dt
CAAM 452 Spring 2005
Dissipation in Action
• Consider the right difference version:
du
 c u
dt
 uˆm  uˆm  0  e
i
2 m
ct
L
e
1 2 t
 mt  O
dx
dx
 3m  
3
• We are going to experimentally determine how
much dissipation the solution experiences.
CAAM 452 Spring 2005
Testing Methodology
1)
Reduce the time-stepping error to a secondary effect by choosing a 4th
order Runge-Kutta (JST) time-integrator and a small dt.
2)
Fix the method (choose one of the difference operators for the spatial
derivative)
3)
Do a parameter study in M, i.e. we ask the questions: how does
increasing the number of data points change the error.
4)
We need to understand what questions we are asking:
1) Is the computer code stable (as predicted by the theory)?
2) Does the computer code converge (as predicted by the theory)?
3) What order of accuracy in M do we achieve?
4) We hypothesize that if the theory holds then we should achieve:
u T , xm   u T , xm   Cdx p
5)
What is the actual approximate p achieved?
CAAM 452 Spring 2005
Initial Condition
• Because we do not wish to introduce uncertainty
over the source of errors in the computation we use
an initial condition which is infinitely smooth.
u  x,0   e
4cos 2  x 
with x [0,2)
CAAM 452 Spring 2005
Computing Approximate
Order of Accuracy
• We compute the error by:
 : max u  xm , t  8   u  xm , t  8 
m
• We will compute the approximate order of accuracy by
p
assuming:
L
  Cdx p  C  
M 
• where C is independent of M
• If we measure the error for two different M then we can
p
eliminate the constant:
 L  
 
1  C 

log  1 
p

 M 1   1  M 2 
 2 



p

M 
p

 M2 
 1
2
 L  
log
M 
2  C 


 1
 M2  
• In the case when M2=2M1
 
log  1 
 2 
p
log  2 
CAAM 452 Spring 2005
After 4 Periods
du
 c u
dt
• The numerical pulses are in the right place but
have severely diminished amplitude.
CAAM 452 Spring 2005
After 4 Periods
du
 c u
dt
M
Maximum error at t=8
20
0.69149169495179
40
0.69137534363447
0
80
0.68329408124856
0.01
160 0.63163960492185
0.11
320
0.26
0.52742560650408
~order
After 4 periods the solution is totally flattened in all but the last 2 results. If we bothered
to keep increasing M we would eventually see the error decline as 1/M
CAAM 452 Spring 2005
The Unstable Left Stencil
du
 c u
dt
• I repeated this with everything the same, except the choice
of   instead of  
 1  2 m 2 t

3 
2 m
c
t

O







m
i
ct
du
2!
d
x
M
dx





 c u  uˆm  uˆm  0  e L e
dt
We clearly see that there is initial growth
near the pulses, but eventually the
dominant feature is the highly oscillatory
and explosive growth
(large m in the above red term).
CAAM 452 Spring 2005
Cont (snapshot)
du
 c u
dt
 uˆm  uˆm  0  e
 1  2 m 2 t
3
2 m
i
ct c  2!dx  M  t  dx O  m

L
e

 

CAAM 452 Spring 2005
Dispersion in Action
du
 c 0u
dt
• Consider the central difference version:
du
 c 0u
dt
 uˆm  uˆm  0  e
i

e
 
 i 3 it
2 m
 mt  O  m5
ct c
3!dx
dx
L
• With the same time integrator before, M=100,…
• We note that the remainder terms are dispersive corrections
– i.e. they indicate that modes of different frequency will
travel with different speed.
• Furthermore, to leading order accuracy the higher order
modes will travel more and more slowly as m increases.
CAAM 452 Spring 2005
2nd Order Central Difference
After 4 Periods
du
 c 0u
dt
• We notice that the humps start to shed rear oscillations as
the higher frequency Fourier components lag behind the
lower frequency Fourier components.
CAAM 452 Spring 2005
Convergence Study (t=8)
du
 c 0u
dt
What should we use as an error
Indicator ??
CAAM 452 Spring 2005
2nd
Order Central Difference
du
 c 0u
dt
M
Maximum error at t=8
20
0.94240801474506
40
0.35029899850733
1.42
80
0.38298512776859
0.13
160 0.19118127466878
1.00
320
1.80
0.05471730524482
~order
• We do not see a convincing 2nd order accuracy
• I computed this by log2(error_M/error{M-1})
CAAM 452 Spring 2005
4th
Order Central Differencing
du
 c 4u  uˆm  uˆm  0  e
dt
i

e
du
 c 4u
dt
 
 i 4 5 it
2 m
 mt  O  m7
ct c
dx 5!
dx
L
CAAM 452 Spring 2005
4th
Order Central Difference
du
 c 4u
dt
M
Maximum error at t=8
~order
20
0.44353878542277
40
0.12852952972116
1.79
80
0.02040190886767
2.66
160 0.00134076478259
3.93
320 0.00008319438414
4.01
• We see pretty convincing 4th order accuracy
• I computed this by log2(error_M/error{M-1})
CAAM 452 Spring 2005
6th
Order Central Difference
du
 c 6u
dt
M
Maximum error at t=8
~order
20
0.19409575615520
40
0.04953652526983
1.97
80
0.00169299155603
4.87
160 0.00002891724105
5.87
320 0.00000046041415
5.97
• We see pretty convincing 6th order accuracy
• I computed this by log2(error_M/error{M-1})
CAAM 452 Spring 2005
Summary of Testing Procedure
1)
2)
3)
4)
5)
6)
Understand what you want to test
Keep as many parameters fixed as possible
If possible, perform an analysis before hand
Run parameter sweeps to determine if code agrees with analysis.
Estimate error scaling with a single parameter if possible.
It should be apparent here that the errors computed in the simulations
only approximate the tidy power law (dx^p) in the limit of small dx
(large M).
7) It is quite common to refer to the range of M before the error scaling
applied as the “preasymptotic convergence range”.
8) The preasymptotic range is due to the other factors in the error
estimate which are much larger than the dx^p term (i.e. the p’th
derivatives may be very large, or the constant may be large…)
9) In this case we heuristically set dt small, using a high-order RungeKutta time integrator, and then performed a parameter sweep on M.
We could have been unlucky and the time-stepping errors may have
been dominant which would mask the dx behavior.
10) To be careful we would try this experiment with even smaller dt and
check if the scalings still hold.
CAAM 452 Spring 2005
Summary of our Approach to Designing
Finite Difference Methods
• We have systematically created finite difference methods by
separating the treatment of space and time derivatives.
• Then designing a solver consists of choosing/testing:
– a time integrator (so far we covered Euler-Forward, LeapFrog,
Adams-Bashford(2,3,4), Runge-Kutta)
– a discretization for spatial derivatives
– a discrete differential operator which has all eigenvalues in the
left half of the complex plane (assuming the PDE only admits
non-growing solutions).
– Possibly using Gerschgorin’s theorem to localize the
eigenvalues of the discrete differentiation operator
– dt so that dt*largest eigenvalue (by magnitude) of the
derivative operator sits inside the stability region  stability.
– (small dx) spatial truncation analysis  consistency and
accuracy.
– Fourier analysis for classification of the differential operator.
– Writing code and testing.
CAAM 452 Spring 2005
Some Classical Finite Difference Methods
• However, there are a number of classical methods
which we have not discussed and do not quite fit
into this framework.
• We include these for completeness..
CAAM 452 Spring 2005
Leap Frog (space and time)
• We use centered differencing for both space and time.
 umn 1  umn 1 
umn1  umn1
 c

2dt
2
dx


• We know that the leap frog time stepping method is only
stable for operators with eigenvalues in the range:
dt  i  1,1
• However, we also know that the centered difference
derivative matrix is a skew symmetric matrix with
eigenvalues: ic
2 m


sin 

dx  M 
dx
cdt
• So we are left with a condition:
 1 i.e. dt  c
dx
CAAM 452 Spring 2005
cont
• We can perform a full truncation analysis (in space
and time):
u  xm , tn1   u  xm , tn1 
 u  xm1 , tn   u  xm1 , tn  
T 
 cdt 

2
2dx


cdt
3
 O  dt  
O  dx3 
dx
n
m
• We know that dt <= dx so
Tmn  O  dt 3   cO  dx3 
CAAM 452 Spring 2005
Lax-Friedrichs Method
• We immediately conclude that the following method is not
stable: u n1  u n
m
m
dt
 c 0umn
• Because the Euler-Forward time integrator only admits one
stable point (the origin) on the imaginary axis, but the
central differencing matrix has all eigenvalues on the
imaginary axis.
• However, we can stabilize this formula by replacing the
second term in the time-stepping formula:
umn1  0.5  umn 1  umn 1 
dt
 c 0umn
CAAM 452 Spring 2005
cont
• This formula does not quite fit into our constructive
process (method of lines approach).
umn1  0.5  umn 1  umn 1 
 c 0umn
dt
• We have admitted spatial averaging into the
discretization of the time derivative.
• We can rearrange this:
1 n
dt n
n 1
n
um   um1  um1   c
um1  umn 1 

2
2dx
1
dt  n
1
dt  n
 1  c  um1  1  c  um1
2
dx 
2
dx 
CAAM 452 Spring 2005
cont
u
n 1
m
1
dt  n
1
dt  n
 1  c  um1  1  c  um1
2
dx 
2
dx 
• We can immediately determine that this is a stable method as
long as c*dt/dx <=1
• Given, this condition we observe that the time updated
solution is always bounded between the values of the left and
right neighbor at the previous time because this is an
interpolation formula.
CAAM 452 Spring 2005
cont
• A formal stability analysis would involve:
1
dt  n
1
dt  n
u  1  c  um1  1  c  um1
2
dx 
2
dx 
1
dt 
1
dt 
 u n1  1  c  Ru n  1  c  Lu n
2
dx 
2
dx 
1
dt
1
dt
 umn1  1  c  eim umn  1  c  e  im umn
2
dx 
2
dx 
dt
 cos  m   c i sin  m  umn
dx
n 1
m


• Which gives stability for each mode if:
dt
cos  m   c i sin  m   1
dx
CAAM 452 Spring 2005
cont
2 m
• Thus considering all the possible modes  m 
M
• We note that the middle mode requires:
  dt



cos    c i sin    1
dx
2
2
dt
 c 1
dx
• This condition gives a sufficient condition for all modes to be
bounded.
• By the invertibility and boundedness of the Fourier transform
we conclude that the original equation is stable.
CAAM 452 Spring 2005
cont
• We can recast the Lax-Friedrichs method again
1
dt  n
1
dt  n
n 1
um  1  c  um1  1  c  um1
2
dx 
2
dx 

u
n 1
m
1 2 2 n
 1  cdt 0  u  dx  um
2
n
m
• This method consists of Euler Forward, central
differencing for the space derivative and an extra
dissipative term (i.e. a grid dependent advection
diffusion equation: u u dx 2  2u )
t
c
x

2 x 2
CAAM 452 Spring 2005
CFL Number
• The ratio  
dt
appears repeatedly, in particular in
dx
the estimates for the maximum possible time step.
• We refer to this quantity as the CFL (CourantFriedrichs-Lewy) number.
dt

 C
• Bounds of the form:
which result from a
dx
stability analysis are frequently referred to as CFL
conditions.
CAAM 452 Spring 2005
Lax-Wendroff Method
• We are fairly free to choose the parameter in the stabilizing
term:
u
n 1
m
1 2 2 n
 1  cdt 0  u  cdt  um
2
n
m
• The artificial viscosity term acts to shift the eigenvalues of
the spatial operator into the left half-plane.
• Recall – the Euler-Forward stability region is the unit circle
centered at -1 in the complex plane. So pushing the
eigenvalues off the imaginary axis allows us to choose a dt
small enough to push the eigenvalues of the discrete spatial
operator into the stability region…
CAAM 452 Spring 2005
Class Exercise
• Please provide a CFL condition for:
u
n1
m
u
n 1
m
4
1
n
n
n
n
 c  um1  um1   c  um2  um2 
3
6
• In terms of:  
dt
dx
• What is the truncation error?
CAAM 452 Spring 2005
Von Neumann Analysis
• By stealth I have introduced the classical Von
Neumann analysis.
• The first step is to Fourier transform the finite
difference equation in space.
• A short cut is to make the following substitutions:
u u
n
m
n
m
umn 1  eim umn
u
n
m 1
e
 i m
umn
umn  2  e 2im umn
umn 2  e 2im umn
......
CAAM 452 Spring 2005
Analysis cont
• We can also make the substitutions for the
difference operators:

i m

1
1
 m 
i m
2
 
1  e
 2ie sin   
dx
dx 
 2 




i m

1
1
 m 

 i m
2
i m n
n
1 e
 2ie sin   
um1  e um    
dx
dx 
 2 
umn 1  e  im umn 
   1 eim  e  im  1 i sin  

0
m
umn 2  e 2im umn 
2dx
dx

1 
umn 2  e 2im umn 
2
 
2  m 
     2 4sin   
dx 

......
 2 
2
 2 1  cos  m 
dx
CAAM 452 Spring 2005
umn  umn




Example
• The Lax-Wendroff scheme
u
n 1
m
1 2 2 n
 1  cdt 0  u  cdt  um
2
n
m
• becomes
i m
 i m
2i m
2i m





 n
e

e
1
e

2

e
n 1
n
2
um  1  cdt 
um  cdt 
um

2


2
dx
 2dx  



CAAM 452 Spring 2005
Example cont
• So the Lax-Wendroff scheme becomes a neat one
step method for each Fourier coefficient:
i m
 i m
2i m
2i m





 n
e

e
1
e

2

e
n 1
n
2
um  1  cdt 
um  cdt 
um

2


2
dx
 2dx  



• We have neatly decoupled the Fourier coefficients
so we are back to solving a recurrence relation in n
n
for each um:
m
• Hence we need to examine the roots of the stability
polynomial for the above

 eim  e  im   1 2  e 2im  2  e 2im  
z  1  cdt 
 cdt 


2


dx
 2dx   2



CAAM 452 Spring 2005
Example cont
• This case is trivial as we just need to bound the
single root by 1

 eim  e  im   1 2  e 2im  2  e 2im  
z  1  cdt 
 cdt 


2


dx
 2dx   2



 1  ci sin  m   c 2  cos  2   1
• Plot it this as a function of theta…
CAAM 452 Spring 2005
Final Von Neuman Shortcut
• We can skip lots of steps by making the direct
substitution:
n
n im
um  g e
• Here g is referred to as the amplification factor and
i*m*theta is our previous theta_m
CAAM 452 Spring 2005
Next Time
• Definition of consistency
• Lax-Richtmyer equivalence theory
• Treating higher order derivatives
CAAM 452 Spring 2005
Homework 3
u
u
c

t

x
Q1a) use your Cash-Karp Runge-Kutta time integrator from HW2 for time stepping
Q1) Build a finite-difference solver for
Q1b) use the 4th order central difference in space (periodic domain)
Q1c) perform a stability analysis for the time-stepping (based on visual inspection of
the CK R-K stability region containing the imaginary axis)
Q1d) bound the spectral radius of the spatial operator
Q1e) choose a dt well in the stability region
u  x,0   e
4cos 2  x 
Q1f) perform four runs with initial condition
(use M=20,40,80,160) and compute maximum error at t=8
, x   0,2
Q1g) estimate the accuracy order of the solution.
Q1h) extra credit: perform adaptive time-stepping to keep the local truncation error
from time stepping bounded by a tolerance.
CAAM 452 Spring 2005