Transcript Numerical Methods for Partial Differential Equations

Numerical Methods for Partial
Differential Equations
CAAM 452
Spring 2005
Lecture 2
Instructor: Tim Warburton
Note on textbook for finite difference methods
• Due to the difficulty some students have experienced in
obtaining Gustafasson-Kreiss-Oliger I will try to find
appropriate and equivalent sections in the online finitedifference book by Trefethen:
http://web.comlab.ox.ac.uk/oucl/work/nick.trefethen/pdetext.html
CAAM 452 Spring 2005
Recall Last Lecture
• We considered the model advection PDE:
u
u
c  0
t
x
• defined on the periodic interval [0,2pi)
• We recalled that any 2pi periodic, C1 function could be represented as a
uniformly convergent Fourier series, so we considered the evolution of
the PDE with a single Fourier mode as initial condition. This converted
the above PDE into a simpler ODE for the time-dependent coefficient
(i.e. amplitude) of a Fourier mode:
duˆ
 icuˆ  0
dt
CAAM 452 Spring 2005
Recall: the Advection Equation
• We wills start with a specific Fourier mode as the initial
condition:
1) Find 2 -periodic u  x, t  such that x   0, 2  , t  0, T 
u
u
c  0
t
x
given
1 i x ˆ
u ( x,0)  f  x  =
e f   x   0, 2 
2
where f is a smooth 2 -periodic function of one frequency 
• We try to find a solution of the same type:
1 i x
u  x, t  
e uˆ  , t 
2
CAAM 452 Spring 2005
cont
• Substituting in this type of solution the PDE:
u
u
c  0
t
x
• Becomes an ODE:
u
u  
   1 i x
1 i x  duˆ

 c    c 
e uˆ  , t   
e   icuˆ 
t
x  t
x   2
2
 dt


duˆ

 icuˆ  0
dt
• With initial condition
uˆ ,0   fˆ  
CAAM 452 Spring 2005
cont
• We have Fourier transformed the PDE into an ODE.
• We can solve the ODE:
duˆ
 icuˆ  0 
ict
ict ˆ

u

,
t

e
u

,0

e
f  
ˆ
ˆ
dt
 
 

uˆ  ,0   fˆ   
• And it follows that the PDE solution is:



1 i  x ct  ˆ

u
x
,
t

e
f    f  x  ct 



2

1 i x ˆ
initial condition: f  x  
e f   
2

1 i x
e uˆ  , t 
2
solution : uˆ  , t   eict fˆ  
ansatz : u  x, t  
CAAM 452 Spring 2005
Note on Fourier Modes
• Note that since the function should be 2pi periodic
we are able to deduce:  
• We can also use the superposition principle for the
more general case when the initial condition
contains multiple Fourier modes:
1   i x ˆ
f  x 
e f  

2  
1   i  x ct  ˆ
 u  x, t  
e
f    f  x  t 

2  
CAAM 452 Spring 2005
cont
• Let’s back up a minute – the crucial part was when
we reduced the PDE to an ODE:
u
u
duˆ
c  0
 icuˆ  0
t
x
dt
• The advantage is: we know how to solve ODE’s
both analytically and numerically (more about this
later on).
CAAM 452 Spring 2005
Add Diffusion Back In
• So we have a good handle on the advection
equation, let’s reintroduce the diffusion term:
u
u
 2u
c  d 2
t
x
x
• Again, let’s assume 2-pi periodicity and assume the
same ansatz:
• This time:
u
u
u
c  d 2
t
x
x
2
1 i x
u
e uˆ  , t 
2
duˆ
 icuˆ   d 2uˆ
dt
duˆ
  ic  d 2  uˆ
dt
CAAM 452 Spring 2005
cont
• Again, we can solve this trivial ODE:
duˆ
  ic  d 2  uˆ
dt
uˆ  uˆ  ,0  e
u  uˆ  ,0  e
ic d 2 t
i  ct  x   d 2t
e
CAAM 452 Spring 2005
cont
• The solution tells a story:
u  uˆ  ,0  e
i  ct  x   d 2t
u  f  ct  x  e
e
 d  2t
• The original profile travels in the direction of
decreasing x (first term)
• As the profile travels it decreases in amplitude
(second exponential term)
CAAM 452 Spring 2005
What Did Diffusion Do??
• Advection:
• Diffusion:
u
u
c  0
t
x
duˆ
 icuˆ  0
dt
u
u
 2u
c  d 2
t
x
x
duˆ
  ic  d 2  uˆ
dt
• Adding the diffusion term shifted the multiplier on the
right hand side of the Fourier transformed PDE (i.e.
the ODE) into the left half plane.
• We summarize the role of the multiplier…
CAAM 452 Spring 2005
Categorizing a Linear ODE
Im   
Increasingly
oscillatory 
duˆ
 uˆ  uˆ  uˆ  0  e t
dt
Exponential growth 
Increasingly
 oscillatory
Exponential decay
Re   
Here we plot the dependence of the solution to the top left ODE
on mu’s position in the complex plane
CAAM 452 Spring 2005
Categorizing a Linear ODE
Im   
Re   
Here we plot the behavior of the solution for 5 different choices of mu
CAAM 452 Spring 2005
Summary
duˆ
 uˆ  uˆ  uˆ  0  e t  uˆ  0  e Re  t  cos  Im    t   i sin  Im    t  
dt
• When the real part of mu is negative the solution decays
exponentially fast in time (rate determined by the magnitude
of the real part of mu).
• When the real part of mu is positive the solution grows
exponentially fast in time (rate determined by the magnitude
of the real part of mu).
• If the imaginary part of mu is non-zero the solution oscillates
in time.
• The larger the imaginary part of mu is, the faster the solution
oscillates in time.
CAAM 452 Spring 2005
Solving the Scalar ODE Numerically
• We know the solution to the scalar ODE
duˆ
  uˆ
dt
• However, it is also reasonable to ask if we can solve it
approximately.
• We have now simplified as far as possible.
• Once we can solve this model problem numerically, we will
apply this technique using the method of lines to
approximate the solution of the PDE.
CAAM 452 Spring 2005
Stage 1: Discretizing the Time Axis
du
 f u 
dt
• It is natural to divide the time interval [0,T] into
shorted subintervals, with width often referred to as:
dt , t or k
• We start with the initial value of the solution u(0)
(and possibly u(-dt),u(-2dt),..) and build a
recurrence relation which approximates u(dt) in
terms of u(0) and early values of u.
CAAM 452 Spring 2005
Example
• Using the following approximation to the time
derivative:
du u  t  dt   u  t 
dt

dt
• We write down an approximation to the ODE:
du
 f u 
dt
u  dt   u  0 
 f u  0
dt
CAAM 452 Spring 2005
Example cont
• Rearranging:
u  dt   u  0 
 f u  0
dt
u  dt   u  0   dtf  u  0  
• We introduce notation for the approximate solution
after the n’th time step: u n
• Our intention is to compute un  u  ndt 
• We convert the above equation into a scheme to
compute an approximate solution:
u0  u  0 
u1  u0  dtf  u0 
CAAM 452 Spring 2005
Example cont
• We can repeat the following step from t=dt to t=2dt and so
on:
u0  u  0 
u1  u0  dtf  u0 
u2  u1  dtf  u1 
un1  un  dtf  un 
• This is commonly known as the:
– Euler-forward time-stepping method
– or Euler’s method
CAAM 452 Spring 2005
Euler-Forward Time Stepping
• It is natural to ask the following questions about this time-stepping
method:
u0  u  0 
un1  un  dtf  un 
• Does the answer get better if dt is reduced? (i.e. we take more timesteps between t=0 and t=T)
• Does the numerical solution behave in the same way as the exact
solution for general f?
(for the case of f(u) = mu*u does the numerical solution decay and/or
oscillate as the exact solution)
• How close to the exact solution is the numerical solution?
• As we decrease dt does the end iterate converge to the solution at t=T
CAAM 452 Spring 2005
Experiments
• Before we try this
analytically, we can
code it up and see what
happens.
• This is some matlab
code for Euler forward
applied to:
u0  u  0 
un1  un  dt un
http://www.caam.rice.edu/~caam452/CodeSnippets/EulerForwardODE.m
CAAM 452 Spring 2005
What dt can we use?
• With dt=0.1, T=3, uo=1, mu=-1-2*i
Matlab
The error is quite small.
CAAM 452 Spring 2005
Larger dt=0.5
• With dt=0.5, T=3, uo=1, mu=-1-2*i
The error is visible but the trend is not wildly wrong.
CAAM 452 Spring 2005
Even Larger dt=1
• With dt=1, T=3, uo=1, mu=-1-2*i
The solution is not remotely correct – but is at least bounded.
CAAM 452 Spring 2005
dt=2
• With T=30, dt=2, uo=1, mu=-1-i
Boom – the approximate solution grows exponentially fast, while the true solution decays!
CAAM 452 Spring 2005
Observations
• When dt is small enough we are able to nicely
approximate the solution with this simple scheme.
• As dt grew the solution became less accurate
• When dt=1 we saw that the approximate solution
did not resemble the true solution, but was at least
bounded.
• When dt=2 we saw the approximate solution grew
exponentially fast in time.
CAAM 452 Spring 2005
cont
• Our observations indicate two qualities of time
stepping we should be interested in:
– stability: i.e. is the solution bounded in a similar way
to the actual solution?
– accuracy: can we choose dt small enough for the
error to be below some threshold?
CAAM 452 Spring 2005
u0  u  0 
Geometric Interpretation:
un1  un  dtf  un 
• We can interpret Euler-forward as a shooting method.
• We suppose that un is the actual solution, compute the
actual slopef  un  and shift the approximate solution by dtf  un 
u n 1
un
dtf  un 
dt
• Note:
t
n
tn 1
– the blue line is not the exact solution, but rather the same
ODE started from the last approximate value of u computed.
– in this case we badly estimated the behavior of even the
approximately started solution in the interval dt
 two kinds of error can accumulate!.
CAAM 452 Spring 2005
Interpolation Interpretation:
u0  u  0 
un1  un  dtf  un 
• We start with the ODE: du  f  u 
dt
• Integrate both sides from over a dt interval:
tn 1
tn 1
du
t dt dt  t f  u  dt
n
n
• Use the fundamental theorem of calculus:
tn 1

tn
du
tn 1
dt  u t  un1  un
n
dt
• Finally, replace f with a constant which interpolates
f at the beginning of the interval…
CAAM 452 Spring 2005
cont
tn 1

tn
tn 1
du
tn 1
dt  u t  u  tn1   u  tn 
n
dt
tn 1
tn 1
 f  u  t   dt   f  u  t   dt  f  u t    dt  f u t   dt
n
tn
tn
n
n
tn
 u  tn1   u  tn   f  u  tn   dt
Again, we have recovered Euler’s method.
CAAM 452 Spring 2005
Stability
Let’s choose f(u) = mu*u
u0  u  0 
un1  un  dt un  1  dt   un
We can solve this immediately:
u0  u  0 
un1  1  dt   un  1  dt  
n 1
u 0
Next suppose Re(mu)<=0 then we expect the actual
solution to be bounded in time by u(0).
For this to be true of the approximate solution we
require: 1  dt   1
CAAM 452 Spring 2005
Stability Condition for Euler Forward
• Since mu is fixed we are left with a condition which
must be met by dt
1  dt   1
• which is true if and only if:
2
2
1  Re   dt     Im   dt    1
• The region of the complex plane which satisfies this
condition is the interior of the unit circle centered at
Im
-1+0*i
x
Re
CAAM 452 Spring 2005
cont
• i.e.  dt will not blow up if it is located inside the
unit disk:
Im
x
-1+0i
Re
• Notice: the exact solutions corresponding to
Re(mu)<0 all decay. However, only the numerical
solutions corresponding to the interior of the yellow
circle decay.
CAAM 452 Spring 2005
Disaster for Advection!!!
• Now we should be worried. The stability region:
Im
x
Re
-1+0i
• only includes one point on the imaginary axis (the origin)
but our advection equation for the periodic interval has mu
which are purely imaginary!!!. duˆ
dt
 icuˆ  0    ic
• Conclusion – the Euler-Forward scheme applied to the
Fourier transform of the advection equation will generate
exponentially growing solutions.
CAAM 452 Spring 2005
Nearly a Disaster for
Advection-Diffusion!!!
duˆ
2
  ic  d  uˆ
dt
• This time the stability region must contain:
 dt   ic   2d  dt
Im
x
-1+0i
Re
• Because the eigenvalues are shifted into the left half
part of the complex plane, we will be able to choose
dt small enough to force the mu*dt into the stability
region.
• We can estimate how small dt has to be in this case
CAAM 452 Spring 2005
cont (estimate of dt for advection-diffusion)
• For stability we require:
 dt  1  1
  ic   d  dt  1  1
2
 1  dt d   cdt   1
2
2
2
 2dt d   dt d   cdt   0
2
2
2
2
 dt  0 or
  0 or
dt  d  c
2
2
2

2d
CAAM 452 Spring 2005
cont
• The only interesting condition for stability is:
2d
dt  2 2 2
 d c
• What does this mean?, volunteer?
• i.e. how do  , d , c influence the maximum stable dt?
• Hint: consider
i) d small, c large
ii) d large, c small
CAAM 452 Spring 2005
Multistep Scheme (AB2)
• Given the failure of Euler-Forward for our goal
equation, we will consider using the solution from 2
previous time steps.
• i.e. we consider un1  un  dt  af  un   bf  un1  
• Recall the integral based interpretation of EulerForward:
tn 1
u  tn1   u  tn    f  u  t   dt  f  u  tn   dt
tn
CAAM 452 Spring 2005
cont
• We interpolate through f(un-1)and f(un) then
integrate between tn and tn+1
Linear interpolant
of the f values
f
Ifn+1
f(un)
f(un-1)
tn-1
tn
tn+1
Approximate integral
CAAM 452 Spring 2005
cont
• This time we choose to integrate under the
interpolant of f which agrees with f at tn and tn-1
• The unique linear interpolant in this case is:
 t  tn1 
 t  tn 
I1 f  
f  un   
f  un1 


 tn  tn1 
 tn1  tn 
• We need to compute the following integral of the
interpolant..:
tn
I
1
fdt
tn
CAAM 452 Spring 2005
cont
tn 1

tn
 t  tn1 
 t  tn 
I1 fdt   
f  un   
f  un1  dt


t  tn1 
 tn1  tn 
tn  n
tn 1
f  un 

dt
f  un1 
t  t  tn1  dt  dt
n
tn 1
tn 1
  t  t  dt
n
tn
f  un   t
f  un1   t



 tn1t  
 t nt 


dt  2
dt  2
 tn
 tn
2
tn 1
2
tn 1
3
1

 dt  f  un   f  un1  
2
2

Last step: homework exercise
CAAM 452 Spring 2005
Adams-Bashford 2
• The resulting AB2 scheme requires the solution at
two levels to compute the update:
un1  un 
tn 1
 If dt
1
tn
3
1

 un   f  un   f  un1  
2
2

CAAM 452 Spring 2005
Next Time
• Stability region for AB2
• Generalization to AB3, AB4 (use more historical data to compute
update)
• Accuracy/consistency of Euler-Forward, AB2,AB3,AB4
• Convergence.
• Runge-Kutta schemes (if time permits)
• READING:
http://web.comlab.ox.ac.uk/oucl/work/nick.trefethen/1all.pdf
• p10-p55 (do not do exercises)
CAAM 452 Spring 2005