Transcript Numerical Methods for Partial Differential Equations

Numerical Methods for Partial
Differential Equations
CAAM 452
Spring 2005
Lecture 3
AB2,AB3, Stability, Accuracy
Instructor: Tim Warburton
Recall: Euler-Forward Time Stepping
• By several different approximations we arrived at the
following approximate ODE time-stepping scheme:
u0  u  0 
du
 f u 
dt
un1  un  dtf  un 
• Assuming: f  u   u for some  
• We established the following stability region 1  dt   1
Im
x
-1
Re
• i.e. if mu*dt is in the yellow region the approximate solution
will be bounded.
CAAM 452 Spring 2005
Recall: Adams Bashford 2
• We went on to derive a similar method which uses
f  u  tn   and f  u  tn1  
3
1

• Namely AB2: un1  un  dt  f  un   f  un1  
2
2

• This was based on an interpolation of f(u(t))
• We next ask which values of mu*dt is this stable for?
CAAM 452 Spring 2005
Derivation of AB3
•
Suppose we wish to use three previous values of f
•
f  u  tn   ,f  u  tn1   and f  u  tn2  
As previously we look for a fit through this data,
although this time we choose a different path.
(i) Consider dt = 1 (will rescale later)
(ii) Will form the coefficients of the interpolant using a
Vandermonde matrix and a monomial basis in t
(iii) Finally integrate the integrand (easy with this
choice of basis).
CAAM 452 Spring 2005
Simpler
Approach
If 2  a  bt  ct 2
f 2  a  b 2  c 4
(i) form quadratic f 1  a  b  c
interpolant: f 0  a
(ii) Compute
coefficients:
 a   1 2 4 
  b    1 1 1 
  

 c  1 0 0 
  

1
 f 2 
f 
 1 
 f 
 0
1 2 4 
 If 2  1 t t 2  1 1 1 


1 0 0 


1
 f 2 
f 
 1 
 f 
 0
1 2 4 
1
1
  If 2 dt   1 t t 2  1 1 1 


0
0
1 0 0 


(iii) Integrate
interpolant:
1 2 4 
 1 1 

 1
 1 1 1 
 2 3 

1
0
0


1
1
 f 2 
 f  dt
 1 
 f 
 0
 f 2 
 f   1  23 f  16 f  5 f 
0
1
2
 1  12
 f 
 0
AB3
• We now have an interpolant for dt=1 based on
three time levels: 1
1
0 If2dt  12  23 f0  16 f1  5 f2 
• We rescale and shift:
tn 1

tn
dt
If 2 dt   23 f n  16 f n1  5 f n2 
12
• And finally have AB3:
dt
un1  un   23 f n  16 f n1  5 f n2 
12
CAAM 452 Spring 2005
Absolute Stability of AB2
• We again assume: f  u   u for some  
3
1

• AB2 becomes: un1  un  dt  un  un1 
2
2

• Tidying up:
3 
1

un1   1    un   un1  0   dt 
2 
2

• To start this method we require:
u1  u  dt 
u0  u  0 
CAAM 452 Spring 2005
Absolute Stability of AB2 cont
• We define the absolute stability region of this
method as the set of nu in the complex plane for
which the sequence {un} is bounded as n  
3 
1

un1   1    un   un1  0
2 
2

• We will make continuous reference to the stability
polynomial of this recurrence relationship (obtained
n
by setting : un  z and canceling excess powers of
z)
3 
1

2
  z   z   1    z  
2 
2

CAAM 452 Spring 2005
Note on Solutions of Recurrence Relations
• Given an s+1 term recurrence relationship (concisely written as):
j s
 u
• With stability polynomial:
j 0
j n j
0
j s
  z    j z j
j 0
• which (we suppose) has s distinct roots
rm m1
m s
m s
• Then any solution to the recurrence relation can be written as: un
  rmn cm
m 1
• where the cm depend on the s initiating values for the recurrence.
n
u

r
• This follows because n
m is a solution, and there is a limited rank of
possible solutions because of the linearity of the recurrence in the initial
values.
n 0
un n1s
Board demo
CAAM 452 Spring 2005
cont
• Since all solutions under these conditions can be
m s
written as:
un   rmncm
m 1
• We immediately observe that this is bounded if and
only if rm  1
• See Trefethen p43-46 for details of the case where
the stability polynomial has 1 or more roots with
multiplicity greater than 1.
CAAM 452 Spring 2005
Root Condition for Absolute Stability
See: Trefethen p58-60
• A linear multistep formula is absolutely stable for a
particular value of nu if and only if all the zeros of
the stability polynomial   z  satisfy z  1 and any
zero with z  1 is a simple root.
CAAM 452 Spring 2005
Drawing the Stability Region
• We need to find the roots of the stability polynomial:
3 
1

  z   z   1    z  
2 
2

2
• i.e. r 
such that   r   0
• Which lie on the margin of stability r  1
i
r

e
for  [0,2 )
• The method we adopt is to set
and find the value of nu such that r is a root of the
stability polynomial.
CAAM 452 Spring 2005
cont
• We can rearrange the root equation to find nu in
terms of z:
3  1

2
  r   r   1    r    0
2  2

r 2  r e 2i  ei
 

3 1 3 i 1
r
e 
2
2 2
2
• We can now ask matlab to plot this curve in the
complex plane…
CAAM 452 Spring 2005
Matlab Plot of Margin of Stability for AB2
CAAM 452 Spring 2005
AB2 v. AB3 v. AB4
• I repeated this for AB2, AB3, AB4
• Notice how the overall region of absolute stability gets
smaller as we use more and more previous values of f
• This translates into stricter restrictions on dt to force
nu=mu*dt into the stability regions.
CAAM 452 Spring 2005
Zooming Into The
All Important Imaginary Axis
• Recall, we discounted Euler-Forward as useless for timestepping the advection equation because it had zero
coverage of the imaginary axis.
• This picture does not quite tell us if
the imaginary axis is inside the absolute
stability region for AB2, AB3 and AB4.
• So we had better check the root condition.
CAAM 452 Spring 2005
Root Condition for AB2
• We need to check that the roots of:
3 
1

2
z   1    z    0
2 
2

• Are bounded above by 1, and if any roots are on
the unit circle we must make sure they are simple.
• We set nu=i*alpha
3 
1

z   1  i  z  i  0
2 
2

2
CAAM 452 Spring 2005
Always Experiment First
• So rather than waste time figuring out where the
roots are, why not let matlab do the hardwork:
CAAM 452 Spring 2005
If at first..
• Since the plot is not quite clear, we plot:
sign(maximum|root|-1)*log10(|maximum|root|-1|)
In summary: the only stable point on the imaginary axis for AB2 is the origin.
AB3 is stable for nu in i*[-.723…,.723…]
CAAM 452 Spring 2005
Accuracy & Consistency
Section 1.3 p21-29 Trefethen
• We already established experimentally that stability
is not enough for a good solution – we also need to
examine how accurate a solution will be.
CAAM 452 Spring 2005
Consistency
• Linear Multistep schemes (like Euler-Forward,
AB2,AB3..) can be written as:
j s
 u
j 0
j n j
j s
 dt   j f  un j 
j 0
• It is natural to ask what is the local discretization
error (or local truncation error) for the above
formula when we substitute u  tn  for u n
• We define the residual operator
j s
j s
j 0
j 0
Lu  tn  :   j u  tn j  dt   j f  u  tn j  
CAAM 452 Spring 2005
Consistency cont
• If the linear multistep method is accurate the
truncation error should be small.
j s
j s
j 0
j 0
Lu  tn  :   j u  tn j  dt   j f  u  tn j  
• Definition: if Lu  tn   O  dt p1  as dt  0 then we say
the scheme is p’th order accurate.
CAAM 452 Spring 2005
Example
• Euler-Forward has: s  1, 0  1,1  1,  0  1,  0  1
• So:
j s
j s
j 0
j 0
Lu  tn  :   j u  tn j  dt   j f  u  tn j  
 u  tn1   u  tn   dtf  u  tn  
 u  tn1   u  tn   dt
du
 tn 
dt
• The next step is to estimate this in terms of dt and
time derivatives of the solution.
CAAM 452 Spring 2005
Recall: Taylor’s Theorem
• Assuming that f is sufficiently smooth (has n bounded
derivatives):
• Then we can evaluate f at x+delta by a linear sum of f at x
and its derivatives:
df
 2 d2 f
 3 d3 f
 n1 d n1 f
f t     f t    t  
t 
t  ... 
t  Rn
2  
3  
n 1  
dt
2! dt
3! dt
 n  1! dt
• The residual takes various forms:
n 1
u t 
  u  d n f
Rn  
u du
n  
 n  1! dt
u t
(Lagrange remainder)
Rn 
 n dn f
n! dt
n
t  for some t  t , t   
*
*
•http://mathworld.wolfram.com/TaylorsTheorem.html
•http://mathworld.wolfram.com/LagrangeRemainder.html
CAAM 452 Spring 2005
Back to Forward-Euler
• We can now express the local truncation error for
Forward-Euler in terms of:
du
d 2u
u  tn  ,  tn  , 2  tn  ,...
dt
dt
• expanding about tn we obtain:
du
Lu  tn   u  tn1   u  tn   dt  tn 
dt

du
dt 2 d 2u * 
du
 u  tn   dt  tn  
t  u  tn   dt  tn 
2  n 
dt
2! dt
dt


dt 2 d 2u *

t
2  n
2! dt
CAAM 452 Spring 2005
cont
dt 2 d 2u *
2
Lu  tn  
t

O
dt


2  n
2! dt
This indicates that the Euler-Forward method
is “first order in time”.
So what does this mean for solution accuracy?
CAAM 452 Spring 2005
Existance and Uniqueness for
the Initial Value Problem
Theorem (p13 Trefethen):
Let f(u,t) be continuous with respect to t and uniformly
Lipschitz continuous with respect to u for t   0, T  .Then there
exists a unique differentiable function u(t) that satisfies the
initial-value problem:
u  0   u0
du
 f  u  t  , t  for all t   0, T 
dt
The important condition here is that the right hand side f is not too non-linear.
With these conditions – we know we have a single solution to aim for.
CAAM 452 Spring 2005
Definition: Convergence
• A linear multistep method is convergent if, for all
IVP’s satisfying the existance/uniqueness theorem
conditions on an interval [0,T] and all starting
values v0 , v1 ,..., v1s  for the method tend to u(0) as
dt  0
• the solution vn satisfies
v  t   u  t   o 1 as dt  0
• uniformly for all t  0, T 
The notation small ‘o’ of 1 as dt tends to zero, requires that the solution drops to zero as dt tends to 0.
CAAM 452 Spring 2005
Dahlquist Equivalence Theorem
p54 Trefethen
Theorem:
A linear multistep formula is convergent if and only
if it is consistent and stable.
Proving convergence of a method
directly is difficult, however proving
consistency and stability is
relatively straight forward.
CAAM 452 Spring 2005
Finally: an estimate for solution accuracy
• Theorem (global p’th order accuracy):
• Consider an IVP that satisfies the existance/uniqueness
theorem conditions on an interval [0,T] and in addition that
the right hand side function f(u,t) is p times continuously
differentiable with respect to u and t. Let an approximate
solution be computed by a convergent linear multistep
formula of accuracy >= p with starting values as in the
convergence statement. Then this solution satisfies:
v  t   u  t   O  dt p  as dt  0
• uniformly for all t   0, T 
CAAM 452 Spring 2005
Summary
• Stability + consistency  convergence
• convergence + method accuracy  global solution
accuracy
CAAM 452 Spring 2005
Sketch of Solution Accuracy Proof
(not entirely rigorous)
• Given the way we constructed the AB schemes it is
instructive to write down the following three
equations:
T
 i  u T   u  0    f  u  t   dt
0
T
n N
 ii  u T   u  0    I p f  u  t   dt   L  u  tn  
0
N=T/dt
n 1
T
 iii  u T   u  0    I p f  u  t   dt
0
CAAM 452 Spring 2005
Sketch cont
T
 i  u T   u  0    f  u  t   dt
0
T
n N
0
n 1
 ii  u T   u  0    I p f  u  t   dt   L  u  tn  
T
 iii  u T   u  0    I p f  u  t   dt
0
(i) Represents the actual solution (which we know exists and is unique
by the Lipschitz assumptions on f)
(ii) Represents acknowledges the impact of using an interpolant for f on the
ODE (i.e. the local truncation errors represent the residuals)
(iii) Represents the numerical solution u tilde.
CAAM 452 Spring 2005
Sketch cont
• Part 1: we used (i)-(ii) (with T=dt) to estimate L for
each time step.
 ii  -  i  L  u  tn    O  dt p1 
• Part 2: we look at (ii)-(iii)
T
n N
0
n 1
 ii  u T   u  0    I p f  u  t   dt   L  u  tn  
T
 iii  u T   u  0    I p f  u  t   dt
0
T
n N
0
n 1
 ii  -(iii) u T   u T   u  0   u  0    I p  f  u  t    f  u  t    dt   L  u  tn  
CAAM 452 Spring 2005
Sketch cont
• We can now try to bound each term:
T
u T   u T   u  0   u  0    I p  f  u  t    f  u  t    dt 
n N
 L u t  
0
T
 u  0   u  0    I p  f  u  t    f  u  t    dt 
0
Bounded by assumption on
Initial conditions
n 1
n
n N
 L u t  
n 1
n
(i) Bounded by assumption on
consistency.
(ii) Estimated by accuracy study.
Bounded by assumption on
Stability of the Ip interpolation
operator
CAAM 452 Spring 2005
Euler-Forward Solution Accuracy
Consider:
du
u  0   u0 ,   u
dt
1) the right hand side function is certainly Lipschitz!
2) It is certainly p=1 times differentiable
3) If we keep mu*dt in the absolute stability region (unit circle
centered at -1 in the complex plane) then it is stable.
4) By Taylor’s theorem it is a p=1 order accurate method.
5) All of the above indicate global solution accuracy which is
first order in dt, i.e. the error will be O(dt) as dt0
CAAM 452 Spring 2005
AB2 Global Accuracy
• We already established the stability region for AB2
– so as long as we keep mu*dt in that region we
have stability.
• Now we have to consider accuracy:
1 du
 3 du
Lu  tn   u  tn1   u  tn   dt 
 tn  
 tn1  
2 dt
 2 dt


du dt 2 d 2u dt 3 d 3u* 
 u  dt 

u
2
3 
dt 2! dt
3! dt 

 3 du 1  du
d 2u dt 2 d 3u**  
 dt 
   dt 2 
3 
2
dt
2
dt
dt
2!
dt



dt 3  d 3u *  dt 3  d 3u ** 

t 
t


3  n 
3  n 
3!  dt
 2!2  dt

 O  dt 3 
AB2 is a 2nd order
acurate method
CAAM 452 Spring 2005
cont
• If the AB2 stability condition is met then the global
accuracy result indicates that for a solution at fixed
time the solution will be
O  dt 2 
• Thus we characterize the scheme as 2nd order
method accurate, and 2nd order solution accurate.
CAAM 452 Spring 2005
AB3 Convergence
Expand all terms in the local truncation error about tn
Lu  tn   u  tn 1   u  tn   dt
1  du
du
du

23
t

16
t

5
t







n
n 1
n2 
12  dt
dt
dt


du dt 2 d 2u dt 3 d 3u dt 4 d 4u* 
 u  dt 


u
2
3
4 
dt 2! dt
3! dt
4! dt 

 du
 du
d 2u dt 2 d 3u dt 3 d 4u**  

 23  16   dt 2 

3
dt
dt
dt
2!
dt
3! dt 4  
1

dt
12   du
d 2u 4dt 2 d 3u 8dt 3 d 4u***  

 5   2dt 2 

3
4 
dt
dt
2!
dt
3!
dt
 
 
dt 4  d 4u *  2dt 4  d 4u **  5dt 4  d 4u *** 

t 
t

t



4  n 
4  n 
4  n 
4!  dt
9  dt
9  dt



 O  dt 4 
 The AB3 scheme is 3rd order accurate as long as the
solution is sufficiently smooth.
CAAM 452 Spring 2005
Homework 1
Read chapter 1 of Trefethen’s book.
Q1) derive AB4 (use Matlab or Mathematica)
Q2) plot the stability region for AB4
Q3) indicate which part of the complex plane is inside
the AB4 region of absolute stability (you may need
to test at least one point in each closed loop of the
marginal stability curve)
continued on next slide
CAAM 452 Spring 2005
Homework 1 cont
Q4a) estimate upper bounds for dt to solve the following with AB2,AB3, and
AB4
u  0  1
du
 u
dt
  .1  .5i
t  0,5
Q4b) write a code which uses Euler-Forward for the first time step, AB2 for
the second time step and AB3 for the third time step (use a suitable dt
for AB3 for all time steps)
Q4c) use this code and plot the error as a function of time for dt, dt/2, dt/4,
dt/8 (i.e. run the simulation 4 times).
Q4d) estimate the order of accuracy of the code based on the error at t=5
Q4e) (extra credit) compare the code with an AB3 method started with
exact initial data.
CAAM 452 Spring 2005