Numerical Methods for Partial Differential Equations
Download
Report
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ˆ
icuˆ 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 icuˆ
t
x t
x 2
2
dt
duˆ
icuˆ 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ˆ
icuˆ 0
ict
ict ˆ
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 eict 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
icuˆ 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ˆ
icuˆ d 2uˆ
dt
duˆ
ic d 2 uˆ
dt
CAAM 452 Spring 2005
cont
• Again, we can solve this trivial ODE:
duˆ
ic 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ˆ
icuˆ 0
dt
u
u
2u
c d 2
t
x
x
duˆ
ic 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
un1 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
un1 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
un1 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:
un1 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
un1 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 un1 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 tn1 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 tn1 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
un1 un dt un 1 dt un
We can solve this immediately:
u0 u 0
un1 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
icuˆ 0 ic
• 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
ic d uˆ
dt
• This time the stability region must contain:
dt ic 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
ic 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 un1 un dt af un bf un1
• Recall the integral based interpretation of EulerForward:
tn 1
u tn1 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 tn1
t tn
I1 f
f un
f un1
tn tn1
tn1 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 tn1
t tn
I1 fdt
f un
f un1 dt
t tn1
tn1 tn
tn n
tn 1
f un
dt
f un1
t t tn1 dt dt
n
tn 1
tn 1
t t dt
n
tn
f un t
f un1 t
tn1t
t nt
dt 2
dt 2
tn
tn
2
tn 1
2
tn 1
3
1
dt f un f un1
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:
un1 un
tn 1
If dt
1
tn
3
1
un f un f un1
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