Numerical Methods for Partial Differential Equations

Numerical Methods for Partial Differential Equations CAAM 452 Spring 2005 Lecture 9 Instructor: Tim Warburton

Today • Change of plan • I will go through in detail how to solve homework 3 step by step.

Homework 3 Q1) Build a finite-difference solver for Q1a) use your Cash-Karp Runge-Kutta time integrator from HW2 for time stepping Q1b) use the 4th order central difference in space (periodic domain) Q1c) perform a stability analysis for the time-stepping (based on visual inspection of the CK R-K stability region containing the imaginary axis) Q1d) bound the spectral radius of the spatial operator Q1e) choose a dt well in the stability region Q1f) perform four runs with initial condition (use M=20,40,80,160) and compute maximum error at t=8 Q1g) estimate the accuracy order of the solution. Q1h) extra credit: perform adaptive time-stepping to keep the local truncation error from time stepping bounded by a tolerance.

Q1a – Use Cash-Karp RK for Time-Stepping • We will use a vector version of Cash-Karp: • Where f is a vector valued function – in our case it will be a linear operator acting on a vector argument.

Where the b,c coefficients arehttp://www.library.cornell.edu/nr/bookcpdf/c16-2.pdf Provided in cashkarp.m on the web site. Original paper at: http://portal.acm.org/citation.cfm?id=79507&coll=GUIDE&dl=GUIDE&CFID=38381912&CFTOKEN=60548034

Recall: 4th Order Central Difference Scheme The fourth order central difference derivative acts on any vector and gives the following value for each entry of a result vector: Where the increments on the indexing is done modulo M. Since we will use this operator repeatedly I made a function

Code • Each time step consists of evaluating the 6 intermediate vectors k1,k2,..,k6 • And piecing them together. • Notice – here I set up u at the start as a vector, and delta4 accepts a vector (and dx) as arguments and returns a vector. • i.e. each of k1,k2,..,k6 is a vector.

Alternate Code • I could also have used loops • Note: u is a row vector of length M so the k’s are a matrix of size 6 x M

Cash-Karp • The Cash-Karp Runge-Kutta integrator in our case is: • Notice, we do not need the time component on the right hand side, because our finite-difference right hand side is independent of time.

Linear Stability Analysis • We achieve the linear stability analysis by assuming f is linear in u:

cont • We replace mu*dt with nu

cont • We wish to remove all the intermediate variables and express the scheme in a one-step form like: • Where the multiplier function pi(nu) is a 6th order polynomial in nu. • It is certainly possible to find all the coefficients of this polynomial by hand but a little messy. • We can take a short cut – assuming Cash-Karp is actually a 5th order scheme as claimed we know that the first 6 terms of the polynomial multiplier must be the first 6 terms of the Taylor series for

cont • So we know that: • Where C is to be determined. • However, looking at the definition of the stages we see that there is only one contribution to the 6th order term: • So

Stability Condition • We need to plot the stability region, so we determine the margin of stability by finding nu such that • The curve of points in the (root) complex plane with magnitude 1 can be parameterized in theta by • So for each theta in [0,2pi) we seek the corresponding nu such that:

Matlab roots Command • The matlab roots command can be used to find roots of a polynomial. • If the polynomial is say: • Then construct a vector of coefficients (highest order first): • And invoke: roots(a) • A vector of the roots (if any) of the polynomial will be returned.

Adapting RKabstab.m • So I took the routine from the class web page and modified it slightly • I hard coded it to use the multiplier polynomial for the Cash-Karp • I changed the plotting to use a scatter plot

Absolute Stability for Cash-Karp • This is a classic picture. • What is interesting is that it does not convincingly cover the imaginary axis!.

On The Imaginary Axis Here I used Matlab’s polyval to evaluate themultiplier polynomial for nu on the imaginary axis. Conclusion – the absolute stability region is just to the left of the imaginary axis (not a big issue here)

Q1d) Bound the spectral radius of the derivative operator • All the eigenvalues of the 4th order central difference are on the imaginary axis, with values (recall from Lecture 6): • So the gotcha is that there is no dt such that the eigenvalues*dt are all exactly inside the stability region. • However, we will estimate the largest eigenvalue and make an assumption..

Q1d) continued We can estimate the largest magnitude eigenvalue directly.

Q1e) • Since we are only integrating for 8 periods let’s choose a reasonable maximum |nu|<=1 • So for close to absolute (linear) stability • Implies that we can choose: • i.e. the time step restriction is: • We can further reduce the right hand side of the inequality to reduce marginal growth.

Q1f: Test Runs • I already set up this initial condition for the web demos: • I modified the testrig.m file to call the centraldifference4CKRK54.m file • It looks like reasonable 4th order convergence.. • We examine the error ratios to be more certain T=4

Q1g) Estimate Order of Accuracy of Solution Cool – asymptotically it looks like a fourth order code (in space)

Q1h) Estimating the Local Truncation Error • Given the k vectors we can also construct a lower order (4th order in time) approximation to the updated solution: • Then we can ball park what the time step should be:

Q1h) Code Modify last time step to landon EndTime Form alternative embeddedfourth order approximation Estimate reasonable dt Check to see if we shouldreduce time step.

Q1h) Testing • I started with a fairly large time step and let the code adjust itself These are pretty much unchanged from before.

Q1h) Adaptive Time Stepping(started with dt=100*dx) Q) Why does it take moretime steps for the Lower resolution runs?

Note on the Cash-Karp RK • Check out: • http://www.math.sfu.ca/~cbm/mscthesis/cbm-mscthesis.ps.gz • There is an interesting treatment of RK schemes, with the idea of parameterizing embedded RK schemes (n+1 stage, n’th order) by choosing the coefficient in the multiplier directly. • Increasing this parameter a little will make sure that at least part of the imaginary axis is inside the absolute stability region.

Numerical Methods for Partial Differential Equations