Transcript Numerical Methods for Partial Differential Equations
Numerical Methods for Partial Differential Equations
CAAM 452 Spring 2005 Lecture 5 Summary of convergence checks for LMM and one-step time stepping Instructor: Tim Warburton
Today
• I will summarize the steps we took to examine the stability, consistency and thus convergence of the AB and RK schemes
CAAM 452 Spring 2005
AB
• The general AB scheme reads:
u n
1
u n
dt
m
1
m n m
• To check linear stability we assume that
f
is linear in
u
(we can also achieve a similar goal by assuming
f
is Lipschitz) • Assuming:
n
u n
u n
1
u n
dt m
1
m u n m
• We know that all the solutions of this recurrence relationship are bounded provided the corresponding stability polynomial:
z s
z s
1
dt
m z m
1 has roots which satisfy the root condition (all roots bounded by 1, and if they have magnitude 1 must be simple roots).
CAAM 452 Spring 2005
AB
cont
• We must then check the roots of:
z s
z s
1
dt
m
1
m z
• For roots with |z|<1 we are sure that the method is stable. What we wish to know is what values of mu*dt guarantee this condition. We can think of this formula defining a map between the (complex) plane of roots to the (complex) plane of mu*dt • We must determine the image of the unit circle in the root plane in the mu*dt plane.
CAAM 452 Spring 2005
AB
cont
• Given the root equation:
z s
z s
1
dt
m
1
m z
• We determine the map in a trivial manner:
z s
z s
1
dt
m
1
m z
dt
z s
z s
1
m
1
m z
• Example AB3: : (23
z
2
z
3 -16
z z
2 5) /12
CAAM 452 Spring 2005
AB3 Stability
z
re i
:
z
3 (23
z
2 -16
z z
2 5) /12
CAAM 452 Spring 2005
AB Accuracy
• Next we perform a local truncation error analysis.
• That is we estimate what correction at each step of the solution is required for the exact solution of the ODE to be a solution to the numerical method:
n
1
n
dt
m
1
m
n m
T n
CAAM 452 Spring 2005
cont
• We now use Taylor expansions to express the truncation error in terms of u and its derivatives at tn:
n
1
n
dt m
1
m
n m
T n
• We found that Tn scaled as:
T n
Cdt s
1
d dt s s
1
u
1
for some constant and
t n
*
n n
1 • Finally we apply Dahlquist’s equivalence theorem, which requires consistency and stability convergence and corollary which gives us global error estimates:
u n
C n
, ,...
dt s
CAAM 452 Spring 2005
RK Analysis
• A similar analysis is required for the Runge-Kutta schemes. • The stability analysis was straightforward as we found that the recurrence for a linear
f
was:
u n
1 • For stability we required:
m
0
m m
!
u n m
0
m m
!
• We examined this by setting
m
0
m m
!
1
e i
CAAM 452 Spring 2005
RK3 Stability
and let Matlab do the work.
re i
1 2 2 3!
3
re i
CAAM 452 Spring 2005
RK Accuracy
• The error estimate process for RK schemes is a little more complicated, but still boils down to dt expansions of the one-step operator.
• Again, once we have consistency and stability we have convergence.
CAAM 452 Spring 2005
Implicit Schemes
• Later on in the semester we will discuss implicit schemes. • Then we will also discuss the Butcher block formulation of RK methods.
CAAM 452 Spring 2005
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
t
c
u
x
• And not just for one mode.
CAAM 452 Spring 2005
Sample Points in Time Space
u
t
c
u
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.
t n x x m
3
x m
2
x m
1
x m x m
1
x m
2
x m
3 • We follow this choice with the dilemma of how to reconstruct an approximation to the derivative of discretely sampled data.
CAAM 452 Spring 2005
t n t n t n t n t n
Choice of Stencil
x x m
3
x m
2
x m
1
x m x m
1
x m
2
x m
3
x m
3
x m
2
x m
1
x m x m
1
x m
2
x m
3
x m
3
x m
2
x m
1
x m x m
1
x m
2
x m
3
x m
3
x m
2
x m
1
x m x m
1
x m
2
x m
3
x m
3
x m
2
x m
1
x m x m
1
x m
2 And many more combinations
x m
3
x x x x
CAAM 452 Spring 2005
Use Taylor’s Series
• In much the same spirit as time-domain we can expand the solution at each grid point in terms of u and its derivatives as a power series in
dx (the space between spatial grid points)
• Examples:
m
1
m
dx du dx
m
2 2
dx d u
m
2!
dx
2 3 3
dx d u
3!
dx
3
m
1
u x m
dx du dx
m
2 2
dx d u
m
2!
dx
2 3 3
dx d u
3!
dx
3
CAAM 452 Spring 2005
Left Differencing
t n x m
3
x m
2
x m
1
x m x x m
1
x m
2
x m
3 • Using the data points on this stencil we wish to compute an approximation to the gradient of u at
x m
m
1
u x m
dx du dx
m
2!
2
dx
2
m
3 3!
dx
3 • Dividing and rearranging we find:
du dx
dx m
1
2
dx d u
2!
dx
2 2 3
dx d u
3!
dx
3
CAAM 452 Spring 2005
Left Differencing
t n x m
3
x m
2
x m
1
x m x m
1
x m
2
x m
3
m u x m
1
dx
2
dx d u
2!
dx
2 2 3
dx d u
3!
dx
3
x du dx
This gives us an obvious one-sided linear interpolation formula of the derivative. We will later on decide whether it is suitable for time-stepping with AB or RK.
CAAM 452 Spring 2005
Left Difference Operator
• We define the following difference operator (appalling but fairly standard notational choice):
u m
u m
u m
1
dx u
m’th data point.
• The difference formula shows
u m
du
m
dx
2
dx d u
m
2!
dx
2 2 3
dx d u
3!
dx
3 i.e. delta- is a first order approximation of the derivative operator
CAAM 452 Spring 2005
Right Difference Operator
t n x m
3
x m
2
x m
1
x m x m
1
x m
2
x m
3
u m
u m
1
u m dx
u m
du
m
dx
2
dx d u
m
2!
dx
2 2 3
dx d u
3!
dx
3
x
i.e. delta+ is a first order approximation of the derivative operator
CAAM 452 Spring 2005
Central Differencing Operator
0
t n x m
3
x m
2
x m
1
x m x x m
1
x m
2
x m
3 • Subtracting the forward and backward expansions:
m
1
m
1
u x m u x m
dx du dx
dx du dx
m
2!
2
dx
2 2 2
dx d u
2!
dx
2
m
3 3!
dx
3 3 3
dx d u
3!
dx
3 4 4!
dx
4 4 4
dx d u
4!
dx
4 • We obtain:
m
1
m
1
2
dx du dx
m
2 3 3
dx d u
3!
dx
3 5 5!
dx
5 5 5
dx d u
5!
dx
5
CAAM 452 Spring 2005
cont
0
u m
:
u m
1 2
dx u m
1
du dx
m
2 3
dx d u dx
3 accurate in dx • Note: for the error term to be bounded we require that the third derivative be bounded.
CAAM 452 Spring 2005
Approximating the Second Derivative
• We can create an approximation of the second 2
u m
derivative: :
u m
1
u m
u m
u m
u m
u m
1
dx
u m
1
u m
1
dx
2 2
u m
u m
1
dx
2
dx
2 1
dx
2
m
m
dx du dx dx du dx
m
m
2 2
dx d u
2!
2
m
dx
2
m
3 3
dx d u
3!
dx
3 2 2
dx d u
2!
dx
2
m
3 3
dx d u
3!
dx
3
m
m
4 4
dx d u
4!
4 4
dx d u
4!
dx
4
dx
4
CAAM 452 Spring 2005
cont
2
u m
1
dx
2 2
d u
m
dx
2
dx du dx
2 2
dx d u
2!
2
dx
2 3 3
dx d u
3!
dx
3 4 4
dx d u
4!
dx
4
m
dx du
m
dx
2 2
dx d u
m
2!
dx
2 3 3
dx d u
m
3!
dx
3 4 4
dx d u
4!
dx
4
CAAM 452 Spring 2005
Delta Operators
u m
du dx
m
2
dx d u
m
2!
dx
2 2 3
dx d u
3!
dx
3
u m
du dx
m
2
dx d u
m
2!
dx
2 2 3
dx d u
3!
dx
3
2
u m
2
d u
m
dx
2
0
u m
du dx
m
2 3
dx d u
3!
dx
3
x m
CAAM 452 Spring 2005
We use:
Using the delta’s: Example 1
0
u m
du
dx
m
2 3
dx d u
2
dx d
2 3!
dx
2 3!
dx du dx
3
m
0
u m
du dx
m
dx
2 3!
2 3
dx d u
3!
dx
3
2
du dx
m
CAAM 452 Spring 2005
Example 2
• We can use the difference operators to eliminate the third order spatial derivative:
0
u m
2
dx d
2 3!
dx
2
du dx
m
du
m
dx dx
2
2 3!
du
m
dx dx
2 3!
0
u m
0
u m
du dx
0
dx
2 3!
0
u m
CAAM 452 Spring 2005
Example 2
cont
• We can expand out the difference operator
du
m dx
0
dx
2 3!
u m
1
dx
2 3!
2
u m
1 2
dx
0
u m
u m
1
u m
1 2
dx
1 3!2
dx
u m
2
u m
2
u m
1 2
u m
1
u m
u m
2 1 12
dx
u m
2 8
u m
1 8
u m
1
u m
2 • And we have a super difference formula for the derivative which is 4 th order accurate!
• What is the stencil for this formula ? (volunteer)
CAAM 452 Spring 2005
Example 2
cont
• We can think of this difference formula as a matrix operator acting on the point data values aligned in a vector:
du dx
m
1 12
dx
u m
2 8
u m
1 8
u m
1
u m
2 • Becomes
CAAM 452 Spring 2005
For 10 Points on a Periodic Interval
d dx
u u u
1 2 3
u
10 1 12
dx
8 0 1 8 1 8 0 8 1 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 1 8 0 8 1 8 0 8 1
u u u
1 2 3
u
10 Notice: the matrix is skew-symmetrix imaginary eigenvalues
CAAM 452 Spring 2005
Diagonalize
• We now use this formula in the PDE which we evaluate at the 10 spatial sample points:
d dt
u u u u u u u u u
1 2 3 4 5 6 7 8 9
u
10
d c dx
u u u u u u u u u
1 2 3 4 5 6 7 8 9
u
10
c
12
dx
0 8 1 8 1 8 0 8 1 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 1 8 0 8 1 1 8 0 8
u u u u u u u u u
1 2 3 4 5 6 7 8 9
u
10 This is an identity (although we are a little fuzzy about the O(dx^4) term)
CAAM 452 Spring 2005
Some Notation
• We will represent the set of data values at the spatial points as a vector:
u
u u u u u u u u u
1 2 3 4 5 6 7 8 9
u
10
CAAM 452 Spring 2005
cont
• The system matrix is an approximation of the derivative operator so we will denote it as
D D
1 12
dx
8 0 1 8 1 8 0 8 1 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 1 8 0 8 8 1 1 8 0
CAAM 452 Spring 2005
PDE
System of ODE’s > Decoupled ODE’s
• So we can express the PDE as:
d dt
u
c
Du
• Let’s make a semi-discrete scheme:
d dt
u
c
Du
• As a theoretical exercise we can diagonalize the skew symmetric
D D
Λ
diagonal matrix. • Then:
d dt
u
d dt
v
1
c
S S u v
where
v
1
CAAM 452 Spring 2005
cont
• The system of 10 ODE’s:
d dt
v
v
• are now totally decoupled, so we can apply our extensive knowledge of treating scalar ODE’s to each of the components of vhat.
• All we have to do is make sure that the eigenvalues of the c
D
matrix are squeezed into the time stepping region by choosing small enough dt
CAAM 452 Spring 2005
cont
• By linearity of the right hand side we can be sure that time stepping (say using AB or RK) the original ODE will be fine region…
cdt
max 1
i
10
i
• We now need to solve:
d dt
u
c
Du
• Noting our studies of the stability regions, the fact that all the eigenvalues sit on the imaginary axis
rules out
– Euler-Forward (unconditionally unstable) – AB2 and RK2 (marginally stable)
CAAM 452 Spring 2005
cont
• In this case we can use Gerschgorin’s theorem to bound the magnitude of the eigenvalues.
D
1 12
dx
0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 0 8 1 1 8 8 1 1 8 0 • The Gerschgorin disk for each row is centered at zero with radius:
12
c dx
8 1 8
3
c
2
dx
CAAM 452 Spring 2005
Cont (dt estimate)
• The skew symmetry tells us the eigenvalues are on the imaginary axis so we know that if:
3
cdt
2
dx i
, 3
cdt
2
dx i
• is in the stability region we have sufficient conditions for stability of the ODE.
• This is a pretty good estimate, as the eigenvalues
CAAM 452 Spring 2005
Comment on Gerschgorin Bound
>> [v,dxD] = test4th(10); >> >> v v = 0 - 1.36603956377562i
0 - 1.17011114634479i
0 - 0.94222308910582i
0 - 0.62520425034077i
0 - 0.00000000000000i
0 + 0.00000000000000i
0 + 0.62520425034077i
0 + 0.94222308910582i
0 + 1.17011114634479i
0 + 1.36603956377562i
>> 12*dxD ans = 0 8 -1 0 0 0 0 0 1 -8 -8 0 8 -1 0 0 0 0 0 1 1 -8 0 8 -1 0 0 0 0 0 0 1 -8 0 8 -1 0 0 0 0 0 0 1 -8 0 8 -1 0 0 0 0 0 0 1 -8 0 8 -1 0 0 0 0 0 0 1 -8 0 8 -1 0 0 0 0 0 0 1 -8 0 8 -1 -1 0 0 0 0 0 1 -8 0 8 8 -1 0 0 0 0 0 1 -8 0 • The Gerschgorin estimate was pretty good in this case (we suggested eigenvalues in
i*[-1.5,1.5]
)
CAAM 452 Spring 2005
Note on Those Eigenvalues
• What should the eigenvalues of the 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); >> dx = 2*pi/10; >> sort(dxDeigs/dx) ans = >> [dxDeigs, dxD] = test4th(20); >> dx = 2*pi/20; >> >> e = sort(dxDeigs/dx); >> e(1:10) >> [dxDeigs, dxD] = test4th(80); >> dx = 2*pi/80; >> e = sort(dxDeigs/dx); >> >> e(1:10) ans = 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
-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
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
cont
• From Fourier analysis we know that the field u, if sufficiently smooth has a Fourier representation.
• Taking the first derivative we achieve:
u
ˆ
i
2
L x d dx
e i
2
x L
2
L
ˆ
i
2
x L i
2
x
du dx
2
L
L
• Also, each Fourier mode is an eigenfunction of the differentiation operator with eigenvalue
2
L
for all
CAAM 452 Spring 2005
• However, we see an oddity • More about this later.
cont
2
for all
L
• By the completeness of the Fourier modes, we would hope the eigenvalues of the discrete approximation differentation operator would converge to
0, 1, 2, 3,...
• There is the +/-1.6639 in there.
>> [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
Using RK4
a
dt
b
dt
c
dt
d
dt n n n
cdt
Du
n
a
/ 2
cdt
c b
/ 2
cdt
cdt
n
n n
c
a b
/ 2 / 2
u
n
1
u
n
1 6
a
2
b
2
c
d
To avoid confusion with subindices for the vector, we have used super indices to represent the time level.
CAAM 452 Spring 2005
Using AB4
CAAM 452 Spring 2005
That’s a Scheme
• Without much trouble we have created a stable finite difference scheme for the advection equation – given the local truncation errors are small we might expect that the scheme is high order too.
• We will examine the accuracy and conditions for convergence of this kind of finite difference scheme next lecture.
• However, I did steer the choice of spatial difference discretization towards a skew symmetric operator • I will motivate the choices next lecture.
CAAM 452 Spring 2005
Homework 2 (due 02/03/2005)
Q1) Implement the 6 stage, 5 th order Runge-Kutta from: http://www.library.cornell.edu/nr/bookcpdf/c16-2.pdf
Use equations: 16.2.4 , 16.2.5, 16.2.6 with coefficients at the top of table 717 (Cash-Karp coefficients provided on the website) Q1a) Solve
du u
dt
cos 1 0 1 Q1b) Estimate dtmax by concocting a Lipschitz constant for the right hand side of the ODE.
Q1c) Solve 4 times, using (i) dt=dtmax (ii) dt=dtmax/2 (iii) dt=dtmax/5 (iv) dt=dtmax/6At each time step estimate the error (using 16.2.7) and compare with actual error (plot both on a graph, use log scale if necessary).
CAAM 452 Spring 2005