Numerical Methods for Partial Differential Equations

Numerical Methods for Partial Differential Equations CAAM 452 Spring 2005 Lecture 4 1-step time-stepping methods: stability, accuracy Runge-Kutta Methods, Instructor: Tim Warburton

Today • Recall AB stability regions and start up issues • Group analysis of the Leap-Frog scheme • One-step methods • Example Runge-Kutta methods: • Modified Euler • General family of 2nd order RK methods • Heun’s 3rd Order Method • The 4th Order “Runge-Kutta” method • Jameson-Schmidt-Turkel • Linear Absolute Stability Regions for the 2nd order family RK • Global error analysis for general 1-step methods (stops slightly short of a full convergence analysis) • Warning on usefulness of global error estimate • Discussion on AB v. RK • Embedded lower order RK schemes useful for a posteriori error estimates.

Recall: AB2 v. AB3 v. AB4 • These are the margins of absolute stability for the AB methods: • Starting with the yellow AB1 (Euler-Forward) we see that as the order of accuracy goes up the stability region shrinks. • i.e. we see that to use the higher order accurate AB scheme we are required to take more time steps. Q) how many more?

Recall: Requirements Starting Requirements AB1: AB2: AB3: 1 solution level for start 2 solution levels for start 3 solution levels for start

cont • So as we take higher order version of the AB scheme we also need to provide initial values at more and more levels. • For a problem where we do not know the solution at more than the initial condition we may have to: • Use AB1 with small dt to get the second restart level • Use AB2 with small dt to get the third restart level • March on using AB3 started with the three levels achieved above. AB2 AB1 AB3

Recall: Derivation of AB Schemes • The AB schemes were motivated by considering the exactly time integrated ODE: • Which we approximated by using a p’th order polynomial interpolation of the function f

Leap Frog Scheme • We could also have started the integral at: • And used the mid point rule: • Which suggests the leap frog scheme:

Volunteer Exercise • accuracy: what is the local truncation error? 2) stability: what is the manifold of absolute linear stability (try analytically) in the nu=dt*mu plane? 2a) what is the region of absolute linear-stability?

cont 3) How many starting values are required? 4) Do we have convergence? 5) What is the global order of accuracy? 6) When is this a good method?

One Step Methods • Given the difficulties inherent in starting the higher order AB schemes we are encouraged to look for one-step methods which only require to evaluate • i.e. • Euler-Forward is a one-step method: • We will consider the one-step Runge-Kutta methods. • For introductory details see: • “An introduction to numerical analysis”, Suli and Mayers, 12.2 (p317) and on • Trefethen p75- • Gustafsson,Kreiss and Oliger p241-

Runge-Kutta Methods • The Runge-Kutta are a family of one-step methods. • They consist of s stages (i.e. require s evaluations of f) • They will be p’th order accurate, for some p. • They are self starting !!!.

Example Runge-Kutta Method (Modified Euler) • Modified Euler: • Note how we only need one starting value. • We can also reinterpret this through “intermediate” values: • This looks like a half step to approximate the mid-interval u and then a full step. • This is a 2-stage, 2nd order, single step method.

Linear Stability Analysis • As before we assume that f is linear in u and independent of time • The scheme becomes (for some given mu): • Which we simplify (eliminate the uhat variable):

cont • We gather all terms on the right hand side: • [ Note: the bracketed term is exactly the first 3 terms of the Taylor series for exp(dt*mu), more on that later ] • We also note for the numerical solution to be bounded, and the scheme stable, we require:

cont • The stability region is the set of nu=mu*dt in the complex plane such that: • The manifold of marginal stability can be found (as in the linear multistep methods) by fixing the multiplier to be of unit magnitude and looking for the corresponding values of nu which produce this multiplier. • i.e. for each theta find nu such that

cont • We can manually find the roots of this quadratic: • To obtain a parameterized representation of the manifold of marginal stability:

Plotting Stability Region for Modified Euler

Checking Modified Euler at the Imaginary Axis • As before we wish to check how much of the imaginary axis is included inside the region of absolute stability. • Here we plot the real part of the “+” root

Is the Imaginary Axis in the Stability Region ? • We can analytically zoom in by choosing nu=i*alpha (i.e. on the imaginary axis). • We then check the magnitude of the multiplier: • So we know that the only point on the imaginary axis with multiplier magnitude bounded above by 1 is the origin. • Modified Euler is not suitable for the advection equation.

General 2 stage RK family • Consider the four parameter family of RK schemes of the form: • where we will determine the parameters (a,b,alpha,beta) by consideration of accuracy. • [ Euler-Forward is in this family with a=1,b=0

cont • The single step operator in this case is:

cont • We now perform a truncation analysis, similar to that performed for the linear multistep methods. • We will use the following fact:

cont (accuracy) • We expand Phi in terms of powers of dt using the bivariate Taylor’s expansion • where:

cont • We construct the local truncation error as: • Now we choose a,b,alpha,beta to minimize the local truncation error. • Note – we use subindexing to represent partial derivatives.

cont • Consider terms which are the same order in dt in the local truncation error: • Condition 1: • Condition 2: • Under these conditions, the truncation is order 3 so the method is 2nd order accurate. It is not possible to further eliminate the dt^3 terms by adjusting the parameters.

Case: No Explicit t Dependence in f It is easier to generalize to higher order RK in this case when there is no explicit time dependence in f.

Traditional version In terms of intermediate variables: Second Example Runge-Kutta: Heun’s Third Order Formula This is a 3rd order, 3 stage single step explicit Runge-Kutta method.

Again Let’s Check the Stability Region With f=mu*u reduces to a single level recursion with a very familiar multiplier:

Stability of Heun’s 3rd Order Method • Each marginally stable mu*dt is such that the multiplier is of magnitude 1, i.e. • This traces a curve in the nu=mu*dt complex plane. • Since we are short on time we can plot this using Matlab’s roots function…

Stability Region for RK (s=p) rk3 rk2 rk4

Heun’s Method and The Imaginary Axis • This time we consider points on the imaginary axis which are close to the origin: • And this is bounded above by 1 if rk3

Observation on RK linear stability • For the s’th order, s stage RK we see that the stability region grows with increasing s: • Consequently we can take a larger time step (dt) as the order of the RK scheme increase. • On the down side, we require more evaluations of f

Popular 4th Order Runge-Kutta Formula • Four stages: see: http://web.comlab.ox.ac.uk/oucl/work/nick.trefethen/1all.pdf p76 for details of minimum number of stages to achieve p’th order.

Imaginary Axis (again) • With the obvious multiplier we obtain: • For stability we require: rk4

Imaginary Axis Stability Summary 2.83 for the 4th Order “Runge-Kutta” method 1.73 for Heun’s 3rd Order Method 0 for modified Euler

Bounding the Global Error in Terms of the Local Truncation Error Theorem: Consider the general one-step method • and we assume that Phi is Lipschitz continuous with respect to the first argument (with constant ) • i.e. for • Then assuming it follows that

cont Proof: we use the definition of the local truncation error: to construct the error equation: we use the Lipschitz continuity of Phi: tidying:

proof cont

proof summary • We now have the global error estimate: • Broadly speaking this implies that if the local truncation error is h^{p+1} then the error at a given time step will scale as O(h^p): • Convergence follows under restrictions on the ODE which guarantee existance of a unique C1 solution and stable choice of dt.

Warning About Global Error Estimate • It should be noted that the error estimate is of almost zero practical use. • In the full convergence analysis we pick a final time t and we will see that exponential term again. • Convergence is guaranteed but the constant can be extraordinarily large for finite time:

A Posteriori Error Estimate • There are examples of RK methods which have embedded lower order schemes. • i.e. after one full RK time step, for some versions it is possible to use a second set of coefficients to reconstruct a lower order approximation. • Thus we can compute the difference between the two different approximations to estimate the local truncation error committed over the time step. • google: runge kutta embedded • Numerical recipes in C: • http://www.library.cornell.edu/nr/bookcpdf/c16-2.pdf

My Favorite s Stage Runge-Kutta Method • There is an s stage Runge-Kutta method of particular simplicity due to Jameson-Schmidt-Turkel, which is of interest when there is no explicit time dependence for f

RK v. AB • When should we use RK and when should we use AB? rk3 rk2 rk4

Class Cancelled on 02/17/05 • There will be no class on Thursday 02/17/05 • The homework due for that class will be due the following Thursday 02/24/05

Numerical Methods for Partial Differential Equations