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

Im x Re Recall: Euler-Forward Time Stepping • By several different approximations we arrived at the following approximate ODE time-stepping scheme: • Assuming: • We established the following stability region • i.e. if mu*dt is in the yellow region the approximate solution will be bounded. -1

Recall: Adams Bashford 2 • We went on to derive a similar method which uses • Namely AB2: • This was based on an interpolation of f(u(t)) • We next ask which values of mu*dt is this stable for?

Derivation of AB3 • Suppose we wish to use three previous values of f • As previously we look for a fit through this data, although this time we choose a different path. • Consider dt = 1 (will rescale later) • Will form the coefficients of the interpolant using a Vandermonde matrix and a monomial basis in t • Finally integrate the integrand (easy with this choice of basis).

Simpler Approach • form quadratic interpolant: • Compute coefficients: • Integrate interpolant:

AB3 • We now have an interpolant for dt=1 based on three time levels: • We rescale and shift: • And finally have AB3:

Absolute Stability of AB2 • We again assume: • AB2 becomes: • Tidying up: • To start this method we require:

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 • We will make continuous reference to the stability polynomial of this recurrence relationship (obtained by setting : and canceling excess powers of z )

Note on Solutions of Recurrence Relations • Given an s+1 term recurrence relationship (concisely written as): • With stability polynomial: • which (we suppose) has s distinct roots • Then any solution to the recurrence relation can be written as: • where the cm depend on the s initiating values for the recurrence. • This follows because is a solution, and there is a limited rank of possible solutions because of the linearity of the recurrence in the initial values. Board demo

cont • Since all solutions under these conditions can be written as: • We immediately observe that this is bounded if and only if • See Trefethen p43-46 for details of the case where the stability polynomial has 1 or more roots with multiplicity greater than 1.

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 satisfy and any zero with is a simple root.

Drawing the Stability Region • We need to find the roots of the stability polynomial: • i.e. • Which lie on the margin of stability • The method we adopt is to setand find the value of nu such that r is a root of the stability polynomial.

cont • We can rearrange the root equation to find nu in terms of z: • We can now ask matlab to plot this curve in the complex plane…

Matlab Plot of Margin of Stability for AB2

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.

Zooming Into The All Important Imaginary Axis • Recall, we discounted Euler-Forward as useless for time-stepping the advection equation because it had zero coverage of the imaginary axis. • This picture does not quite tell us ifthe imaginary axis is inside the absolutestability region for AB2, AB3 and AB4. • So we had better check the root condition.

Root Condition for AB2 • We need to check that the roots of: • 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

Always Experiment First • So rather than waste time figuring out where the roots are, why not let matlab do the hardwork:

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.

Consistency • Linear Multistep schemes (like Euler-Forward, AB2,AB3..) can be written as: • It is natural to ask what is the local discretization error (or local truncation error) for the above formula when we substitute • We define the residual operator

Consistency cont • If the linear multistep method is accurate the truncation error should be small. • Definition: if then we say the scheme is p’th order accurate.

Example • Euler-Forward has: • So: • The next step is to estimate this in terms of dt and time derivatives of the solution.

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: • The residual takes various forms: (Lagrange remainder) • http://mathworld.wolfram.com/TaylorsTheorem.html • http://mathworld.wolfram.com/LagrangeRemainder.html

Back to Forward-Euler • We can now express the local truncation error for Forward-Euler in terms of: • expanding about tn we obtain:

cont This indicates that the Euler-Forward methodis “first order in time”. So what does this mean for solution accuracy?

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 .Then there exists a unique differentiable function u(t) that satisfies the initial-value problem: 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.

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 for the method tend to u(0) as • the solution satisfies • uniformly for all The notation small ‘o’ of 1 as dt tends to zero, requires that the solution drops to zero as dt tends to 0.

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.

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: • uniformly for all

Summary • Stability + consistency  convergence • convergence + method accuracy  global solution accuracy

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: N=T/dt

Sketch cont • Represents the actual solution (which we know exists and is unique by the Lipschitz assumptions on f) • Represents acknowledges the impact of using an interpolant for f on the ODE (i.e. the local truncation errors represent the residuals) • Represents the numerical solution u tilde.

Sketch cont • Part 1: we used (i)-(ii) (with T=dt) to estimate L for each time step. • Part 2: we look at (ii)-(iii)

Sketch cont • We can now try to bound each term: Bounded by assumption on Initial conditions (i) Bounded by assumption on consistency. (ii) Estimated by accuracy study. Bounded by assumption on Stability of the Ip interpolationoperator

Euler-Forward Solution Accuracy Consider: • the right hand side function is certainly Lipschitz! • It is certainly p=1 times differentiable • If we keep mu*dt in the absolute stability region (unit circle centered at -1 in the complex plane) then it is stable. • By Taylor’s theorem it is a p=1 order accurate method. • 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

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: AB2 is a 2nd orderacurate method

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 • Thus we characterize the scheme as 2nd order method accurate, and 2nd order solution accurate.

AB3 Convergence Expand all terms in the local truncation error about tn  The AB3 scheme is 3rd order accurate as long as thesolution is sufficiently smooth.

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

Homework 1 cont Q4a) estimate upper bounds for dt to solve the following with AB2,AB3, and AB4 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.

Numerical Methods for Partial Differential Equations