Lecture 2Monday, April 4, 2005Supplementary Reading: Osher and Fedkiw, Sections 3.3 and 3.5;Leveque, Sections 6.7, 8.3, 10.2, 10.4For a reference on Newton polynomial interpolation via divided differ-ence tables, see Heath, Scientific Computing, Section 7.3.3.In the previous lecture we started with a PDE and discretized it to con-struct a numerical method. For convergence, we require that the numericalmethod be both consistent and stable, possibly imposing conditions on thetime ste p (e.g., 4t < α4x) to achieve stability.In this lecture we look at some higher order accurate discretizations.Our goal in using higher order accurate methods is to improve upon consis-tency. The order of accuracy of the method refers to the order of the localtruncation error resulting from the discretization. For example,O(4x) 1st order accurateO(4x2) 2nd order accurateO(4x3) 3nd order accurate......The local truncation error may also have terms such as O(4t4x).1 TVD Runge-KuttaTo achieve higher order accuracy in the temporal discretization, one can useTotal Variation Diminishing (TVD) Runge-Kutta (RK) methods. Thesemethods guarantee that the total variation of the solution does not increase,so that no new extrema are generated. For example, if your solution rep-resents a temperature profile, spurious oscillations caused by the numericalmethod may trigger a chemical reaction in your simulation. Using a TVDmethod would ensure that this situation does not occur. A related conceptis that of Total Variation Bounded (TVB) methods, where the growth in1the total variation is bounded. Both concepts are related to stability. Belowwe consider 1st order, 2nd order, and 3rd order TVD RK.• 1st order, O(4t) errorThe 1st order TVD RK is identical to forward Euler and 1st orderRK. It is given byφn+1− φn4t+~Vn· ∇φn= 0• 2nd order, O(4t2) errorThe 2nd order TVD RK method is also known as 2nd order RK, themidpoint rule, modified Euler, and Heun’s predictor-corrector method.First, an Euler step is taken to advance the solution to time tn+ 4tφn+1− φn4t+~Vn· ∇φn= 0 (1)followed by a second Euler step to advance the solution to time tn+24tφn+2− φn+14t+~Vn+1· ∇φn+1= 0 (2)followed by an averaging stepφn+1=12φn+12φn+2(3)that takes a convex combination of the initial data and the result oftwo Euler steps. The final averaging step produces the se cond orderaccurate TVD (or TVB for HJ ENO and HJ WENO) approximationto φ at time tn+ 4t.• 3rd order, O(4t3) errorFirst, an Euler step is taken to advance the solution to time tn+ 4tφn+1− φn4t+~Vn· ∇φn= 0 (4)followed by a second Euler step to advance the solution to time tn+24tφn+2− φn+14t+~Vn+1· ∇φn+1= 0 (5)2followed by an averaging stepφn+12=34φn+14φn+2(6)that produces an approximation to φ at time tn+124t. Then anotherEuler step is taken to advance the solution to time tn+324tφn+32− φn+124t+~Vn+12· ∇φn+12= 0 (7)followed by a second averaging stepφn+1=13φn+23φn+32(8)that produces a third order accurate approximation to φ at timetn+ 4t. This third order accurate TVD RK method has a stabil-ity region that includes part of the imaginary axis. Thus, a stable(although ill-advised) numerical method results from combining thirdorder accurate TVD RK with central differencing for the spatial dis-cretization.2 Hamilton-Jacobi ENORecall that in the first order accurate upwind differencing, we approximatethe spatial derivative as follows.if ui> 0, use φ−x≈ D−φ =φi−φi−14xif ui< 0, use φ+x≈ D+φ =φi+1−φi4xif ui= 0, then uiφx= 0, so do nothingThis scheme can be improved upon by using a more accurate approxima-tion for φ−xand φ+x. The velocity, u, is still used to decide whether φ−xor φ+xis used, but the approximations for φ−xor φ+xcan be improved significantly.The simplest way to approximate the spatial derivative is to assume thatthe function is piecewise linear. This is what the approximations D+andD−do. However, we could also pass a parabola or higher order polyno-mial through the data points. We can then differentiate the interpolatedpolynomial to get the derivative.The general idea behind ENO is to use a higher order accurate polyno-mial to reconstruct φ and then to differentiate it to get the approximationto φx. Such a polynomial is constructed at each grid point. The key to the3algorithm is to choose the neighboring points for the interpolation so thatwe are not interpolating across steep gradients.Below we describe the implementation of the HJ ENO method in detail.We do Newton polynomial interpolation via a divided difference table. Givena set of n data points, (x1, φ1), (x2, φ2), ..., (xn, φn), the Newton polynomialinterpolating these points has the formPn−1(x) = α0+ α1(x − x1) + . . . + αn−1(x − x1) . . . (x − xn−1)The coefficients, αj, are the entries in the divided difference table, whichwe describe below. We use the notation, Dkito represent the kthdivideddifference at grid point i. For example, D0i= φi. The first few levels of thedivided difference table are constructed as follows.0thD0iφ = φi1stD1i+12φ =D0i+1φ−D0iφ4x2ndD2iφ =D1i+12φ−D1i−12φ24x3rdD3i+12φ =D2i+1φ−D2iφ34x......Note that the zeroth level is defined at grid nodes, the first level is definedmidway in between grid nodes, the third level at grid nodes, etc. Also fromthe above table we can see that D1i+12φ = (D+φ)iand D1i−12φ = (D−φ)i.The divided differences are used to reconstruct a polynomial of the formφ(x) = Q0(x) + Q1(x) + Q2(x) + Q3(x) (9)that can be differentiated and evaluated at xito find (φ−x)iand (φ+x)i. Thatis, we useφx(xi) = Q01(xi) + Q02(xi) + Q03(xi) (10)to define (φ−x)iand (φ+x)i. Note that the constant Q0(x) term vanishes upondifferentiation.To find φ−xwe start with k = i − 1, and to find φ+xwe start with k = i.Then we defineQ1(x) = (D1k+1/2φ)(x − xi) (11)so thatQ01(xi) = D1k+12φ (12)4implying that the contribution from Q01(xi) in equation 10 is the backwarddifference in the case of φ−xand the forward difference in the case of φ+x.In other words, first order accurate polynomial interpolation is exactly firstorder upwinding. Improvements are obtained by including the Q02(xi) andQ03(xi) terms in equation 10 leading to second and third order accuracyrespectively.Looking at the divided difference table and noting that D1k+1/2φ waschosen for first order accuracy, we have two choices for the second orderaccurate correction. We could include