Transcript Numerical Methods for Partial Differential Equations

Numerical Methods for Partial
Differential Equations
CAAM 452
Spring 2005
Lecture 6
Various Finite Difference Discretizations for the
Advection Equations – Phase and Dissipation Error
Instructor: Tim Warburton
Note on AB3
Stability
z  rei

z3  z 2
 :
(23 z 2 -16 z  5) /12
CAAM 452 Spring 2005
Review: After Time Stepping… Back to PDE’s
• We have now proposed two quite distinct but
reliable ways to advance an ODE approximately in
time.
• But our original goal was to solve a PDE
u
u
c
t
x
• And not just for one mode.
CAAM 452 Spring 2005
Review: Sample Points in Time Space
u
u
c
t
x
• We chose to sample the ODE at discrete points on the time axis.
• We can also make the same choice for the spatial representation
at each discrete time level.
tn
xm3 xm2 xm1 xm
xm1 xm2 xm3
x
• We follow this choice with the dilemma of how to reconstruct an
approximation to the derivative of discretely sampled data.
•Assume c is positive
CAAM 452 Spring 2005
Review: Choice of Stencil
tn
tn
tn
tn
tn
xm3 xm2 xm1 xm
xm1 xm2 xm3
x
xm3 xm2 xm1 xm
xm1 xm2 xm3
x
xm3 xm2 xm1 xm
xm1 xm2 xm3
x
xm3 xm2 xm1 xm
xm1 xm2 xm3
x
xm1 xm2 xm3
x
xm3 xm2 xm1 xm
And many more combinations
CAAM 452 Spring 2005
Today: Delta Operators
um1  um
  um 
dx
um  um1
  um 
dx
um1  um1
 0um 
2dx
We will consider the use of each of these difference formula
(and more sophisticated versions) as approximations for the spatial
derivatives.
We consider stability… by which we mean we will examine the
eigenspectrum of the difference operator treated as a matrix.
We will use Gershgorin’s theorem frequently:
http://mathworld.wolfram.com/GershgorinCircleTheorem.html
CAAM 452 Spring 2005
Right Difference Matrix
dum
 c  um
dt
• Consider a case with 10 points on a periodic interval:
• The semi-discrete scheme for the advection equations is:
 u0 
 1 1
  u0 
u 

 u 

1
1
1
 

 1 
 u2 

  u2 
1 1
 

 

1
1
 

 
 
1 1
d   c 
  
 
1 1
dx   dx 
 
 

 
1 1
 

 

1
1
 

 
 

1 1   
 

  
1  u9 
1
 u9 
Note: I used indexing from zero to match with classical definitions coming up.
CAAM 452 Spring 2005
Discrete
  Operator Matrix
• By Gerschgorin’s theorem we conclude that the matrix has all
eigenvalues contained in the circle centered at -c/dx with radius
c/dx
• Thus the eigenvalues are all in the left half plane, with the
possible exception of one or more on the origin  good candidate
for time-stepping.
 u0 
 1 1
  u0 
u 

 u 

1
1
 1

 1 
 u2 

  u2 
1 1
 

 

1
1
 

 
 
1 1
d   c 
  
 
1 1
dx   dx 
 
 

 
1 1
 

 

1
1
 

 
 

1 1   
 

  
u
1

1

  u9 
 9
CAAM 452 Spring 2005
cont
• We can be more precise about the location of the
eigenvalues of  
• We first introduce the Fourier transform for the data
vector:
M 1
uˆn   ume
 im
n 2
M
m 0
• This is simply a linear, map applied to the vector of
m 2
M 1

in
values with inverse given by:
1
M
um 
uˆ e

M
n 0
n
CAAM 452 Spring 2005
Fourier Transform in Matrix Notation
• Introducing the following notation for the Fourier
transform matrix and a shift operator:
uˆn 
m  M 1

m 0
Fnmum where Fnm : e
 im
n 2
M
 un1 if n  M  1
R nmum  
u1 if n  M  1
• We consider the following eigenvalue problem:
k M
c
c un    R nk  I nk  uk  un
dx k 1
or
c
 R  I  u  u
dx
CAAM 452 Spring 2005
Right Shift Matrix
• We write down explicitly what the right shift matrix
is:
R nm   n  1,m
• Where the box around the addition operator
indicates that this is addition modulo M
• This time the delta indicates the Kronecker Delta:
• http://mathworld.wolfram.com/KroneckerDelta.html
CAAM 452 Spring 2005
cont
• We perform a similarity transform using the Fourier
c
matrix:
 R  I  u  u
dx
• Becomes: c
1
F  R  I  F v  v
dx
• We now consider the first term: FRF 1
• First note (basic properties of the Fourier
transform):
n 2
 im
Fnm  e
M
and  F 1 nm
1 imnM2
 e
M
CAAM 452 Spring 2005
cont
1
FRF

nm
1

M
1

M

e
e
M 1 M 1  i 2 nj
M
 e
j 0 k 0
 j  1,k e
M 1  i 2 nj i 2 m j 1
M
M
e
e
i 2 mk
M
What we discover after some algebra
(with slight liberties taken with the
notation) is that the Fourier transform
matrix diagonalizes the right shift matrix.
j 0
i 2 m
M M 1  i 2 nj i 2 mj
M
M
M
i 2 m
M
e
e
j 0
 nm
CAAM 452 Spring 2005
cont
• Finally we return to the   and can now find its
c
eigenspectrum:
1
dx
• Using
F  R  I  F v  v
1
FRF

nm  e
i 2 m
M
 nm
• We find that the   matrix is diagonalized by the
Fourier matrix and that
c
  e
dx 
i 2 m
M

 1 for m  0,..., M  1

CAAM 452 Spring 2005
Important Result: Fourier Transform of
Right Shift Matrix
 FRF 
1
nm
e
i 2 m
M
 nm
i.e. the Fourier (similarity) transform of the
right shift matrix is a diagonal matrix
CAAM 452 Spring 2005
Left Difference Matrix
dum
 c um
dt
• Consider a case with 10 points on a periodic interval:
• The semi-discrete scheme for the advection equations is:
1  u0 
 u0 
1
u 
 1 1
 u 
1
 

 1 
 u2 

  u2 
1 1
 

 

1
1
 

 
 
1 1
d   c 
  
 
1 1
dx   dx 
 
 

 
1 1
 

 

1
1
 

 
 

 
1 1
 

  
1 1   u9 

 u9 
CAAM 452 Spring 2005
cont
• Applying Gerschgorin this time we find that the
eigenvalues are all in a disk centered at +c/dx with
radius c/dx.
• This is bad news !.
• All the eigenvalues are in the right half plane, apart
from the possiblity of a subset which are at the
origin.
• i.e. any non zero mode will be unstable for time
stepping operator (and clearly not physical)
CAAM 452 Spring 2005
Detailed Study of the
Left Difference Operator Spectrum
• This time we introduce the left shift operator:
Lnm   n,m  1
• This time the eigenvalue problem to consider is:
c
c u   I  L  u  u
dx
• Through a similar analysis (since the Fourier transform also
diagonalizes the left shift matrix) we find:
i 2 m



c
M
  1  e
for m  0,..., M  1

dx 

• This confirms the Gerschgorin estimate.
CAAM 452 Spring 2005
Important Result: Fourier Transform of
Left Shift Matrix
 FLF 
1
nm
e
 i 2 m
M
 nm
i.e. the Fourier (similarity) transform of the
left shift matrix is a diagonal matrix
CAAM 452 Spring 2005
Polynomials of Shift Matrices
• We can express most finite difference operators as
polynomials of shift operators: du  r p
j k
  a jk R L  u
dt  j 0 k 0

• We can apply the Fourier similarity transform to diagonalize
p
r
this:

dv
1
1 j
1 k 
 F  a jk  FRF   FLF   Fv
dt
 j 0 k 0

• Where v  F1u is the vector of Fourier coefficents for u.
• Hence each Fourier coefficient satisfies:
p
r

dv
1
1 j
1 k 
 F  a jk  FRF   FLF   Fv
dt
 j 0 k 0

i 2 m j k 
i 2 mj  i 2 mk
p
r



dvm  r p
M
M
M
  a jk e
e
 vm   a jk e
 vm
dt  j 0 k 0

 j 0 k 0

CAAM 452 Spring 2005
Centered Differencing
dum
 c 0um
dt
• We recall the centered difference operator can be
expressed in terms of the left and right difference
operators:
1
0 
     
2dx
c
 c 0u 
 R  L  u  u
2dx
c
1

F  R  L  F v  v
2dx
 im 2
 ic
c  imM2
m 2 

M
 
e
e
 sin 
 m  0,.., M  1


2dx 
 dx  M 
CAAM 452 Spring 2005
cont
• So the  0 operator has purely imaginary eigenvalues
ic
m2 

  sin 
 m  0,.., M  1
dx  M 
• As before we must resort to a time-stepping
method which will in this case accommodate a
portion of the imaginary axis.
CAAM 452 Spring 2005
The 4th Order Central
Difference Formula
 dx 2 2 
dum
  4um :  0 1 
  um
dx
3!


• Recall the fourth order central difference formula:
dum
c

 um2  8um1  8um1  um2 
dt
12dx
du
c

R 2  8R  8L  L2  u

dt
12dx
duˆm
c
2 m

ei 2m  8eim  8e im  e 2im uˆm where  m 
dt
12dx
M
ic

sin  2 m   8sin  m   uˆm

6dx


• As for the  0 central differencing formula all the eigenvalues
are imaginary.
• We have used the diagonalizing property of the Fourier
transform on polynomials of the shift matrices
CAAM 452 Spring 2005
Summary of Mode Evolution
• Using: m 
2 m
M
• We write down the ODE satisfied by each discrete
Fourier mode of the transformed data:
du
dt
du
dt
du
dt
du
dt
duˆm c im
 c  u

  e  1 uˆm
dt
dx
duˆm c
 c u

 1  e  im  uˆm
dt
dx
duˆm ic
 c 0u

 sin  m  uˆm
dt
dx
 dx 2 2 
duˆm
ic
 c 0 1 
 u 

8sin  m   sin  2 m   uˆm
6
dt 6dx


CAAM 452 Spring 2005
Interpreting the ODEs
• We can assert an interpretation of multipliers in the
ODE for each discrete Fourier mode.
• We know that the m’th mode has physical shape e
where L is the periodic length (L=M*dx)
i
2 mx
L
• So (reverting back to our approach for turning a
PDE into an ODE analytically) we would ideally be
solving: duˆm
2 m
dt
 ci
L
uˆm
• Thus we can perform a small theta analysis of the
multipliers in each of the cases we just derived.
CAAM 452 Spring 2005
(1) Right Difference
duˆm c im
  e  1 uˆm
dt dx
c im
c
 m2  m3  m4
4 
e  1   i m   i   O  m  

dx
dx 
2!
3! 4!

c  i 2 m  c   m2  m3  m4
4 
 
     i   O  m  
dx  M  dx  2!
3! 4!

ic 2 m c   m2  m3  m4
4 

    i   O  m  
L
dx  2!
3! 4!

 uˆm  uˆm  0  e
 1 2 m 2 1  2 m 3 

2 m
i
ct ct  dx  M   dx O  M  


L
e
The small theta analysis (low wave number) implies the solution will decay
CAAM 452 Spring 2005
(2) Left Difference
duˆm c
 i m
 1  e  uˆm
dt dx
2
3
4

c
c



 i m
4 
m
m
m
1  e    i m   i   O  m  

dx
dx 
2!
3! 4!

c  i 2 m  c   m2  m3  m4
4 
 
    i   O  m  
dx  M  dx  2!
3! 4!

ic 2 m c   m2  m3  m4
4 

   i   O  m  
L
dx  2!
3! 4!

 uˆm  uˆm  0  e
 1  2 m 2 1
3
2 m
i
ct ct  2!dx  M   dx O  m

L
 
e



The small theta analysis shows that the solution propagates and grows exponentially fast
CAAM 452 Spring 2005
(3) Center Difference
duˆm ic
 sin  m  uˆm
dt dx
ic
ic 
 m3
5 
sin  m     m   O  m  
dx
dx 
3!

c  i 2 m  ic   m3
5 
 
     O  m  
dx  M  dx  3!

ic 2 m ic   m3
5 

    O  m  
L
dx  3!

 uˆm  uˆm  0  e
i

e
 
 i 3 it
2 m
 mt  O  m5
ct c
3!dx
dx
L
The small theta analysis shows that the solution propagates and does not grow
CAAM 452 Spring 2005
(4) 4th Order Center Difference
duˆm
ic

8sin  m   sin  2 m   uˆm
dt 6dx
1 3 1 5
 

7 
8






O



m
m
m
m 

3!
5!
ic  8sin  m   ic

 




6dx   sin  2 m   6dx  
1 3 32 5
7 

2


8




O


  m
m
m
m 
3!
5!

 
ic 2 m ic  24 5
7 





O



m
m 
L
6dx  5!

 uˆm  uˆm  0  e
i

e
 
 i 4 5 it
2 m
 mt  O  m7
ct c
dx 5!
dx
L
CAAM 452 Spring 2005
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
Terms
• The dominant error is dispersive if it causes different modes
to travel at different speeds as with  0
 2 m  1  3 
2 m 1
i
3
m
i

 ct ct O 





3!dx
m 
uˆm  uˆm  0  e L 3!dx e dx
 c  c L

c
1

O
 
 mdx  
2 m



 
3
m
5
m


L
• To leading order accuracy ctilde is the actual numerical
wavespeed, compared with the physical advection speed
• The dominant error is dissipative if it reduces the wave
amplitude over time as with  
uˆm  uˆm  0  e
i

 
1 2 1
2 m
3
ct ct  m  O  m
dx
dx
L
e
CAAM 452 Spring 2005
Comparison of Multipliers
6
CAAM 452 Spring 2005
The Imaginary Part
• We now see that for each of
the methods has a dispersive
nature (wave speed changes
with wave number).
• Also, the imaginary part of the
multiplier is not monotone and
linear over [0, pi] or [pi,2*pi]
• Clearly there will be spurious
modes (eigenvalues) which
show up in the discrete
spectrum as we saw last time.
Spurious modes
CAAM 452 Spring 2005
Recap Via a New Example
• There is a 6th order central difference scheme:
 dx 2
dx 4 2 2 
 6   0 1 
   
 
6
30


• We Fourier transform this:
 dx 2
dx 4 2 2  1
F 6 F  F 0  1 
   
    F
6
30


1
 dx 2 1

dx 4 1
1
 F 0  1 
F F  F F  
F F  F 1F F 1F F 1F   F 1
6
30


 dx 2

dx 4
1
1
 F 0F  1 
F  F F  F  
F  F 1F  F 1F  F 1F  F 1


6
30


2
2
i sin  m 
1
1
  nm
1  1  eim 1  e  im    1  eim  1  e  im 
dx
6
30
1


CAAM 452 Spring 2005
cont
• We simplify this to:

sin  m 
1
1
i m
 i m
i m 2
 i m 2
F 6 F   nmi
1  1  e 1  e    1  e  1  e 
dx
6
30
sin  m   2 2   m  8
4  m 
  nmi
1

sin

sin

 
 
dx  3
2
15
 
 2 
1

• Compare with GKO p83.
• The eigenvalues are all imaginary
• Tedious Taylor expansion computations will reveal
an error O(theta^7)
• Next we plot this dispersion relationship against the
other examples:
CAAM 452 Spring 2005
Phase Error Plot
• It should be clear that the 6th order central difference operator is more
accurate over a wider range of theta – but still prone to spurious
eigenvalues.
Should be
(1/6)*dx^2
CAAM 452 Spring 2005
Summary So Far
• We have established that for the advection equation
u
u
it is not appropriate to use   for the spatial
c
t
x
derivative term (non-negative real part for multiplier).
• The other difference formulae are stable in this
context:   ,  0 ,  4 ,  6
• In this context   is referred to as an upwind derivative
– it uses a dissipative mechanism to remain stable.
• The given central difference operators are carefully
designed so that their eigenvalues sit on the imaginary
axis.
CAAM 452 Spring 2005
cont
• We have established that in each of the schemes we have
replaced the exact ODE multiplier
i 2 m
c
L
• with an approximation, in this lecture denoted by:
i m
c
 error terms
dx
• We found that this approximation is not uniform in wave
number – i.e. there is a non monotonic dependence in m
(embedded in the error terms).
• This explains our observation of the spurious eigenvalues in
the computer evaluations of the spectrum of the discrete
operator.
CAAM 452 Spring 2005
Recall Those Spurious Eigenvalues
• What should the eigenvalues of the 4th order derivative operator
ideally tend to in the limiting of decreasing dx?
• Example (10, 20, 80 points with the periodic length 2*pi):
>> [dxDeigs, dxD] = test4th(10); >> [dxDeigs, dxD] = test4th(20);
>> dx = 2*pi/20;
>> dx = 2*pi/10;
>>
>> sort(dxDeigs/dx)
>> e = sort(dxDeigs/dx);
>> e(1:10)
ans =
0.0000 - 0.0000i
0.0000 + 0.0000i
0.0000 - 0.9950i
0.0000 + 0.9950i
0 - 1.4996i
0 + 1.4996i
0 - 1.8623i
0 + 1.8623i
-0.0000 - 2.1741i
-0.0000 + 2.1741i
ans =
-0.0000
-0.0000
-0.0000 - 0.9997i
-0.0000 + 0.9997i
0.0000 - 1.6233i
0.0000 + 1.6233i
-0.0000 - 1.9901i
-0.0000 + 1.9901i
0.0000 - 2.9290i
0.0000 + 2.9290i
>> [dxDeigs, dxD] = test4th(80);
>> dx = 2*pi/80;
>> e = sort(dxDeigs/dx);
>>
>> e(1:10)
ans =
0.0000 - 0.0000i
0.0000 + 0.0000i
-0.0000 - 1.0000i
-0.0000 + 1.0000i
-0.0000 - 1.6639i
-0.0000 + 1.6639i
-0.0000 - 2.0000i
-0.0000 + 2.0000i
0.0000 - 2.9997i
0.0000 + 2.9997i
CAAM 452 Spring 2005
Upwinding
• Recall: we referred to the   differencing operator
as the “upwind derivative”.
• This is due to the physical placement of the
sampled data relative to the advection direction:
Advection direction
xm3 xm2 xm1 xm xm1 xm2 xm3
• i.e. the differencing only uses data from the
upstream (upwind) direction.
x
• More on this next time.
CAAM 452 Spring 2005
Next Lecture
•
Introduction of some standard finite-difference methods.
•
Define (for finite difference schemes):
1)
2)
3)
4)
•
Consistency
Stability
Accuracy
Courant-Friedrichs-Lewy condition for stability
We will examine convergence:
1) Continuous time+discrete space (as seen today)
2) Discrete time and discrete space
3) Lax-Richtmyer equivalence theory (consistency+stability
 convergence)
CAAM 452 Spring 2005