PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 5/08 KIProj/Accl-PrelelRC/Dateouejndd INTEGRAL DEFERRED CORRECTION METHODS FOR SCIENTIFIC COMPUTING By Maureen Marilla Morton A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Applied Mathematics 2010 ABSTRACT INTEGRAL DEFERRED CORRECTION METHODS FOR SCIENTIFIC COMPUTING By Maureen Marilla Morton Since high order numerical methods frequently can attain accurate solutions more efficiently than low order methods, we develop and analyze new high order numerical integrators for the time discretization of ordinary and partial differential equations. Our novel methods address some of the issues surrounding high order numerical time integration, such as the difficulty of many popular methods’ construction and han- dling the effects of disparate behaviors produce by different terms in the equations to be solved. We are motivated by the simplicity of how Deferred Correction (DC) methods achieve high order accuracy [72, 27]. DC methods are numerical time integrators that, rather than calculating tedious coefficients for order conditions, instead con- struct high order accurate solutions by iteratively improving a low order preliminary numerical solution. With each iteration, an error equation is solved, the error de- creases, and the order of accuracy increases. Later, DC methods were adjusted to include an integral formulation of the residual, which stabilizes the method. These Spectral Deferred Correction (SDC) methods [25] motivated Integral Deferred Cor- rections (IDC) methods. Typically, SDC methods are limited to increasing the order of accuracy by one with each iteration due to smoothness properties imposed by the gridspacing. However, under mild assumptions, explicit IDC methods allow for any explicit rth order Runge-Kutta (RK) method to be used within each iteration, and then an order of accuracy increase of r is attained after each iteration [18]. We extend these results to the construction of implicit IDC methods that use implicit RK methods, and we prove analogous results for order of convergence. One means of solving equations with disparate parts is by semi—implicit integra- tors, handling a “fast” part implicitly and a “slow” part explicitly. We incorporate additive RK (ARK) integrators into the iterations of IDC methods in order to construct new arbitrary order semi-implicit methods, which we denote IDC-ARK methods. Under mild assumptions, we rigorously establish the order of accuracy, finding that using any rth order ARK method within each iteration gives an order of accuracy increase of 7‘ after each iteration [15]. We apply IDC-ARK methods to sev- eral numerical examples and present preliminary results for adaptive timestepping with IDC-ARK methods. Another means of solving equations with disparate parts is by operator splitting methods. We construct high order splitting methods by employing low order split- ting mcthods within each IDC iteration. We analyze the efficiency of our split IDC methods as compared to high order split methods in [77] and also note that our construction is less tedious. Conservation of mass is proved for split IDC methods with semi-Lagrangian WENO reconstruction applied to the Vlasov—Poisson system. We include numerical results for the application of split IDC methods to constant advection, rotating, and classic plasma physics problems. This is a preliminary, yct significant, step in the development of simple, high order numerical integrators that are designed for solving differential equations that display disparate behaviors. Our results could extend naturally to an asymptotic preserving setting or to other operator splittings. Copyright by MAUREEN MARILLA MORTON 2010 DEDICATION T 0 my father, in loving memory. ACKNOWLEDGMENT If I were to mention everyone who contributed to my success, the list would be longer than my dissertation. Within the limitations of space, I wish to acki'iowledge a few of these amazing people. The first thanks go to my advisor, Dr. Andrew Christlieb, who has infinite faith in the capabilities of his students, and whose support of me is not only academically excellent, but also extraordinarily benevolent and personal. I also am indebted to the members of our research team and our collaborators, in particular to Drs. Benjamin Ong and Jing—Mei Qiu, who have provided valuable insight to my research and indispensable friendship. I wish to thank my comprehensive and defense committee members, Drs. G. Bao, D. Liu, K. Promislow, J. Qian, Y. Wang, and Z. Zhou, for their wisdom, help, and general good humor. Everyone knows that the secretaries and the computer technicians are the most important members of any department, so I cannot fail to thank them as well for making all the details go smoothly. My fellow graduate students were tremendously instrumental in sustaining me through the entire process, and I am particularly grateful for my Chinese friends, who allowed me the privilege of sharing in their yéuyz‘. Most of all, I thank my family: my mother, for teaching me to work hard, to be independent, and how to understand the academic culture; my father, for teaching me that relationships are more important than achievements and for inspiring me with his own unrealized dream of going to college (to become a math teacher); and my brother and sisters, for their delightful encouragement despite their not-so—secret belief that I am insane for studying so much. I am ever grateful to Jesus, my God, for sustaining me here at MSU, just as he promised, and without whose help I would have failed, “For the Lord gives wisdom, and from his mouth come knowledge and understanding” (Proverbs 2:6, NIV). vi TABLE OF CONTENTS List of Tables ................................. x List of Figures ................................ xi Introduction 1 1 Physical Motivation 8 1.1 Chemical Rate Equations ......................... 8 1.2 Convection-Diffusion-Reaction Equations ................ 9 1.3 Vlasov—Poisson Equations ........................ 11 1.4 Hyperbolic Conservation Laws ...................... 12 2 Deferred Correction Methods I 16 2.1 Defect Correction Methods ........................ 16 2.1.1 DC Algorithm for an IVP .................... 16 2.1.2 DC Algorithm for a General Setting ............... 19 2.1.3 Differential Formulation of the DC Algorithm ......... 20 2.1.4 Order Theorems for DC Methods ................ 22 2.1.5 DC Methods and Collocation .................. 23 2.1.6 Continuing Developments in DC Methods ........... 24 2.2 Spectral Deferred Correction Methods ................. 24 2.2.1 Formulation of SDC Methods .................. 26 2.2.2 Order Theorems for SDC With Euler Base Schemes ...... 29 2.2.3 Stability for SDC Methods With Euler Base Schemes ..... 30 2.2.4 Further Developments for SDC Methods ............ 30 2.2.5 Semi-implicit SDC Methods Constructed with Euler Methods 31 Order of Accuracy of SISDC With Euler Base Schemes . . . . 34 Stability of SISDC With Euler Base Schemes .......... 34 2.3 Integral Deferred Correction Methods .................. 35 2.3.1 Explicit IDC-RK Methods .................... 35 2.3.2 Implicit IDC-RK Methods .................... 39 Implicit Runge—Kutta Methods ................. 39 General Formulation of IDC with implicit RK ......... 40 Truncation Error for IDC-BE .................. 43 Truncation Error for Implicit IDC-RK ............. 56 Vii 3 Semi—implicit IDC-ARK Methods 75 3.1 Semi—implicit Integration ......................... 75 3.2 Additive Runge—Kutta Methods ..................... 76 3.3 Semi—implicit IDC—ARK Methods .................... 79 3.3.1 General Formulation of IDC-ARK ................ 79 3.3.2 Example with Second Order ARK ................ 85 3.4 Theoretical Analysis of IDC-ARK Methods ............... 87 3.5 Stability .................................. 98 3.5.1 Stability for Convection-Diffusion Type Problem ........ 99 3.5.2 Stability for Other Types of Problems ............. 101 3.6 Numerical Results ............................. 105 3.6.1 ODE Test Problems ....................... 106 Test Problems ........................... 106 Order of Accuracy ........................ 107 Order Reduction ......................... 108 Efficiency ............................. 112 3.6.2 Advection-Diffusion Example .................. 115 3.7 Concluding Remarks ........................... 117 4 Adaptive IDC Methods 118 4.1 Introduction ................................ 118 4.2 Integral Deferred Correction ....................... 120 4.2.1 Error Equation .......................... 120 4.2.2 Implementation of an IDC Method ............... 121 4.2.3 Adaptive Time Step Control ................... 122 4.3 Numerical Comparisons ......................... 124 4.3.1 Van der Pol Oscillator ...................... 126 4.4 Possible Efficiency Simplifications .................... 132 4.5 Conclusions ................................ 133 5 Split IDC Methods 134 5.1 Introduction and Motivation ....................... 134 5.1.1 General Motivation: Modeling and Simulating Plasmas . . . . 134 5.1.2 Specific Motivation: Improved Vlasov Solvers ......... 136 5.2 Method of Lines .............................. 141 5.3 WENO ................................... 142 5.3.1 Reconstruction .......................... 144 5.3.2 ENO ................................ 145 5.3.3 WENO ............................... 146 1-D Scalar Case .......................... 147 l-D Systems ............................ 152 5.4 Semi—Lagrangian Methods ........................ 155 5.5 Splitting Methods ............................. 157 5.6 IDC with Splitting Methods ....................... 160 viii 5.6.1 Overview of IDC methods .................... 160 5.6.2 Prediction Loop .......................... 162 5.6.3 Splitting for Correction ...................... 163 5.7 Efficiency Analysis ............................ 168 5.8 Conservation of Mass ........................... 173 5.9 Numerical Results ............................. 177 5.9.1 Constant Advection ........................ 178 5.9.2 Rotation .............................. 180 5.9.3 Drifting Rotation ......................... 181 5.9.4 Vlasov—Poisson .......................... 186 Two Stream Instability ...................... 187 Landau Damping ......................... 190 5.10 Conclusions ................................ 197 I‘hture Work 199 6.1 Asymptotic Preserving Methods ..................... 199 6.2 IDC with AP-ARK Theory ........................ 200 6.2.1 Moving Towards an AP IDC Proof ............... 200 6.3 AP IDC Methods for P1 Equations ................... 205 6.3.1 The P1 Problem and the Splitting ................ 205 6.3.2 Asymptotic Limit ......................... 207 6.3.3 Splitting Error .......................... 208 6.3.4 Split Scheme with one IDC Correction ............. 209 6.3.5 Analytic Solutions of Splitting Parts .............. 212 6.3.6 Anticipated Numerical Tests ................... 214 Bibliography .................................................... 215 LIST OF TABLES 4.1 Test characteristics from [59] for Van der Pol’s oscillator solved with various semi-implicit methods when 6 = 10’6. Note here that all characteristics are for steps of size At, including for the IDC methods, since we accept or reject each At step, not each 6t step. ....... 130 5.1 Number of split solves required for one step 'r of each splitting method. Yoshida’s methods are from [77], and r : At. IDC-S2 refers to split IDC methods constructed with second order splitting, and IDC-Sl refers to split IDC methods constructed with first order splitting. 7' 2 61$ for IDC methods. ........................ 170 3.1 3.2 3.3 3.4 3.5 LIST OF FIGURES Stability regions for 3rd, 6th, 9th, and 12th order IDC constructed using (a) forward and backward Euler with 2, 5, 8, and 11 correction loops (and 3, 6, 9, 12 uniform quadrature nodes), respectively, and (b) 3rd order ARK3KC with 0, 1, 2, and 3 correction loops (and 3, 6, 9, and 12 uniform quadrature nodes), respectively. The crosses denote the 3rd order stability regions, diamonds for 6th order, circles for 9th order, and triangles for 12th order. The stability regions are scaled by the number of implicit function evaluations ((#stages)(M +1)M, where FEBE has 1 stage, ARK3KC has 4 stages). ........... 100 Stability region for fourth order IDC with ARK2A1 in prediction and correction loops. A N is plotted, consistent with Definition 3.5.3, and a = 7r / 2 ................................... 103 Stability regions for IDC with 2nd order DIRK (from [2]) in prediction and correction loops. The labels are the number of correction loops, and are placed outside of the stability regions. From 0, 1,. . . ,4 cor- rection loops, the IDC methods have corresponding order of accuracy 2, 4, . . . , 10. Stability regions are standard, according to Definition 2.2.2.105 Convergence study of absolute error at T = 4 vs H, for 3rd, 6th, and 9th order IDC constructed using 3rd order ARK3KC with 0, 1, and 2 correction loops (and 3, 6, and 9 quadrature nodes), respectively. The order of accuracy is clearly seen, as the dotted reference lines (with slopes of 3, 6, 9) indicate. .................... 109 Convergence study of absolute error at T = 4 vs H, for 4th and 8th order IDC constructed using 4th order ARK4KC with 0 and 1 correction loops (and 4 and 8 quadrature nodes), respectively. The order of accuracy is clearly seen in (a), as the dotted reference lines (with slopes of 4, 8) indicate. ...................... 110 xi 3.6 3.7 3.8 3.9 4.1 4.2 4.3 Convergence study of absolute error at T = 4 vs H, for 3rd, 6th, and 9th order IDC constructed using 3rd order ARK3KC with 0, 1, and 2 correction loops (and 3, 6, and 9 quadrature nodes), respectively. The dotted reference lines (with slopes of 3, 6, 9) indicate the expected order of the methods, but note order reduction. ............ 111 Efficiency studies of # implicit function evaluations vs absolute error at T = 4, using (a) ARK methods of second (2A1 [56], 2ARS [3]), third (3BHR1 [8], 3KC [45]), and fourth (4A2 [56], 4KC [45]) orders (no IDC). (b) for 6th and 9th order IDC constructed using 3rd order ARK3KC and 3rd order ARK3BHR1 with 1 and 2 correction loops (and 6 and 9 quadrature nodes), respectively. ............. 113 Efficiency study of # implicit function evaluations vs absolute error at T = 4, for 3rd, 6th, and 9th order IDC constructed using 3rd order ARK3KC with 0, 1, and 2 correction loops (and 3, 6, and 9 quadrature nodes), respectively). .................... 114 CF L study showing absolute error at T = 0.1 vs H = At = 0.5Asr. FFT is used for the spatial discretization, and time-stepping is done with 3rd and 6th order IDC constructed with 3rd order ARK3KC with 0 and 1 correction loops of IDC, respectively, where the advection term is treated explicitly and the diffusion term is treated implicitly. This choice of methods ensures that the error is dominated by the time-stepping. Since the expected order of accuracy in time is seen, as the dotted reference lines (with slopes of 3, 6) indicate, it appears that the CF L condition is At S cAa: ................... 116 Preliminary results comparing (4.1a) an adaptive ARK method (em- bedded order 4(3) from [45]) and (4.1b) an adaptive IDC method con- structed with a 3rd order ARK method from [45] (both recursively implemented). Both plots show solution to Van der Pol’s oscillator (c = 1) for the same recursion error tolerance, 10—5. ......... 127 Efficiency study comparing semi—implicit ARK and IDC—ARK meth- ods used to solve Van der Pol’s oscillator. Plots show # imp f vs atol, as described in the text. ...................... 128 Efficiency study comparing semi-implicit ARK and IDC-ARK meth- ods used to solve Van der Pol’s oscillator. Plots show # imp f vs atol, as described in the text. ...................... 129 xii 5.1 5.2 5.3 5.4 5.5 5.6 Convergence study for (5.149) using (5.1a) IDC constructed with lst order split methods, each split equation is solved via forward Euler, combined with 5th order WENO. The order of accuracy is clear, as reference lines (with slopes of 1, 2, 3, 4) indicate. (5.1b) IDC constructed with 2nd order split methods, each split equation is solved via 2nd order Runge-Kutta, combined with 9th order WENO. The order of accuracy is clear, as reference lines (with slopes of 2, 4, 6) indicate. In both figures, the error is in the norm ll + 100, compared to a reference solution computed on a finer time mesh. ........ 179 Convergence study for (5.150). (5.2a) IDC constructed with lst or- der split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th order WENO. The order of accuracy is clear, as reference lines (with slopes of 2, 3, 4) indicate. (5.2b) IDC constructed with 2nd order split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th order WENO. The order of accuracy is clear for 2nd and 4th order, but unclear for 6th order, as reference lines (with slopes of 2, 4, 6) indi— cate. In both figures, the error is in the norm ll + ’00. comparing successive solutions ............................. 182 The order of accuracy is clear, as reference lines (with slopes of 1, 2) indicate. In all figures, the error is in the norm 11 + [00- ....... 184 The order of accuracy is clear, as reference lines (with slopes of 2, 3) indicate. In all figures, the error is in the norm [1 + 100- ....... 185 Convergence study for (5.88) with ICs as in (5.155). IDC constructed with lst order split methods, each split equation is solved in a semi- Lagrangian setting with conservative 5th order WENO. (5.5a) Here three IDC nodes are used, with 0, 1, and 2 correction loops. The order of accuracy is clear, as reference lines (with slopes of 1, 2, 3) indicate. (5.5b) Here four IDC nodes are used, with 0, 1, 2, and 3 correction loops. The order of accuracy is clear, as reference lines (with slopes of 1, 2, 3, 4) indicate. In both figures, the error is in the norm ll + loo, comparing successive solutions. ............. 189 Solution for (5.88) with ICs as in (5.155). IDC constructed with lst order split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th order WENO, N3; = N0 = 400, At = 1/300. ................................... 190 xiii 5.7 5.8 5.9 5.10 5.11 5.12 Physical quantities for (5.88) with ICs as in (5.155). IDC constructed with lst order split methods, each split equation is solved in a semi— Lagrangian setting with conservative 5th order WENO, N3; = 64, NU = 128, At = 1/40, and M + 1 = 4. The labels lst, 2nd, 3rd, and 4th in the legends denote first order splitting (prediction only), 2nd order split IDC (prediction, 1 correction), 3rd order split IDC (prediction, 2 corrections), and 4th order split IDC (prediction, 3 corrections) methods, resp. ....................... Convergence study for (5.88) with ICs as in (5.156). In Figure 5.8a, a = 0.01, and in Figure 5.8b, a = 0.5. IDC constructed with lst order split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th order WENO. Here four IDC nodes are used, with 0, 1, 2, and 3 correction loops. The order of accuracy is clear, as reference lines (with slopes of 1, 2, 3, 4) indicate. The error is in the norm [1 + 100, comparing successive solutions. ........ Solution for (5.88) with ICs as in (5.156), a = 0.01. IDC constructed with lst order split methods, each split equation is solved in a semi- Lagrangian setting with conservative 5th order WENO, N3; = NU = 400. (5.9a) At = 1/300. (5.9b) At = 1/300. .............. Physical quantities for (5.88) with ICs as in (5.156), a = 0.01. IDC constructed with lst order split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th order WENO, N3; = 64, N0 = 128, and M + 1 = 4. The labels lst, 2nd, 3rd, and 4th in the legends denote first order splitting (prediction only, At = 1/40), 2nd order split IDC (prediction, 1 correction, At = 1/40), 3rd order split IDC (prediction, 2 corrections, At = 1/80), and 4th order split IDC (prediction, 3 corrections, At = 1 / 80) methods, resp. Electric field for (5.88) with ICs as in (5.156) and a = 0.01. IDC constructed with lst order split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th order WENO, N3; = 64, NU = 128, M +1 2 4. For lst order: At = 1/80, 2nd order: At = 1/80, 3rd order: At = 1/80, 4th order: At = 1/80. Solution for (5.88) with ICs as in (5.156), a = 0.5. IDC constructed with lst order split methods, each split equation is solved in a semi- Lagrangian setting with conservative 5th order WENO, N3; = NU = 191 192 193 . 194 . 195 400, M +1 2 4. (5.12a) At = 1/80. (5.12b) At = 1/300. ....... 195 xiv 5.13 Physical quantities for (5.88) with ICs as in (5.156), a = 0.5. IDC constructed with lst order split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th order WENO, N3; = 64, NU = 128, and M + 1 = 4. The labels lst, 2nd, 3rd, and 4th in the legends denote first order splitting (prediction only, At = 1/40), 2nd order split IDC (prediction, 1 correction, At = 1/40), 3rd order split IDC (prediction, 2 corrections, At = 1 / 80), and 4th order split IDC (prediction, 3 corrections, At = 1 / 80) methods, resp. Electric field for (5.88) with ICs as in (5.156), a = 0.5. IDC con- structed with lst order split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th order WENO, N3; = 64, N0 = 128, 111 + 1 = 4. For lst order: At = 1/80, 2nd order: At = 1/80, 3rd order: At = 1/80, 4th order: At = 1/80. XV . 196 . 197 Introduction Throughout scientific history, experiments, physical and mathematical models, and, more recently, scientific computing have been partners in generating new knowledge. Although all are vital components of discovery, we focus on the role that scientific computing plays. In many cases, experiments may not be feasible due to expense or physical constraints, and models may be too complex to find analytic solutions. In these situations, further insight may be found through computationally approx- imating solutions to models, such as those models given by ordinary differential equations (ODEs), partial differential equations (PDEs), differential algebraic equa- tions (DAEs), or others. The daily changes in technology and scientific discoveries means that the need for innovative numerical methods always exists. In particular, it is desirable to obtain accurate solutions at a low computational cost and in real time. Frequently, higher order numerical methods can attain such accuracy more ef- ficiently than lower order methods, although care must be taken to maintain certain properties; e.g., for some physical problems, one must be careful that a numerical method conserves mass when the physical theory asserts that mass is preserved. In this dissertation, we focus on higher order numerical integrators for the time discretization of ODEs and PDEs. Some of the issues surrounding high order numerical time integration include: difficulty of the method’s construction (for example, the order conditions for con- structing higher order Runge-Kutta methods and the number of coefficients one 1 must find for some general high order splitting methods increase exponentially with the order of the method [37, 77]); the equations exhibit multi—scale behaviors and / or different operators and terms in the equations produce disparate but coupled (fre quently nonlinear) behaviors, leading to severe timestep restrictions or difficulties in numerical implementation; and treatment of boundary conditions at high order. In this dissertation, we contribute to the work that tackles some of these issues by further developing Integral Deferred Correction (IDC) methods, described below. One way to overcome the first issue of a. complicated construction of high order methods was introduced and later expanded upon in articles such as [78, 65, 28, 27, 72]. This class of methods, called Deferred Correction (DC) methods (or Defect Correction), involves taking a prediction step to form an approximate solution using a numerical method of choice, then iteratively forming an error equation, solving for the error using any numerical method of choice, and updating the approximate solution. The main principle is that, at each iteration, the order of accuracy of the DC method increases by the order of the numerical method used at each iteration. The prediction step can be considered one iteration, and each step that solves the error equation (also called the correction step) can be considered one additional iteration. The beauty of DC methods is their simplicity in constructing higher order methods. One only needs to use low order integrators within the prediction and correction steps, and the iterative process attains arbitrary high order without needing to consider complicated order conditions. A key component of DC methods is the treatment of the residual within the error equation. For simplicity, suppose we are solving an initial value problem (IVP) y’(t) = f(t,y), y(0) = 310- (1) The “differential form” of the residual, 7‘, is 7'0) = 77(6) — [(1, 72(0), where 7) is a provisional numerical solution to (1). Approximating the derivative, n’(t), within the DC correction step is numerically unstable (see, e.g., [32]). More recently, Greengard and Rokhlin introduced a variant of DC methods, called Spec- tral Deferred Correction (SDC) methods, by incorporating the integral formulation of the residual, ft <>d — - (tn—[tn ())d OTT T—TIA T] 0 T,T)T T, to stabilize the correction step [25]. This development sparked renewed interest in DC methods. Several methods influenced by SDC methods include Krylov Deferred Correction methods, applied to ODEs and DAEs [41, 42]; semi-implicit and multi- implicit SDC methods (SISDC and MISDC, respectively) applied to ODEs whose right hand sides include stiff terms (i.e., terms with disparate sizes) [60, 47, 9]; and Integral Deferred Correction (IDC) methods, applied to ODEs and PDEs, with a special choice of nodes and integrators used at each step [18, 17], including adjust- ments to allow parallel in time implementation of IDC methods [16]. Typically, SDC methods are limited to increasing the order of accuracy by one after each correction step due to certain smoothness properties that are imposed by the gridspacing. Some works that adjusted SDC methods to incorporate integrators other than Euler methods in the prediction and correction loops include [48, 47, 18, 17]. In particular, [18, 17] present IDC methods that allow for any explicit RK methods of arbitrary order to be used within the prediction and correction loops, and that allow for an order of accuracy increase greater than one after each correction 3 loop. They rigorously established under certain assumptions on the gridspacing and the smoothness of the exact solution that, when an rth order explicit RK method is used to predict and correct the numerical solution, the order of accuracy of the IDC method increases by r for each prediction and correction loop. Thus fewer correction loops are required to attain a certain order of accuracy, when compared to SDC methods. They also found that such high order IDC methods’ stability and efficiency compared favorably to RK methods alone. In this dissertation, we extend the construction and theory of IDC methods using explicit RK integrators to the construction of arbitrary order IDC methods that use implicit RK integrators and a rigorous presentation of analogous results to the order of convergence theory. The second issue of high order numerical integrators, that of handling equations with distinct parts, multi—scale, and / or nonlinear behaviors, is quite broad. We focus on methods that handle different terms in the equations via separate means. First consider semi-implicit, or implicit-explicit, integrators. Typically, a semi-implicit integrator is intended to solve an equation of the form 31,0) = F(t,y) + C(tay), (2) with appropriate initial or boundary conditions. For simplicity of presentation, we suppress any non-time dependence that may be in F and G, which could be functions or operators. Suppose F contains only nonstiff (not particularly large) terms and G contains any stiff (could be large) terms. Semi-implicit integrators can be used to solve the stiff part implicitly and the nonstiff part explicitly, thereby gaining the benefit of the implicit method for the stiffness but potentially saving computational effort by handling some terms explicitly. Popular integrators include additive Runge- Kutta (ARK) methods and integrators constructed from linear multi-step methods, such as Adams or BDF methods. [60, 47] introduce semi-implicit SDC methods that 4 may attain orders of accuracy higher than ARK or BDF methods, but in [60], the increase in accuracy at each correction loop is limited to lst order, while [47] only presents specific examples of a 2nd order ARK and a 2nd order BDF method that may be successfully used in the correction loops. In this dissertation, we expand the results in [47, 18], and our IDC with implicit RK methods to construct a general formulation that clearly presents the ability to incorporate any order ARK scheme into the IDC framework. We denote the new construction as IDC-ARK methods. We rigorously establish under certain assumptions that, when an rth order ARK method is used to predict and correct the numerical solution, the order of accuracy of the IDC method increases by r for each prediction and correction loop. We also include numerical examples of IDC-ARK methods applied to IVPs that substantiate the theory and to an advection-diffusion equation, verifying that the semi-implicit structure allows for an improved CFL condition. To further handle multi-scale situations, especially equations whose solutions contain initial, inner, or boundary layers, we also present preliminary results for adaptive implementation of IDC-ARK methods. Since IDC methods form an approximation to the error at each time step, they fit naturally into an adaptive setting, where the size of the time step is adjusted as necessary according to some estimation of the size of the error at that step. Equations of the form (2) may also be solved via operator splitting methods. Splitting methods may be useful when, for example, y depends on more than one spatial variable, say on (271,332), in addition to time, t, but F depends only on .731 and G depends only on 1172. In this case, we may approximate the solution by successively solving the following two subproblems, 8w = F, Bty = G', which now are simple one-dimensional problems rather than a more complicated 5 two-dimensional problem. Sometimes also an Operator splitting may be chosen such that the subproblems are linear although the original problem (2) is nonlinear. Many high order splitting methods are costly both to construct and to implement [77], although some reasonable methods have been developed in [49, 50], including high order splitting methods constructed via DC methods and applied to Vlasov- Maxwell systems. We wished to seek an improved high order splitting method by using the IDC framework instead. We present these novel split IDC methods in this dissertation. They incorporate simple first order splitting methods and second order Strang splitting methods within the prediction and correction loops to obtain arbitrary order splitting methods. The construction is specified for application to the Vlasov-Poisson equations, which are frequently used to describe the behavior of particles in some plasma physics settings. We show that a new formulation of the error equation in the correction step is required for proper application to the Vlasov- Poisson system. As required by the physical theory, we also prove that split IDC methods for Vlasov equations conserve mass when conservative semi-Lagrangian WENO reconstruction is employed in conjunction with the IDC splitting method. The third issue that causes difficulties in numerical time integration, treating boundary conditions (BCs) in a way that maintains the high order of the numerical method, will not be discussed in this dissertation. However, this area is essential if the discussed methods are to be applied to realistic problems. In particular, this field is mainly unexplored for deferred correction methods. The application of split DC methods to the Vlasov-Maxwell system in [49, 50] only utilizes periodic BCs. Essentially, in this dissertation, we also only apply IDC methods to differential equa- tions with periodic BCs. Motivated by the high order treatment of BCs presented by LeVeque in [52, 53], we expect that equivalent or alternative means of handling BCs could be designed to take advantage of IDC’s unique framework of extending low order methods to high order methods. This design may circumvent the extensive 6 approximations of higher order derivatives that are necessary in LeVeque’s work. Chapter 1 briefly introduces several equations that describe certain physical pro- cesses that motivate some of the numerical issues we wish to address via IDC meth- ods. Chapter 2 presents a history and literature review of DC methods and more recent developments. It also includes a description of our new arbitrary order IDC methods that incorporate implicit RK methods in their construction and a rigorous proof of the theoretical order results for such IDC methods. Chapter 3 lays out the construction of our new semi-implicit IDC-ARK methods, order of accuracy theoret- ical results, how the semi-implicit proof differs from the implicit proof in Chapter 2, and application to numerical examples such as Van der Pol’s oscillator, an initial layer problem, and an advection-diffusion equation. Chapter 4 is a short study on the implementation of IDC methods in an adaptive setting, with some preliminary results for adaptive IDC-ARK methods. Chapter 5 shows the construction of new high order split IDC methods, an analysis of efficiency as compared to high order splitting methods in [77], a proof of conservation of mass, and numerical results for its application to constant advection, rotating, and Vlasov-Poisson equations. Chapter 6 suggests directions for future research, such as extending asymptotic preserving methods to the IDC framework. Chapter 1 Physical Motivation Many physical problems involve both fast and slow events, or can be modeled by equations that have operators whose numerical solutions can be simplified by op- erator splitting. Modeling these situations may result in equations that contain widely separated time scales or equations whose parts may be computed sepa- rately. These equations include but are not limited to: chemical rate equations, convection-diffusion equations, Vlasov-Poisson equations, and hyperbolic conserva- tion laws with stiff relaxation. We present some of these problems as motivation for studying semi—implicit numerical integrators that handle two different time scales and splitting methods that allow for simpler numerical computations, although we have not necessarily tested all the problems described below. 1.1 Chemical Rate Equations Chemical reactions are often modeled by rate equations that are formed by using the mass action law. An example of a system of rate equations is given by the ROBER 8 problem, which describes an autocatalytic reaction [59]. The system is d ayl = —0.04y1+104y2y3. (1.1) (1 fit” = 0.04y1 — 104y2y3 — 3 - 107313, (1.2) d 7 2 _ = .1 , (1,313 3 0 yg, (13) y1(0) = 1, y2(0) = 0, 313(0) = 0, (1.4) where the y,- represent the concentrations of the chemical species involved, and the coefficients on the right hand sides are the rate constants, which describe the rate of certain component reactions. The widely disparate rate constants contribute to the stiff nature of this problem, and the solution to this system changes quickly initially (although it varies smoothly as time progresses). Typically, there are far more species involved in a chemical reaction, thus rendering the system of rate equations much larger than the ROBER problem; e.g., the chemical Akzo Nobel problem or the pollution problem, both also in [59]. However, the ROBER problem sufficiently demonstrates that widely separated time scales may arise in chemical rate equations. 1.2 Convection-Diffusion—Reaction Equations Multiscale effects also play a role in the solution of convection-diffusion-reaction equations. Typically, the diffusion term dominates the other terms due to the dis- cretization of the spatial derivatives, and the reaction term is significantly smaller 9 than the other terms. An example of such an equation is given in [45], in the form p p’u ] a rm 6 m2+p 960 (960+Plu -pYz- . ...y, _ 0 0 4 62a 87:2 0 + 2 2 3y. + 8 T 4 Bu ('9 nos , . Am+'§(7fi) +835< i=1pDzhflfi7‘) 0 DOBQY- “,2 - ,0 281:5 _ ' ' (1.6) This system is a one-dimensional version of the simplified, gas-phase, multicompo— nent, compressible NavierStokes equations with chemical reaction. ()1, A, pDi) are the transport coefficients, where p is molecular viscosity, A is thermal conductivity, p is fluid density, and Di is the effective Fickian diffusion coefficient. The species is designated by i = 1, . . . ,nes, where ncs is the number of chemical species, the fluid velocity is represented by u, T is temperature, Y,- are the species mass fraction, p is pressure, e0 is total specific internal energy, h is the partial specific enthalpy of species i, w,- is the reaction rate of species 2'. The first term on the right hand side is the convection term, the second term is the diffusion term, and the last term is the reaction term. In addition to the differences in sizes of the terms that could arise based on the spatial derivatives, there are also potential separations of time scales based on the sizes of the transport coefficients and the reaction rates. 10 1.3 Vlasov—Poisson Equations Many natural and industrial processes, such as solar weather and integrated circuits manufacturing, involve plasmas. A plasma is sometimes called the fourth state of matter, since adding sufficient energy to a solid produces a liquid, adding energy to a liquid leads to a gas, and increasing the energy yet again results in a plasma. A plasma can also be defined as an electrically neutral collection of charged particles. The sizes, speeds of, and forces acting on or by these charged particles may vary by several orders of magnitude. Hence it is of interest to consider computational issues in models of plasmas. As a special case, the electron motion in a one-dimensional cold plasma with a stationary ion background may be modeled by the following Vlasov-Poisson equations, where the physical constants have been normalized to one. where f (:1:,v,t) is the electron probability density function, x is space, 1) is the velocity, t is time, E is the electric field, d) is the potential, and p is the charge density [54]. If we choose to rewrite (1.7) in a Lagrangian framework, from the point of view of an electron within the flow (as opposed to an Eulerian framework, which has a fixed spatial grid), and discretize the problem in (.r, v)-phase space, we obtain N ODE systems with varying stiffnesscs: cage) = v,(i), (1.10) v’.(t.) = E(t,:r,-), 2': 1,. . .,N, (1.11) where the electric field on the right hand side of (1.11) contains contributions from the electrons, ions, and free space [14]. The velocities and forces may vary sig- nificantly in magnitude between different values of i and among the terms on the right hand side of (1.11) for each i. Alternatively, one may use a semi-Lagrangian framework to solve (1.7) for the distribution, in which case splitting methods can be applied. Firrther details are presented in Chapter 5. 1.4 Hyperbolic Conservation Laws The Vlasov equation is one type of a class of problems known as hyperbolic con- servation laws. Much of the following information about such problems is compiled from [66, 51, 70]. For simplicity of explanation, we explain hyperbolic conservation laws via the general scalar one-dimensional case, solving “t + ffulsr = 0, (1-12) where f (u) is the flux. Corresponding to the physics involved, such an equation satisfies conservation of mass. Consider the case where 21(13, t) = p(:r, t), the density, with flux f(p(:r,t)), over the interval Ii = [:13 1,.7:_+1], where the flow is to the z— 2 z 2 right: f(/)(¢L‘,l)) MN) f(p(I,t)) l l (1.13) $Z—% V x2+$ I The integral form of conservation of mass over the interval 1,- is: mass at [n+1 = / p(r,tn’+l)d.r (1-14) 2 tn+1 = fl.- po, mm + [n f(p(xi_%, z» — we”; mm (1.15) = mass at tn + flow in — flow out. (1.16) We sketch how the integral form (1.14) leads to the differential form of the conser- vation law. Looking at a small piece in space 651: and in time (it, we have (p($»t"+1)- pt‘r- £70) 51‘ = - (f(p($.+1~t)) - f(p($._1=t))) 515. (1-17) 2 2 z 2 which, upon rearranging, gives 1' ,n+1 — .13, Tl p( ,t (lst p( ,t )= 1 (f(p(:ti+1at))_f(p(x- 17t)))' (1.18) (5.1: From this equation, we deduce the differential form of the conservation law Pt + ffplcr = 0, (1-19) which in a more general case is given by (1.12). In a physical setting, the integral form (1.14) is more natural and meaningful. An understanding of some basic analytic properties is needed for considering a numerical solution. We describe the solutions of two types of equations in terms of characteristics, the lines along which solutions are constant. First we consider a linear problem, such as the linear advection-diffusion equation, at + 111: = 0, u(.r, 0) = u0(:r). (1.20) 13 The solution is u(:z:, t) = -u0(;r — t). This solution is constant along characteristics described by griff- = 1, since (in. 2) dr (‘9 ('9 (9 dt (at dt 8:13) (at 8:13) ( ) In this situation, the characteristics will not cross. For a nonlinear problem, the sit- uation is much more complicated, and characteristics may cross after a certain time, resulting in a discontinuity of the solution. As an example, we consider Burger’s equation “t + (223) = 0, u(:L‘,0) = u0(a:). (1.22) a: The solution is constant (u = 11.0) along characteristics described by 931% = u, but since the slope of the characteristics is determined by %, then in many cases, the characteristics will cross after a certain time; e.g., when u0(:z:) = sin :22. Where the characteristics come together, a shock forms. This discontinuity means the differential form of the conservation law does not work, so the integral form is needed. Only a classical strong solution works for the differential equation, but the integral equation can be solved by a weak solution. A weak solution is not unique, but a weak solution satisfying a property called entropy is called the entropy solution, which is unique. For example, in the Riemann problem, the initial conditions are given by at, a: < 512* note) = , (1.23) u-r, a: >:r* where “l and ur are constant functions. If the entropy condition is satisfied, then “l > ur means the solution will form a shock, and “l < ur means the solution will form a rarefaction. 14 The numerical solution of hyperbolic conservation laws is described in Sec- tion 5.3, with an emphasis on a numerical method called WENO, which is designed to solve problems with piecewise smooth solutions that contain discontinuities. Sometimes hyperbolic conservation laws have a stiff relaxation term R(u), such as diffusion, on the right hand side: Ht“. + E)g:f('rt) = 72(11), (1.24) where R(u) is quite large, thus contributing to the multiscale nature of the problem. This type of term shows up in models such as shallow water, granular gas, and traffic flow equations [64] . 15 Chapter 2 Deferred Correction Methods 2.1 Defect Correction Methods Deferred Correction (DC) methods were developed in the 19608 and 1970S. They can be considered methods to solve a system of ODEs (that may arise within an algorithm solving certain PDEs), or in the more general sense, to solve an equation of the form F y = 0. Our current interest lies in the ODE and PDE application. 2.1.1 DC Algorithm for an IVP In particular, suppose we are solving the IVP y’U) = f(t,y). téloaTl, 31(0) = yo- (2.1) DC (sometimes called Iterated Defect Correction) involves taking a prediction step to form an approximate solution via an arbitrary numerical method (e.g. RK), then iteratively forming an error equation, solving for the error, and updating the approx- imate solution. It is of interest to note that the approximate solution approaches 16 the collocation solution to the IVP as the error is improved at each step [28]. We solve the IVP (2.1) on a grid 0=t0 0. If an arbitrary { explicit or implicit) RK scheme of order p (p g M ) is used, and if f satisfies suitable diflerentiability conditions, then the approximate solution after the kth correction satisfies 77%] _ yam) : 0(h’fn'l71(p(k+1),fll)) f0?" h __) 0, (2.24) where y is the exact solution. [27] has a similar version of Theorem 2.1.1, with variant assumptions: Theorem 2.1.2. If the discretization method is of order q, if the approximation k—l] 7)] from the previous DC step{s) gives an 0(hr) error, and if certain smoothness requirements are met, then the result nfkl of a DC step satisfies "[kl — lIhy = 0(h7nin(r+q),J)) for h —> 0, (2.25) 22 where I'lhy is the projection of the exact solution y onto the grid, provided that the maximum attainable order J is larger than p and q. Moreover, further DC steps are possible if J > r + q. In essence, Theorem 2.1.2 tells us that Theorem 2.1.1 holds when p is different at each DC iteration, i.c. when the discretization method has a different order at each iteration. Note the importance of equidistant gridspacing. The same results do not generally hold for nonequidistant grids. 2.1.5 DC Methods and Collocation As aforementioned, the result of the DC solution approaches a collocation solution, or rather the fixed point of the DC method is a collocation solution to the problem (2.13). This collocation solution is determined by the choice of interpolating polyno- mial, for example, and also determines the maximum attainable order of accuracy of the DC method. The following theorem from [28], where the discretization method used at each DC step is backward Euler, illustrates one such claim. Theorem 2.1.3. 0n [0, H], 77* = (yo, 77f, . . . , 777”) is afixed point of the DC method with base scheme of backward Euler if and only if d*(tm) =0, m=0,1,...,M (2.26) where M) = Wm) - rattle». (2.27) and 77(*)(t) interpolates 77* (inc n(*)(t) is the solution of the corresponding collo- cation scheme). This theorem holds for all nodes (2.3), not necessarily equispaced; however, equidistant nodes seem to converge to 77* more readily. 23 2.1.6 Continuing Developments in DC Methods In [6], Stetter is quoted: “The scheme [DC] is simple; the essential point is the defini- tion of the defect [residual].” The treatment of the residual does in fact appear to be the major factor in successful adaptations of DC methods. In the years immediately following DC’s conception, some difficulties became apparent. One of these trou- bles was a lower than expected order of accuracy with nonequidistant nodes (recall that Theorems 2.1.1 and 2.1.2 required equispacing). However, a benefit of some nonequidistant nodes is a higher order collocation solution with fewer gridpoints (e. g. Gaussian quadrature). In 2004, Auzinger et a1 [5] rewrote the residual in several different ways to obtain similar order of accuracy results for nonequispaced nodes. For one modified DC method, the residual and the embedded discretization scheme are both approximated at the unequally spaced gridpoints. For another method, the residual was interpolated at the unequally spaced nodes while the discretization method at each iteration was carried out at equidistant gridpoints. Auzinger et al also provided some numerical examples where their modified DC methods were successfully (and unsuccessfully) applied to stiff IVPs [6]. Another difficulty of early DC methods is the possible instability of the method even when the lower order base scheme is stable, which may be due to the dis- cretization used for forming the residual [40, 35]. The next section suggests how the discretization of the residual could affect the stability and discusses a more stable formulation introduced in 2000. 2.2 Spectral Deferred Correction Methods Recall in Section 2.1.3, the formation of the residual involved the derivative of a polynomial interpolant. This derivative is typically approximated by what is known 24 as spectral differentiation. A spectral differentiation matrix is a linear map DM‘ f(.rm) -> f’(.’L'm) , (2.28) L d h d where D M can be found by representing f by a polynomial interpolant of degree M (or any truncated series expansion), and differentiating the interpolant, so that the i jth entry of ’D M is: d Dij = Ely-(allay (2.29) where M t_tk cj(l)= H t--t' (2.30) k=0,k7éj J k The problem is that D M is increasingly ill-conditioned as M increases [32]. Spectral integration, on the other hand, appears to be inherently more stable. A spectral integration matrix is a linear map If”: ffCEm) _’ fxmf(x)dx , (2.31) _ 1 _ _ Z 1 where I M also can be found by representing f by a truncated series expansion. For example, representing f by its Lagrange polynomial interpolant of degree M, the 25 i jth entry of 1M is: Ilj Z/t Cj(l)dt, (2.32) where cj(t) is as defined in (2.30) above. Perceiving spectral integration’s benefit, Greengard introduced spectral integra- tion to replace spectral differentiation within spectral methods [32]. He proposed recasting the differential equation as an integral equation, where the solution can be stably recovered by integration. For example, he found that for the differentiation matrix PM, the process of differentiation via Chebyshev series can amplify errors by a factor proportional to M 2, whereas for the integration matrix IM, the process of integration amplifies errors by a factor of less than 2.4. In 2000 [25], spectral integration rather than differentiation was applied to the residual in DC, and the new method was designated as Spectral Deferred Correction (SDC). Apparently SDC does not have the instability issues which plague classical DC [60], [40], [35], [32], and it may even cause improved stability in some cases where explicit RK methods are used as the base discretization schemes [18], [17]. 2.2.1 Formulation of SDC Methods SDC methods are deferred correction methods which use spectral integration to approximate the residual and (in the original formulation) Euler timestepping to update the prediction and correction steps [25]. The details of the algorithm may be seen by applying an SDC method to the scalar IVP (2.1). The time interval, [0, T], is discretized into subintervals tn=l1 C and any natural numbers M, Kloop’ SDC with either Euler method as the base discretization scheme with M +1 submeshpoints ( on the interval [0, H] ) and K I oop correction steps converges to the exact solution y = (y0,y1, . . . , ym, . . . , yM) with order of accuracy 0(},mm(M+11Kloop+1))_ Xia, Xu, and Shu [76] also prove a similar theorem. The use of Gauss-Lobatto quadrature nodes in the submeshpoints results in an improved order of accuracy when SDC is combined with one of the Euler methods. With M + 1 Gauss-Lobatto min(2M,K +1)) nodes, the order is 0(h loop [41]. 29 2.2.3 Stability for SDC Methods With Euler Base Schemes Stability regions for SDC methods can be calculated in the traditional sense. Definition 2.2.2. Region of absolute stability. When solving the Dahlquist problem y, = Ay for t E [0, 1] and A E C (2.45) mm = 1, (2.46) the amplification factor Am(A) (or the stability function R(z) for z = Ah 6 C, h is one timestep) is defined by Am(A) = R(::) -—'- 771, (2.47) where 771 is the numerical solution at time t = 1 (ort = h ) The region of absolute stability S is S(Ah) = 3(2) = {A E C s.t. |Am(A)| g 1} = {z E C s.t. |R(z)| S 1}. (2.48) SDC with a base scheme of forward Euler has an interesting property that the size of the stability region increases as the order of the method increases [25]. All of the SDC methods with a base scheme of backward Euler are A(a)-stable. They are not L-stable, but the limit lim|A|—>ooAm(’\) is less than %. Unfortunately, some of the imaginary axis appears to be lost as the order increases [25]. 2.2.4 Further Developments for SDC Methods Several variants of SDC methods have been developed recently. Huang, Jia, and Minion developed an accelerated version of SDC methods for ODE initial value prob- lems and for differential algebraic equation systems using Newton-Krylov methods 30 [41, 42]. In these so-called Krylov Deferred Correction methods, the original SDC method acts as a preconditioner for the system. An iteration of SDC preconditions the equation, a solution for the error is computed by a N ewton-Krylov method, and then the process is repeated, thus accelerating SDC’s convergence to the collocation solution. Also, Minion developed semi-implicit SDC methods which treat the stiff and nonstiff parts of a stiff ODE separately via an implicit and explicit method, respectively [60]. Christlieb, Qiu, and Ong experimented with higher order explicit RK methods as the base scheme in Integral Deferred Correction methods (IDC), which are motivated by SDC methods [18, 17] (see also Section 2.3.1). Also, higher order splitting methods for PDES were constructed by using ADI methods within the differential DC framework, including applications to Vlasov-Maxwell systems with periodic boundary conditions [46, 49, 50]. 2.2.5 Semi-implicit SDC Methods Constructed with Euler Methods When a stiff ODE can be split into a nonstiff part and a stiff part, it is advantageous to treat the nonstiff term explicitly and the stiff term implicitly in order to save computational cost but still handle the stiffness effectively. See the introduction to Chapter 3 for more details on semi—implicit methods. Consider the initial value problem y'(’») = my) = [5012(1) + fNU/ay), '56 [0’7”]. 31(0) = 310- (2.49) The stiff terms are contained in f S(t, y), and f N(t, y) contains only nonstiff terms. Minion [60] developed semi-implicit SDC (SISDC) methods to handle an IVP of the form (2.49). In each predicition and correction step of an SDC method, backward Euler was applied to the stiff term while forward Euler was applied to the 31 nonstiff term. I.e., 1. Prediction step: compute an approximate solution to (2.49) over the grid points (2.86), 0 0 0 0 T][0] = (77!) ],ng 1,. . . Win], . . . ‘77[i})’ (2.50) which is a first order approximation to the exact solution y=(3/01311a-~-~.-33/ma,y/1/[) (2.51) when we apply forward and backward Euler as follows: [0] _ [0] 0 0 nm+1 — 77m + hfS(tm+1v77,[n]+1)+ thUma 77in])- (2.52) 2. Correction step: use the error function to improve the accuracy of the approximate solution at each iteration. For k = 1 to K loop’ where K loop is the number of correction steps: (a) Denote the exact error function from the previous step as e = y(t) — wit—”<0, (2.53) where y(t) is the exact solution and "(k—1)“) is an M m degree poly- nomial interpolating nlk _1l. Notice that the error function e(k—1)(t) is not a polynomial in general. (b) Compute the residual function, Arum : (aw—19W) — f) elk—1)M— f5> — fN(t.n(""1)(t))- 32 (c) Compute the numerical error vector, dik] = (6gk],. . . , (57%,. . . £5611), dis- cretizing the integral form of the error equation, (A's-1w): fs(t~n(k—1)(t) + e(k—1>>— fS> (2.54) + fN(t.n(k"1)(t) + e) — fN(t,n(k_1)(t)) (2.55) _.), Q(k-1)(0) = 0. (2.78) where Q(k_1)(t) = e(k_1)(t) + E(k_1)(t), (2.79) Flt, 6(0) = N, W) + 6(1))- “1277(0), (2-80) t E(k“1)(t) = / Alf—llama (2.81) t0 The rkth order numerical approximation to eri-l] = Q(k—1l(tm) is are] m . Then compute the numerical error vector as 5W = Qlkl _ Elk—ll, (2.82) where Eli—1] = E(k—1l(lm). We apply a pk-stage rkth order implicit RK method to (2.78), giving, 42 for i=1,2,...,pk and m=0,1....,1t’I—1: ,, pk tm cih . k.- = Fem +c.h. all?” +h Z “2'ij - f, + €(k‘1ldr>. i=1 0 Qlkl+1__9[kil + h j::k1bjk 5119] +1 _ —Q[k]+1 _E[k-ll(tm+1). (2.83) In fact, we use It A. tm+1 -_ din]+l :anLl -/to £0” 1)(T)al‘r k:— ~ Qlfi]+1_ nlfi+ll _ 310 — j: aj(t.]f_17[ 1]) a where we have approximated the integral by interpolatory quadrature formulas, as in [25], and multiplied nlk—ll by an interpolation matrix to obtain its value at intermediate stage times tm + cjh, cj 7t 0, 1 [l8]. 4. Update the numerical solution nlkl = n[k_1l + blkl. Truncation Error for IDC-BE An understanding of the proof of the order of accuracy of implicit IDC-RK methods is aided by walking through the case of backward Euler (BE) before the case of higher order implicit RK as the IDC base scheme. The goal is to show that, under _ certain conditions, IDC with a first order BE method used in the prediction and correction steps (IDC—BE) has a local truncation error of 0(hK100p+2). To show this, we need some definitions and propositions about the degree of smoothness and differentiation of functions and of discrete data. Some of these definitions and 43 properties are given below in Section 2.3.2, but also see [18] for further details. The 8 following notation will be used: 3? = partial derivative and 5% = total derivative, e.g. -2léfl@ — [it + a—y at (2.84) d Several analytical and numerical preliminaries are needed to analyze IDC meth— ods. The smoothness of discrete data sets will be established, analog to the smooth- ness of functions; this idea of smoothness is used to analyze the error vectors. Let v5(t) be a function for t E [0,T], and without loss of generality, consider the first interval, [0, H], of the grid (2.71). Denote the corresponding discrete data set, (£4?) = {6640). . . . . (ma/2M», (2.85) where tm are the uniform quadrature nodes given by (2.72). Definition 2.3.5. (smoothness of a function ) A function w(t), t E [0, T], possesses . ‘ s, , ._ as , - __ 0 degrees of smoothness if Hd illloo .— “53de is bounded for s — 0,1,2,...,o, where [llél’lloo := ma‘xte[O,T] [Ii/’20)]. Definition 2.3.6. { discrete difi‘erentiation): Consider the discrete data set, (t, 6:), defined in (2.85), and denote LM as the usual the Lagrange interpolant, an M th degree polynomial that interpolates (t,i,b). An sth degree discrete difierentiation is a linear mapping that maps 1b = (won/21, . ..,ibM) into (1737/), where (dsib)m = 83M aw tion €2le '2 .DS ' 1,5, ’thBTC b5 6 R(M+1)X(M+1) and (133)qu = $671“) (t, wlltztm' This linear mapping can be represented by a matrix multiplica- ltztm, m,n=0,...,M. Given a distribution of quadrature nodes on [0, 1], the differentiation matrices, . 0 1 . . . . . D5 ’ l, s = 1, ..., M have constant entries. If this distribution of quadrature nodes 44 is rescaled from [0,1] to [0, H], then the corresponding differentiation matrices are ~ .01 ~ s~0,1 D1=717D]’]andD3-—-(T1I)D]9 1. Definition 2.3.7. The (S, 00) Sobolev norm of the discrete data set (t: if?) is defined as (231.0 ll} 3 OO 0 [3.06 3: go 0' l =ZIIDs~r 00 3:0 —> where dsil' = Id - (b is the identity matrix operating on 1)). Definition 2.3.8. (smoothness of a discrete date set): A discrete data set, (2.85), possesses a (o S M) degrees of smoothness if “I; is bounded as h —r 0. “3.00 We emphasize that smoothness is a property of discrete data sets in the limit ——-—> as h —> 0. We also impose a g M, because dot/2 :- 6, for o > M. See [18] for a detailed discussion. Example 2.3.9. (A discrete data set with only one degree of smoothness). Consider the discrete data set (Er/7) = {(0.0), (g g) , (g 521-) , (531— 5;) ,(H,O)}. The first derivative ~ ’ 4 10 10 4 d : —— —— __ _ 1w ( 3, 3 ,0, 3 73)? is bounded independent of H, while the second derivative 57/:— 373 16 112 16 272 '2" 3H’ 3H’ 3H’ 3H’3H ’ 45 _,—+ is unbounded as H —» 0. Therefore, (t, it) has one, and only one degree of smoothness in the discrete sense. For further details, definitions, and properties of smoothness in the discrete and continuous sense, we refer the reader to [18]. Now consider the local truncation error of SDC methods with a first order back- ward Euler method in the prediction and correction steps. Analagous to Theorem 4.1 in [18], we have the following theorem: Theorem 2.3.10. Let y(t) be the solution to the IVP (21) and y(t) has at least S (S 2 M + 2) degrees of smoothness in the continuous sense. Consider one time interval of an IDC method with t E [0, H] and uniformly distributed quadrature points O=toh<5-1> (5-1) 5 : y7'l1+1+§yrn+1_m+ (S— 1)! y7n+1 +O(h )+RTTL+1 (2.89) gives S—l (+1) . (-l) 2 ‘ (i) S . S Rm+1 = :2 Try/m +1 + on. )= Tm+1+ ow ) (290) Z: Clearly rm+1 ~ (9(h2) (since y has S degrees of smoothness). Now consider the error: em+i = ym+1 - 77m+1 (2-91) = 31m — 72m + h (f(tm+1vym+l) - f(tm+1i "mi-1)) + Tm+1+ 0013) (2.92) i em + hum+1 + Tm+1 + 0(hS) (2.93) 47 Taylor expanding about ym+1 gives um+l = fum+1vym+1l ' fum+1v71m+ll (2'94) = (jam—Hi ym+1l _ fum+li ym+1 _ em+lll (2'95) S-2 ' ' _1 2+1 . at _ = Z (—_,i)F—(em+ll15;if(tm+lvym+ll + 0((em+l)S 1) (2'96) i=1 ' 1. Now prove part 1. of Lemma 2.3.11 by induction on m, the index of gridpoints on [0, H]. c When m = 0, 60 = 0 ~ 0(112). 0 Assume em ~ 0072), and show em+1 ~ (901.2). Then cm+1 —_— em + hum+1 + Tm+1 + 0(hS), or: 6771+]. — hum+1 = 8n), + 7‘m+1 + 0(h5) (2.97) = on?) + on?) + 0(65) = on?) (2.98) Note that (_1)t+ I: at em+l _ hum+l = em+1_ h( 2 Ti+0<54> and ”aloe ~ (9(h). We are now ready to prove that E has M degrees of smoothness by in- duction. Again, since “aloe ~ 0(h), Ehas at least zero degrees of smoothness in the discrete sense. Assume that E has 3 S M — 1 degrees of smooth- ness. We will show that (dlé) has 3 degrees of smoothness, from which we can conclude that E has (8 + 1) degrees of smoothness. Since fyi has (S — i — 1) degrees of smoothness in the continuous sense, the vector form fyi = [fyi(t0,y0), . . . (”Away/1.1)] has (S — i -— 1) degrees of smoothness in the discrete sense. Consequently, hi_1fyz- (vector) has (S — 2) degrees of smoothness, which implies that E has min (S — 2, 3) degrees of smoothness. Similarly, 723 has (S — 2) degrees of smoothness in the discrete sense. Hence (6115) has 3 degrees of smoothness, which implies that E has (3 + 1) degrees of smoothness. Since this argument holds for S 2 M + 2, we can conclude that E has M degrees of smoothness. Cl Lemma 2.3.12. (Correction steps). Let y(t) be the solution to the IVP (21) and y(t) has at least S (S 2 M + 2) degrees of smoothness in the continuous sense. Consider one time interval of an IDC-BE method with t E [0, H] and uniformly distributed quadrature points (2.86). Suppose in the (k — 1)th (k S M) correction step, elk—1] ~ 0(hk+1) and the rescaled error vector Elk—ll = fielk—ll (now O(h)) has M — k + 1 degrees of smoothness in the discrete sense. Then, after another correction step: 50 1. The updated error vector elk] satisfies Hemlloo ~ out”). (2103) 2. The rescaled error vector Elk—ll = 1 elk—1] has M — k. degrees of smooth— ft—E ness in the discrete sense. Proof. We integrate (e(k-1))’(t) and use that result and blkl to obtain an equation for elk]. Then we prove 1. and 2. Recall the error equation (2.78) at the end of the (k — 1)th correction step, and integrate it over [tm, tm+1]. 451;? = (2)5-” + t:r:n+1F(t,e(k_1)(t))dt— ”Zn-H €(k-1)(t)dt = (3)1571] +/tm+1 F(t,e(k’1)(i))dt —/t F(t, e(k_1)(t))dt tm+1 m — C(k—1)(l)dt tm ._ t - = (2% 1] +/ m+1 F(t,e(k_1)(t))dt _ Utm“ F(t,e(k_1)(t))dt + (-—1)h. r(tm+1,e[f,;]]) 2 2h (1 [ls—1] + (—1) fig (t.m+1,67n+1)+ . . . M+1 M M+1 h d [k-ll M+2 +(——1) (M+ 1)!dil‘(P(tm+1’em+1H001 ) — /tm+1€(k—1)(,)d, tm 51 elk—1l_lk—1l+A§1(_1)i+1Ed’i—1F(elk—ll)+0(hIW+2) em+l _ em . i! dti—l tm+1,e m+1 t —/ m+16(k—1)(t)dt (2.104) Irm Using blkl = (6):]? . . , (5%). . . 2%)) to denote the numerical error vector obtained from the backward Euler scheme at the kth correction step gives: + "1+1”nm+l +6m+1) m+1177m+1l _ ( tm+1€(k—l)(t)dt+0(h/U+2)) , (2105) ,1111_,111+,, (711 11—11,11-11_,(, 11— 11) tm where the last term is due to the spectral integration used to numerically integrate C-(k 1l(1 1) at M + 1 gridpoints. Now subtract equation (2.105) from equation (2.104) in order to obtain em m+1' 111 _ 11— 11_,111 em+1_ em+1 m+1 M+1 ' '.—1 67—11 1 hid7 et—h 1) -/ m+1:(:k1“1)(1)dt—6],€] tm [k— 1] [k— 1]) [k— 1]) — h(f(t "1+1 77m+1 + 6m+1)f(tm+1’nm+1)) + (/t "1+1 C(k_1)(t)dt + O(hM+2)). eifiLl = Bib] +”(WWW 1 "i:+i] +eih+i]lf(tm+1 "i7,:+i] +6ih+i]ll IV 1171—1 d2 , + Z(—1)i(,-+1)1_F(tm+ldti Gilly?) + 0(hAI+2)1 (2.106) 52 e[';1+1=e[f,] + hughl+1+ rlk+ll + 0( 1“”), (2.107) where A. k k 1 E7111: 1111111111111 — 11 141 11.41] + 61.19) k = f(tm+l-ym+ll _ f(tm+11ym+1 — ein1+ll AS'_2 1 (-1 1)2+ 32 1 S—1 = —-f(t 1...+1.ym+111e£,,L11 +O<(111: 1929(111— f(1 461— elk—11(1)) — +1 £212”: flf(ty(1))(e(k ”(11111101101 1>(1))S-1). 3y] Since Elk-’1] has (M — k + 1) degrees of smoothness in the discrete sense, then adZEIk—llu) ~ 011(1) by Proposition 3.19 from [18]. Then ——(——.f<1.99915114191111“ +783;f<1,y<1>>% 3 degrees of smoothness. Since (dl'élkl) has 3 degrees of smoothness, we can conclude that Elk] has (M + 1 — 11') degrees of smoothness. 55 [3 Proof. (of Theorem 2.3.10). Use induction with respect to k, the index of correction steps in the SDC method. 0 By Lemma 2.3.11, Theorem 2.3.10 is valid for k = 0 and the rescaled error vector E’lOl has M degrees of smoothness in the discrete sense. 0 Assume Theorem 2.3.10 is valid for (k — 1) correction loops and the rescaled error vector Elk_1l has (M — k + 1) degrees of smoothness. Then by Lemma 2.3.12, Theorem 2.3.10 is also valid for k correction loops and the rescaled error vector Elk] has (M -— k) degrees of smoothness in the discrete sense. Truncation Error for Implicit IDC-RK Here we provide the local truncation error estimates for IDC methods which use high order implicit RK methods in the prediction and correction steps. Theorem 2.3.13. Let y(t) be the solution to the IVP (2.1) and y(t) has at least S (S 2 M + 2) degrees of smoothness in the continuous sense. Consider one time interval of an IDC method with t E [0, H] and uniformly distributed quadrature points { 2. 86 ) Suppose an rolh order implicit RK method is used in the prediction step and (r1, r2,. . . ’TKloopflh order implicit RK methods are used in K1001) correction steps. Let 3k = 229:0 rj. If SKloop g M + 1, then the local truncation error is of order 3 +1 001 K1001) ). The proof of Theorem 2.3.13 follows from Lemmas 2.3.14 and 2.3.17, below, for the prediction and correction steps, respectively. The following lemma discusses the local truncation error and smoothness of the error for the prediction loop of implicit IDC-RK methods. 56 Lemma 2.3.14. (Prediction step): Let y(t) be the solution to the I VP (2.1) and y(t) has at least S (S 2 M+2) degrees of smoothness in the continuous sense. Con- sider one time interval of an IDC method with t E [0, H] and uniformly distributed quadrature points (2.86). Let nlOl = (ngolmgoly . .,11.,[,9,l,. . ”$811) be the numerical solution on the quadrature points ( 2.86), obtained using an rOth order implicit RK method (see Definition 2.3.4) at the prediction step. Then: 1. The error vector em = y — nlOl satisfies 14011110 ~ 0am“) (2108) and 2. The rescaled error vector ElOl = fielo} has min(S —r0, M) degrees of smooth- ness in the discrete sense. First consider some useful facts, in the exact solution space and the numerical solution space, needed to prove the lemma. Proposition 2.3.15. The exact solution of the I VP (2.1), with S degrees of smooth- ness, satisfies S— 1 hJ _ S 11.....1— — 1111 + 2 711010111) + 001 > (2.109) j: 1 The term y(J J) in (2.109) can be expressed as 0) 1953.1, 1‘.) = wl wnzwfn , )wif wf 2.110 11 ( ,< y<> 1: aqut 1,, Q11U(fqt11y < 1 where f, =Q—7-L, the a l, 1,12, are constant coefficients, and w), tqtyqy atqtayqy qut qy qt w f qt, qt, and n,- are nonnegative integers. 57 Equation (2.110) is called ‘elementary differentiation’ in any standard numerical ODE book, e.g. [37]. and is more easily understood by computing the first few derivatives of y: 11%) = r; 1191(1) = 11+ 111; 1131(1) = f11+ 2 [1111 + 1111!? + 111/1131; ... (2.111) The equation can be easily proved by induction on j. Proposition 2.3.16. Assume that y(t) in I VP (2.1) has S > r degrees of smooth- ness. Suppose that an r I” order RK method is used to solve {21) If h is the size of each sub-interval [tm, tm+1l1 within [0, H], then 7. 77771+1 = 77m+hE1(tm177m)+' ' '+:Er(tm1777n)+Rr(tm,72m)+0(hS), (2.112) where the function E1‘(t, y) is defined in (2.110), and the remainder term S—l Rr(t,r}) = Z N (2.113) '-r+1 ul wni w 2 X31 n- ft 2.1—[Mfg w’Um Umlf wf(trn 77m) 1 q1 (ll ”‘1 Zan __ t (1y) (uh/2191 yr. 311.1— y . wi 111112.111 fn has constant coefi‘lczents 13 1 1 .determined by the formulation of the RK qut" qy (11" method. Proof. The right hand side of (2.112) comes from Taylor expanding the numerical solution of an rth order RK method. Note that the first r + 1 terms coincide with the Taylor expansion of the exact solution (2.109). (2.113) can be proved by induction. [:1 58 Proof. (of Lemma 2.3.14). The proof is identical to that given in [18], but we repeat it here for completeness. The superscript [0] is dropped since there is no ambiguity. 1. First, we prove part 1. of the lemma. From equations (2.109) and (2.112), we see that a Taylor expansion of the error, em+1 = ym+1 —- nm+1, about t = tm gives ern+1 = em + um + 7'1,'m. _ 72m + OMS): where Um = Z 3(Ezltm1yml— Ei(tm177m)l i=1 . 7‘0 'S—l—i ._1 . . hz —1 J (8 )J 8JE- _- 2 F E: ( l ! m -?,(tm,ym) +O((6m)S 7.), - - J By] 221 3:1 S—l hi () 71,777, = 2 7y 2 (Lm) i=T0+1 I S—l . wl...wn.ur T217”: 2 h] 2 fiqldl qgirlfni .-= , , . 7 j r0+1 QE+Q§SJ—1 Jt ”U t ni _.w, 1. “’f m, 11310121193) (15m, Imlf (tmfl? ) N ow we bound llelol “00 by induction. By definition, e0 = 0, so obviously e0 ~ (hm—H). Assume that em ~ (MOT—1). Since um ~ (h7'0'f2),7'1’m ~ (hm—H), and r27", ~ (hm—H), then we have em+1 ~ (hr0+1), which completes the in- ductive proof. 2. Now we prove part 2. of the lemma, the smoothness of the rescaled error vector. We again use an inductive approach, but with respect to s, the degree 59 of smoothness. Note that a divided difference approximation to the derivative of the rescaled error vector gives a 41; . _ (11511,, = w = rm + rm, + 172,", + 0015-70 1), where “i _ Um UILL _ h70+1$ 0 1 1+1~( 1)— 1 . 1) J h 0 J 3);}. ~ :21! SE i-( 1 -(tm1ym)(€m)] 7' j1= J- ByJ +0(h5"’"0—1), ~ 7'1m S—l hi_r0—1 - 1,771 hr0+1 2" y 7n 3 i=T0+1 I b—T0—2 ~ T2m ' /w1~-wn-wf _ 1 _ E J E 2 h -:0 . (lg/qt "’93; qt J qt+q§J <1—1 . (1 q. q.- )wiem. 1111111“? (1m. 11111). 2:1 t t3! y Now we prove that E has rnin(S — r0, M) degrees of smoothness by induction. Since ”aloe ~ (h) is bounded, E has at least zero degrees of smoothness in the discrete sense. Assume that E has .3 < rnin(S —- r0, M) degrees of smoothness in the discrete sense. We will prove that (d1?) has 3 degrees of smoothness, from which we can conclude that E has (3 + 1) degrees of smoothness in the discrete sense. From a similar argument as in the proof of Lemma 2.3.11, a has 3 degrees of smoothness in the discrete sense. Assuming that y(t) has S degrees of smoothness, F1 has (S — r0 — 1) degrees of smoothness. Since 5 has 3 degrees of smoothness, 6 has rnin(s + r0, M) degrees of smoothness, 60 and 77 = y —- e has rnin(s + r0, M) degrees of smoothness. Thus, 13 has rnin(s + r0. S — r0 — 1, M), at least 3, degrees of smoothness in the discrete sense. Therefore, 5 has (5 + 1) degrees of smoothness, and we can conclude that E has rnin(S — r0, M) degrees of smoothness in the discrete sense. [:1 The following lemma discusses the local truncation error and smoothness of the error for the correction loop of implicit IDC-RK methods. Lemma 2.3.17. (Correction step): Let y(t) be the solution to the IVP (2.1) and y(t) has at least S (S 2 M + 2) degrees of smoothness in the continuous sense. Consider one time interval of an IDC method with t E [0, H] and uniformly dis- tributed quadrature points ( 2.86). Suppose elk-1] ~ 0(hsk“1+1) and the rescaled _S£_felk‘1l has M + 1 — sk_1 degrees of smoothness in the kth error vector elk— ll = discrete sense. Then, after the correction step, provided an rkth order implicit RK method is used and k S K1001), 1. The error vector elk] satisfies Helklnoo ~ 001%“) (2.114) and 2. The rescaled error vector Elk] = fielkl has M+1— 5k degrees of smoothness t in the discrete sense. Further, recall the error equation (2.78) >’>— 111101-111» _1(11 : F(t,e(k—1)(t)) — 105—”(1). 61 It can be written as (62(k—1))'(t)= molt—”<11 — law—”<01 2 G>. (2.115) where Qlk-lln) = e(k-1)(1) + E(k—1)(1) and E(k’1)(t) = f6 1(k-1)(1)111. The analysis for the correction steps in this section will rely on this form of the error equation (2.115). Define the rescaled quantities as (it—11(1): LJ’HM) h5k—l ’ ~ k—l _ 1 k—l d M) — mid >111, (311—11,,611—110» = EefiG(o). (2.117) In the correction steps of the IDC method, a p—stage, rth order implicit RK method (2.68) applied to (2.117) gives p 11,: C(k-llno +e,11, egg-1] + h 2 11,511,) for 1' = 1,2,...,p, (2.118) 1:1 elixir—11 1. p11. 1111 l—QO + _ 1.7 l' l 3:1 where the Greek letter fl denotes the rescaled solution in the numerical space. In the actual implementation, we discretize (2.78) as mentioned above in (2.83). 62 Proposition 2.3.18. Ife[k_1] ~ Omsk—1+1), then o(k-1)(i) ~ oasis-1+1) and C(k“1)(t, Q(k-1)(t)) ~ oasis—1+1). (2.120) Proposition 2.3.19. Suppose elk—1] ~ Omsk—1+1) and the rescaled error vector Elk—1] = Egg—_Telk—ll has M + 1 — Sk—l degrees of smoothness in the discrete sense. Then 606—1)“) = 7513—1Q(k—1)(t) satisfies L _ dl “’(k— (HQ 1)(t)~0h(1)iflgM+1—sk_1, foralltE[0,H]. (2.121) Proposition 2.3.20. (“jag-1)“) in (2.117) satisfies AI— —Sk_ 1 ~k—1 ~k—1+l_i ~ k—l ~k—1 M+1_ i=1 or equivalently, A’I—Sk_1 hi k 1 1 1 ' k— 1 1 €£n+11=ein l+h°k—1 Z HG(_1)(tm,Qin-_ l)+0(hM+1) i=1 t —/m+1€(k—l)(r)dr (2.122) tm where ~e_ ~.d2'—1u---ww w.-- G§_11)(t,Q)=;t;:-1-7 1 ”l. fezflfl ”(‘k 1) We“ 1h“? q1q1“q zq _. qi~qz~) 91‘. y t i— tth Q ~ (9h(1), fori g lief — Sk—l' “’lUn' .w As before.7w1 1 "i 7,1. are constant coefficients and wi wf, qtn, q‘~, and "i are qut" (1y qt Q 63 nonnegatiue integers. The proofs of Propostions 2.3.18, 2.3.19, and 2.3.20 are exactly as given in [18]. The proof of the following proposition is the essential difference between the proofs of the truncation error of the explicit and the implicit cases of an RK method used in the prediction and correction loops of IDC. Proposition 2.3.21. The Taylor series for 5(k—1)(t0 + h) = fleas—DUO + [H hsi— 151 sufiiciently smooth error function Elk—1)“) and assuming the exact solution y(t) h) and for in (2.83) coincide up to and including the term hr for a has 5' degrees of smoothness. Proof. First we claim that risk—113,: k, + OMS—1) for a11i=1,2,...,p, (2.123) and, in particular, setting P P .h—Sk :.,1..: we claim that we can write, for all i = 1, 2, . . . ,p and for n = 2,3, . . . , S — 1, hSk-lh—k- S— 23— 2 - p :2 =h" Z Z 5:2 Z Z Cilia 2'11 ll=1l2=1 ln=1j1=1j2= 1 177:1 n . . 5~— ~. _ . S—l H ClinJ'm—la‘Jm—Umih A 1kJn kJn)+O(h )’ (2124) 111:2 _ 1 8’ z z— ,— where Clj — NESTGUO + th, Q0) 21):], A], ’UB; We prove (2.124), and hence (2.123), by induction on n. 64 o n = 2: Using the definitions (2.116) of the rescaled quantities above and Taylor expanding G about Q0 in the y—variable, we obtain hSk—IEZ- — kt -_-_ C(to + C'h,Q0 + At) — C(lo + Cz'h,Q0 ‘1' Bi) 1 81111 ll S—l =:z_:11—1571—y0 (t0+czh./1Qo)(( ) —(Bz-) )+O(h ) 1 l l —u -) 1 = (A,- — 13,-)ZU121A.1 ’1‘ i Factoring (Aflll — (B [33—1, we obtain hsk—IE—k12 5:20111-(24, B.,;_)+O(h5 1) 131:1 _ hs , s—1 -52? 011,11 iaazjl(hk—1hj—l k1)+0(h ), 11= J1=1 and repeating the process for hsk“1l;j1 — kjl’ _ ‘15:? Cllih 2a “(.71 1: J1:1 [221 C h :0 (ask—1 k k)+O(hS—1)- ’2J1 “J1J2( J2 kJ2 12:1 W) S——2$2P #122: Z Z Z C11iCl2J1aiJ1aJ1J2(ma—1722—kJ'2) l1=112=1j1=—1j2=1 3—1 + 0(h ). (2.125) Thus risk-1'13,- — 1e,- ~ (9012). o n+1: Assume (2.124) holds for n < S — 1. Following the same process of Taylor 65 expanding, etc. as in equation (2.125) in the n = 2 case for hsk ‘17:?" — kjn’ we see that hSk—IEZ- — k2- . S—2 10 = h" 2 Z Cilia'ijl ll,...,ln:1]1,...,jn=l n . . 8:. ~, _ , 3—1 H Clmjm_1a]m_1]m(hk 1km km)+0(h ) m=2 S—2 P n —' n . -- . u . —h E . Z Clllalfl H Clme_1aJm—1Jm [1,....ln=1 jl....,jn=1 m=2 S—2 P , . . . Sk—1~. _ . 5—1 2 Cln+13nh . Z GJWJn+1M+125KIOOP= Z rj, J=0 then Kloop S—Sk__1> Z T'jZTk, j=k SOS—Sk_127‘k+l. This proves Proposition 2.3.21. Proposition 2.3.22. Suppose we solve the error equation (2.78) using algorithm 68 { 2. 83). Then the numerical error vector is elk] —o‘[k]+I5k—1+15(k‘1)(im,§[kl)+... m+1— m hSk-1+Tk ~ k—l ~ k -. ~ k—l ~ 1.- + —T'——G7(.k_1)(tm.§2lnl) + hbk-l R£k+1)(tm. olnl) k. t _ / m+1€(k—1)(T)d7+ 0(hM+1), (2.126) tm ~[kl _ 1 where 9m — Wm mid" fom (T)dT) and M—Sk_1 M- U Rrk+1 (1” [Q )— — Z hq q=rk+1 . . w w w with C(k— 1-) (t, 9) ~ 0,1(1)for qll +qb g M — Sk—l' Here (111 “7212 17;. are tqt'Qchj qut qy qt 2 constant coefficients determined by (2.83). Remark 2.3.23. In implementing algorithm ( 2. 83), one needs to evaluate f (t, n(t)) at some intermediate stage, e. g., att = to + cjh. This results in additional function evaluations of f (t,n(t)), which is usually the most expensive part of an algorithm. In our implementation, we avoid this extra cost by performing a polynomial inter- polation. Specifically, given f = (f(t0,770), . . .,f(tm,17m), . . .,f(tM,77M)), we con- struct an M m degree Lagrange interpolant, LM (t, f) and approzcimate f (t, 77(t))| t=t 0 + Cj h using LM(t = to + cjh, f) up to O(hM+1). Remark 2.3.24. T o evaluate the integral term in the right hand side of algorithm k— 1) ( 2 83), we construct a Lagrange interpolant LA! (t c( ) and use 69 fifqflflhm24Nmaa=12.m. to approximate the integral term filo-H: h ((k_1l(t)dt, i = 1, 2,. . . ,p, with an error 0(hM+2). Proof. (of Lemma 2.3.17). The proof is identical to that given in Section 5.2.3 of [18], but we repeat it here for completeness. 1. First we prove part 1. of the lemma. By subtracting the numerical error vector, (2.126), from the integrated error equation, (2.122), we obtain 111 2 1k—11_,1k—11 em+1 em+l m+1 M+1—Sk_1ht “21+th 2 3.70”“ 1)(tm. if. 11) i=Tk+1 hsk— Izhfa “391- _1) elli’l]+E(k‘1)(tm) my hSk—l k 1 _Ga—11,Mlnl+Eklhm) _1 ’ma hsk— 1 ,5-[k 1] (k— 1) s _ ~ +E( (t m) lid-F2 _hk 13,76“ tm, 6’” hsk-l +C’)(h ) = elg] + 21%)] + rgkgl] — rgcln + 0(hM+2), (2.127) 70 i=1 2! dti_1 _512— 1) (1,, 51,1; 1] +E< (’1 U1 m(>)) rk khi (Ii—1 = a jam (f(tm,y(tmll — f (tm,fl(k)(tm))) _zrzklhTZh i— 1L1 . dtil S- 2 ' ' _1]+18] ( (11)" ay£(tm,y(tm))(e (em )1 +O(( elk] )5 )1)) , i=1 71 and since 3k. = Sk—l + rk, k 1 I ~ ~ 1 TlJ—n— ]_ h8k_1 Z Halli—10m» in ]) (2.129) i=Tk+1 [VI—8k . 7‘ ' k~ k— 1 : Sk+l zdz d ~[ ]) h 2:0 (rk+i+1)! dti ((#77661771an k- ll (—k 1) k-1 ~ 6[ +E( (tm T[2,m]=hbl°—1Rrk+l tm. 7” h3k_ 1 (2.130) twirl—8k : hsk+1 2 hi '_ "'t'q 7'_0 qt+qu, 13.8) where y(t ) is the exact solution and n(k‘1)(t) is an filth degree polyno- mial interpolating nlk ‘1]. Note that the error function e(k—1)(t) is not a polynomial in general. 2. Denote the residual function as WW) 2 1n>’1> 1511 n“ 1>11» — 1N1t.n1t>), 81 and compute the integral of the residual. For example, /tm+lc(k_1)(r)dr (3.9) t'0 M g 77%;? — yo ’ 27mfif50jfl7E-k—1A) + fNUjmg-ATIAD, J: where 7m. 7- are the coefficients that result from approximating the in- tegral by interpolatory quadrature formulas, as in [25] (or as described in Section 2.2, w, ”[k—l] = 7nd- equals the quantity given in equation (2.32)), and nAk—1)(tj)- . Compute the numerical error vector, denoted by .. - k k k 5W = (0]) A, . . . ,5[,,], . . . 3%)). (3.10) which is an rkth order approximation to the error el .1-111-11.....i-u,...,ei-”>, where egg—1A = e(k_1)(tm) is the value of the exact error function (3.8) at tn} . To compute 51A] by an rkth order ARK method, we first discretize the 82 integral form of the error equation, 1129—11111) = 1511.,11(k-1)11>+ elk—”11» — 151111-111» + 1N11.n11)+ (AA—”1m — 1N11.wz(k—1)11>> = 1511. n1t) + 1211-111) — E+ 1211-111) - E>- 1N1t.n> 2 F511. elk—”11> — 511—1)(,,)+ FN1t.Q(A’A)(t) — Bit—”10) e 0511. 1294111)) + GN1t.Q(k—1>1t>), Q(k71)(0) = 0. (3.11) where (WC—Uh) = (AA-9(1) + EAk“1)(t), (3.12) ra1t.e)— 10.11. ilk-1111)). (1 = N, S, The fa, a = N, S notation above (with Q and the approximation for E) is used in the numerical computation, while the Ga, 11 = N, S notation is employed in the theoretical explanations and proof in Section 3.4. We apply a pk-stage rkth order ARK method to (3.11) and denote the rkth order numerical approximation to (Jig—1A = QUE-.1) (tin) by Diff] 83 We therefore have, for a- = N, S, i = 1,2, . . . ,pk and m = 0,1,...,M—1: ki =Gaz(tm+c-,hf2]nA(A+h (Aid/V aij kj N+:::aiSJkS) j j=1 N N+ as S E[k— 1 h-(JZCL aI-J-kj+ :1azjkj)— ](t7n+th)), [kl [kl N N 5-5 am +1 _ —o,,, +h 21b, 1,, +1), 1,,- ). (3.14) i=1 From equations (3.12), we compute our numerical error vector (3.10), .5111 z 9111 _ Eve-11, where the components of Elk—ll are Egg—1A = E (A: ‘1)(lm ) In fact, we use 6AA] [k] tm+1 k—l +1— _Qm+1 — to c( )(r) dr (3.9) [k] z Q7n+1 M — (71%;? - :1/0 - Z 7m,j(f5(tjan]A—1A)+fjv(’jUAA 1A»). j=0 where we have approximated the integral by interpolatory quadrature formulas, as in [25]. In computing (3.14), one needs to evaluate fa(t,n(t)), a = N, S, at some intermediate stage, e.g., at t = to + cj h. This results in additional function evaluations of fa(t,r)(t)), which is usually the most expensive part of an algorithm. In our implementation, we avoid this extra cost by 84 performing a polynomial interpolation. Specifically, given fa : (f(.k(l'017l0)1- - - a f(1(t'nli 7’7”), ' ' '1f(X(LA/I:’7A/I))1 we construct an M ”A degree Lagrange interpolant, LM(t, fa) and ap- proximate fa(t177(t))|t=t0+cjh using LM(t = t0+cjh, fa) up to O(hM+1). Again note that fa, Fa are more useful for understanding the numerical implementation, while Ga is employed in the theory in Section 3.4. 4. Update the numerical solution nlkl = n[k_ll + 6m. 3.3.2 Example with Second Order ARK Now we present an example where a second order ARK scheme is employed in the IDC framework. We consider the IMEX RK2 scheme from [3], and illustrate that our formulation is equivalent to Layton’s formulation [47] for this specific IMEX RK2 scheme. Writing the ARK scheme (we denote it by ARK2ARS) in a Butcher table, we have: 101—ggd1—d Ol—ggOI—gg whereg= 1—-‘-é—§ andd= —2—3‘/—§. 85 The prediction loop solving (3.1), for a = N, S and m = 0,1,. . . , M — 1, is: . 0 k? = f(}'(l'"l.77l7An])1 , 0 ,, A? = fora"). + 9”»: Win] + AA9(kiAV + kgllv 0 15,} ___ fa(t,m+1,n[,,l + 1.1.119 + (1 — (1)14" + (1 — g)k28 + 925?». 27.],0,A+1= ADA +h((1 —g)(k{V+k§) +9125;V +23%, The correction loop solving (3.11) is 1.0 = 2011..., 1215;] — Elk—1112...» kg = Fa(tm+1,a[fil + h(dk{V + (1 — (1)2? +11 — g)k§ + 91239) _ Elk 11(f (Am+1)) “[5121 = 1215-”? + 2111 — (1)129 + 259 > + 9129’ + 2:»? >1. Aifih *9i:1+1 ‘ EAAA—IAAA'erll' Letting 1)AA A —(5[:)A+ nAA— 1A an nd QAAA: 6AA] + EAA 1A ,and noting Q%+C(fs(t,n(k"1)(t))+ fN(t,n(k—A)(t))) from [47] is an approximation of fttrrranrCh C(k—1)(T)dT, we see that the formulation of IDC constructed with ARK2ARS in [47] is equivalent to ours, namely: «bfilgwlfilih 9AAiA+2g12N+2§)— BAA—lAltm‘l'ghl, 2231., =n[§;iA+12.n AAA+ 2121’V +11— «112 25V +11—g>2§+g2§>—Elk'”1tm+1>, k k-l "inAH =777An+1A +‘AinA A +1 The notation on the left is as in [47], while the notation on the right is ours. Although 86 the two formulations are equivalent for this type of second order ARK method, the incorporation of arbitrary order ARK methods in the corrections is unclear in [47] but is obvious in equation (3.14). 3.4 Theoretical Analysis of IDC-ARK Methods In this section the analysis for semi-implicit IDC—ARK methods is presented. Recall some mathematical preliminaries related to the smoothness of discrete data sets were introduced in Section 2.3.2. The local truncation error estimates of semi-implicit IDC-ARK are provided, along with a corollary for the local truncation error of im- plicit IDC-RK integrators, as an alternative to the presentation in Section 2.3.2. This paper closely follows the formulation presented in [18], as well as Section 2.3.2. It is recommended that the interested reader who wishes to understand the analysis in this work first reads Section 2.3.2, which presents some mathematical prelimi- naries regarding the smoothness of discrete and continuous data sets,and the basic framework for analyzing the local truncation error. The following theorem states that, with the assumption that the IVP is smooth enough, using an rth order ARK method in the correction loops results in an r-order increase in accuracy of the IDC method, but the overall IDC order of accuracy will not be higher than M + 1, the number of nodes in the subinterval [tn, tn+1l- E.g., if an rth order ARK method is used in each IDC loop, then the IDC method has order rnin(r =2 (1 prediction + K1001) corrections), M + 1). Theorem 3.4.1. Let y(t), the solution to the [VP (3.1), have at least a (and f5 and f N have at least a — 1) degrees of smoothness in the continuous sense, where 87 o 2 M +2. Consider one time interval of an IDC method with t E [0, H] and M +1 uniformly distributed quadrature points ( 3. 7 ) Suppose an r0 th order ARK method )th order ARK methods loop are used in K loop correction steps. For simplicity, assume that cj 9— cN = c5 .7 J the number of stages for the implicit and explicit parts of each ARK method is the (3.5) is used in the prediction step and (r1,r2,. . . ,rK and same. Let 3k = 2320 rj. If SKloop S M + 1, then the local truncation error is of 3 +1 order (9(h K1001? ). The proof of Theorem 3.4.1 follows by induction from Lemmas 3.4.3 and 3.4.4, below, for the prediction and correction steps, respectively, which discuss the 10- cal truncation error and smoothness of the rescaled error when general high order ARK schemes are applied in the prediction and correction loops of IDC methods. Lemma 3.4.3 is the first case, and Lemma 3.4.4 is the induction step. Note that the smoothness results of the rescaled error vectors in both lemmas are essential for the proof of the correction loop in Lemma 3.4.4. The smoothness proofs of Lemmas 3.4.3 and 3.4.4 (analogous to those in [18] and in Section 2.3.2) are quite technical. Since they are presented earlier for the implicit IDC-RK methods, here we simply state and sketch a proof for the smoothness of the rescaled error vector following a forward-backward Euler (FEBE) prediction step in Lemma 3.4.2. The smoothness proof for the more general case is similar in spirit. Once the smooth- ness of the rescaled error vector is proved, the magnitude of the error vector follows directly from Definition 3.2.2. Lemma 3.4.2. ( F EBE prediction step): Consider an IDC method constructed using M + 1 uniformly distributed nodes, and a FEBE integrator for the prediction step. Let y(t), the solution to the I VP (3.1), have at least a Z M +2 degrees of smoothness, and let filo] = (r;([)0], . . . , gig], . . . ,ngefl), be the numerical solution computed after the prediction step. Then the rescaled error vector, 30] = )12e'lol, has M degrees of 88 smoothness in the discrete sense. Proof. We drop the superscript [0] as there is no ambiguity. Since 7lm+1 = 77TH. + th (tm, 77m) + hfS(tm+1a Urn—Fl): (315) and 1’~m.+1 . tm+l 31m+1=ym+/ /N(T,;U(T))dT+/t fs(r,:l/(T))dr. m m 0—2 = ym + thUm, ym) + Z hi+1 JlfN . t . 0+1)! all" (mflm) _2 . - U h2+1 dlfs 1 (H1)! dt'l + hfSUm+1a Urn-+1) + (tm+1a ?/7n+1) + 0070) (316) 2': Subtracting equation (3.15) from equation (3.16) gives 5m+1 = 6m + h/(fNUm, yin) ‘ fN(tma 77ml) + h(fS(tm+1wym+1) — fS(tm+17 77m+1)) 0-2 hi+1 dif/v 0-2 hi+1 difs . (i+1)! dt'i , (1+1)! dii i=1 i=1 (tm. ym) + (tm+1, ym+1) + 0(h0). where em+1 = ym+1 — nm+1 is the error at tm+1. Let um = (fNUm, 31m) - fNUm, 77ml) + (f3(tm+1a ym+1) — fS(tnz+la77m+1))v and 0—2 “‘+1 di 1 fN i=1 (r+1)!( or; am an, gm) + fi>. 7'm = 89 To prove the smoothness of the rescaled error vector, e = e/ h, we will use an inductive approach with respect to s, the degree of smoothness. First, note that a discrete differentiation of the rescaled error vector gives, . é .. 1 — gm —2 d: =—T—ni—————=— — Oh" 3.17 ( 1()7n' h h +T h m+ ( )7 ( ) We are now ready to prove that g has M degrees of smoothness by induction. Since llglloo ~ 0(h), thus, 5 has at least zero degrees of smoothness in the discrete sense. Assume that g? has 3 g M — 1 degrees of smoothness. We will show that (Tl—é has 3 degrees of smoothness, from which we can conclude that g has (3 + 1) degrees of smoothness. First, 0—2 .- 1' um : 2: 562m 6" -)(tm,ym )+ 0:12162'lem-l-1—8y2' (tvm+1 ym—H) i=1 ° y i=1 + 0((6m)0’1)+ (9((em2+1)0_1) a— —2 ' 1~ (9sz( 1 = Z— 2! elnhz 8y z“ (tm gm) +312 7! ém+1hiaz 8y 2' fium+1 ym+ll i=1 2:1 + 0((hé7n,)0--.1) + 0((hém+l)0_l) where we have performed a Taylor expansion of f N(t, 77m) about 3; = gm and of fS(t,77m+1) about 3) = ym+1. Denote fyi to represent either 82f or —;—Sf. 8y _i’ Since fyi has (a — i — 1) degrees of smoothness in the continuous sense, fyi = [fyz'(t0, yo), . . . , fyz-(tM, yM)] has (a — i — 1) degrees of smoothness in the discrete sense. Consequently, hi’lfjfi has (a — 2) degrees of smoothness, which implies that 97113 has min (a — 2, 5) degrees of smoothness. Also 7}; has at least 3 degrees of smoothness from the smoothness property of the IVP (3.1). Hence dig has 3 degrees 90 of smoothness => g has (5 + 1) degrees of smoothness. Since this argument holds for a 2 M + 2, we can conclude that g has M degrees of smoothness. C] Now we present the lemmas required for the proof of Theorem 3.4.1. Lemma 3.4.3. (Prediction step): Let y(t), the solution to the [VP (3.1), have at least a (and f S and f N have at least a — 1) degrees of smoothness in the con- tinuous sense, where a 2 M + 2. Consider one time interval of an IDC method with t E [0, H] and uniformly distributed quadrature points (3.7). Let nlol = (7750] , ngolr . . , 1).,[2] ,. . . , 17%?) be the numerical solution on the quadrature points (3.7), obtained using an rOth order ARK method { as described in Theorem 3.4.1) at the prediction step. Then: I. The error vector em] = y — 7)[0] satisfies Ilelol “00 ~ O(hr0+1). 2. The rescaled error vector elOl = fib—elol has min(o—r0, M) degrees of smooth- ness in the discrete sense. Lemma 3.4.4. (Correction step): Let y(t), the solution to the [VP (3.1), have at least a (and f S and fN have at least a - 1) degrees of smoothness in the continuous sense, where a 2 M+2. Consider one time interval of an IDC method witht E [0, H] and uniformly distributed quadrature points (3.7). Suppose e[k—1] ~ 0(hsk‘1+l) and the rescaled error vector Elk‘ll = fizevc—ll has M + 1 — sk_1 degrees of smoothness in the discrete sense. Then, after the kth correction step, provided an rk th order ARK method (as described in Theorem 3.4.1) is used and k g Kloop’ 1. The error vector elk] satisfies Hewlloo N O(h8k+1). 2. The rescaled error vector elk] = Lelk] has M+1— 3k degrees of smoothness h5k in the discrete sense. Note that most results and the proofs (e.g. smoothness) are similar to the analysis in Section 2.3.2 (implicit IDC-RK) and in [18] (explicit IDC-RK). The main difference between the proof of the truncation error for IDC-RK (explicit or implicit) and IDC-ARK methods lies in a portion of the proof of Lemma 3.4.4, stated below as Proposition 3.4.5, following some notational details. The proposition asserts that the rescaled error. vectors (rescaled by hSk—l), obtained from formally applying ARK methods to a rescaled error equation, are equivalent to the error resulting from the numerical formulation that we implement. Thus, since an rkth order ARK method applied to the rescaled error equations achieves rkth order, the unscaled numerical implementation is 0(h8k‘1+rk), provided the exact solution y(t) and the functions f N, f S are sufficiently smooth. For purposes of clarity, before stating the proposition, we first introduce some notational details related to the rescaling, and recall the error equation (3.11). The analysis for the correction steps in this section will rely on this form of the error equation. Assume that after the (k — 1)th correction, the method is 0(hSk-1). we rescale the error 6(k—1) (3.8), and (QM—1) and Ca from the error equation (3.11), by hSk-l. Then the rescaled equations (3.18) and (3.19) are (9(1) (see Section 2.3.2, or [18]). Define the rescaled quantities as ak—lla) .-= flew—1hr). (3.18a) @(k“1)(t) = fiQlk‘lla), (3.18b) off—Ila, (Elk‘llan = 11313—1 can, hSk—lo(k—1l(t)), a = N, s. (3.18c) 92 Then the rescaled error equation is given by k—l) (Few—Wm = El, H) (t,55(k—1)(t))+ 6f,» (“ya—mm, (3.19) @(k‘llu) = 0. A p—stage, rth order ARK method (3.5) applied to (3.19) gives, for i = 1,2, . . . ,p anda=N,S: ”a_~(k‘1)aS~S)) k1. _Ga (t0+c thok— 11+1; (ZaN a”. k]. azjkj) ~[k1 ~[k—11 7” N~N s~s ol =Q0 with), k, +1), k,- ), (3.20) i=1 where the Greek letter fl denotes the rescaled solution in the numerical space, and égc— 1] equations (3.14) with FN, PS. We use (3.14) with G'N, GS instead of FN, FS for = @(k—1)(t0), In the actual implementation, we discretize (3.11) as in the analysis. The following proposition is where the essential difference between the proof of the truncation error for IDC-RK [18] and IDC-ARK methods lies. Proposition 3.4.5. The Taylor series for the rescaled exact error, Elk—1)(t0+h) = T’s—_TBUC—Uito + h), and for the rescaled numerical error, $6)“ in (3.14), h coincide up to and including the term hr for a sufiiciently smooth error function EVE—1) t , and the exact solution y has a and f , f have a — 1 degrees of N S smoothness. Proof. It suffices to prove the following claim: risk—1'15? = k? + coma—1), Vi=1,2,...,p, a = N, s, (3.21) 93 where k? and ii? are as in equations (3.14) and (3.20). Set 3—1 A.i=h(z a3’h‘k- 1EN+ +2111 c”risk— 11:5) J=1 p N :3 ~N S 8 S = k—l .. k—l hZaUh k] +a )1 k1) (322) F1 F1 N N p N N+ s s B,_h(Zla3kJ +2333): 11:10; aijkj+ (L319) (323) J= = j1= (10:18! (.1, + I ,(o 152/1’ ”BU—1, llj— l—!ayl1 0 (3.] b 0k— where a = N, S, and the equalities in (3.22) and (3.23) hold since az’g = 0 forj 2 i. To prove (3.21), we prove that hSk—IZZQ —- k? = 0(hn) for n == 1, 2, . . . ,a — 1, for i = 1, 2, . . . , p, and Oz = N, S using induction with respect to n; i.e., we show that 9 a a: n aN(h9k— 1N N “3 913—119- S hk— 1"i ki h ZZwlgn kj kJn)+filj,n(h kJn kJn» l,nj.n + who-1), (3.24) where the notation is as follows. 0—20—2 0—2 2= 22 Z Z= 22: Z In 211112 In in J1 J2 J'n 271—711 fllj,Nn22245112210“fllj,n=:: fiCa, where for the )8 coefficients, C represents some C8’ or CIS. , a represents some all; aS or azj, and ON and OS denote that the C, a coefficients that appear in the sums 94 depend on a and N or a and S, respectively. For example, if n = 2, then (IN. 101 N “1N .N ya as S (IN filj,=2 Cllia iJ1Cl231 “J112+Cllza z3161231611132 “S: 10: N N US a a3 S (1.3 ”112 ”112' ”iJ1CI2J1aJ1J'2 +011,a,,-1012,-1a,-1j2 which is clarified in the expansions below. We prove (3.24), and hence (3.21), using induction with respect to n, the power of h in equation (3.24). o n = 1: Using the rescaled quantities (3.18) and Taylor expanding Ga(t,y) about Q(k—1’ in y, we obtain for a = N, S: 0 risk—17.7? — k? (3.25) 5,, 1 (k 1) (k 1) p N~N s~s h C (t0+chQ ( )+h§_:(a,]k] M3113» J'=1 — Ga(t0 + Cih Q“ 1)(tO) + N 201319;] + ais’jkgé’)) j=1 s ~(k 1) p N~N s~s = 0050.0 + C’lh” h k—1(Q (to) + ’1 2:1(02'3'193' + 092'ij ll) J: — Ga(t0 + 3,11, (209—”(10) + 3,) = Ga(t0 + C,,-h, Q(k*1)(l0) + A?) — Ga(t0 + Cih,Q8k—1)+ Bi) 0 2 1 8’1 1 1 — 1: 111' —l1G3y 0(10 + C’lh’k—Q(1)(t0)l((Ail 1 - (Bi) 1) + 00104). (3.26) 95 ll—v l 1- . Factoring (Ai )11 — (Bl-)ll = (Ai — Bi) Zvl— 1 Ai I}; 1, we obtam 8 ~__ in ,(r h k 11.2. —1,. 0—2 —1 = Z cfi,(A,—B,)+ovla ) l1=12 _ 3k—1N- N (,5 5k—175_S h 12:0119 2(0 (13101 kn kJ1’+aiJ1(h' le kn» 71: 1 +C’)1(h:_1) p N N aS 3k 1~S S _. 5.—1 _ 0’ ‘.— . _ . —h :2 Z(Cf1',a,j1( k1 k )+ 0,173 3.101 kn 1:31)) l1=1j1=1 +O(h”—1) _ 1 ON 8k_1TN_ ,N 0‘3 3k—1~S_ S _h 2293,3101 1,1 1.] )+13, 1(1 k] —jk1)) l,1 J,1 +O(}1"_1), Vi=1,2,...,p. (3.27) Note here that SON =CO‘ aN and 13 03— Caa S lj,1 lli iJl’ 'lj,1_llialj1 o n + 1: Assume (3.24) holds for n < a — 1 (for both a = N and S). Following the same process of Taylor expanding as in equations (3.26) and (3.27) in the 96 n = 1 case, we see that for or = N, S: hSk—lha—kq hn (mic—171V- N “S ,—9k 1S .5 Z};Z(filc;{yn( kJn kJn)+/jlj,n(h kin km» ,,an +O(h"‘1) 27) n ON 0—2 p = h ZZWIMUI Z ' Z ’1” 3,11 ln+1=1Jn+1=1 N N 8 N N C , .. llk— 1k —k- ( ln+ljnaJ7Un+1( Jri+l .7n+l) N S 3k— 1 S S 0—1 . . . h k —k- 0 h C’n+117?a-77192n+1( jn+1 jn+1))+ ( )) +filJ,n(h :2 :( l'=n+1 ljn+1;1 S S 5 ~S S 0.__1 +C h k— 1k. .16. h 11.4.1111 11.1.1.3 1,1, ,n,,» + 0< 1) + 0(119‘1) =hn+l Z Z hsk— 1k’-V— kN qrsi+1Jnaj71jn+1( 31 31) l,n+1Ln+1 ((fllajjl’l 111i+1JnaiY1Jn+1 11,310an ’gi-tljnalbjnHXhSk lkli kN’ +< ZNCIIXHJH afnjn+1 +fil1fllcgi+lm aJSanH’ (hsk 119.5; H 4159 +1))+O(h,U—1) —hn+1lnz—I;1in:+—1 131] n(+1 (1,91: lkJ’le-kil’m) +r3”.5+1(h3k— ”5311.31-23“ 1))+O(h,0—1), Vi=1,2,...,p. 97 where aJ’V CIA] (IN 0‘3 ‘S N . Blj n+1 :filelJ/LC [n+1Jn aJan+1 +51]. n61n+1JnaJan+F 'SIJJL+1: IJ n Cln+1jan Jan+1+lJnCln+1jnaJan+1 It follows that (3.21) holds. The rest of the proof of Proposition 3.4.5 is similar to the implicit IDC-RK analysis in Section 2.3.2. Cl Now we state the implicit IDC-RK results again as a corollary to the semi-implicit theorem. Corollary 3.4.6. Let y(t) be the solution to the IVP y’(t) = f(t,y), £60171], y(0) = yo. and y(t) has at least a {o 2 M+2) and f has at least 0—1 degrees of smoothness in the continuous sense. Consider one time interval of an IDC method with t E [0, [I] and uniformly distributed quadrature points (3.7). Suppose an rOth order implicit RK method is used in the prediction step and (r1,r2,. . . ,rKIOOPJth order implicit RK methods are used in K1002) correction steps. Let 3k = 2;:0 rj. If SKloop 3 +1 M + 1, then the local truncation error is of order (9(h loop ). Proof. The result follows from Theorem 3.4.1, if we set f N = 0. E] 3.5 Stability The stability of an explicit (or an implicit) scheme can be measured by the conven- tionally accepted definition of a stability region [37], as presented in Section 2.2.3, 98 Definition 2.2.2. For a semi-implicit numerical method, stability is less easily defined. In the literature, various measures for stability have been proposed for semi-implicit [60, 9] or additive Runge—Kutta methods [62, 56, 11]. In fact, all of these definitions can be applied to any semi-implicit scheme. We focus on two of these definitions. Minion suggests splitting the Dahlquist-type equation in such a way that the stiff part is represented as the real part of the eigenvalue and the nonstiff part is represented by the imaginary part of the eigenvalue, which is useful for seeing whether the numerical method works for a convection-diffusion problem, but may not be reasonable when considering a problem with stiff oscillations. Liu and Zou define an A(a)-stable semi-implicit method (similar to [62], except in [62], they also present a boundary locus-type method), which could cover a wider class of problems [56]. 3.5.1 Stability for Convection-Diffusion Type Problem First consider the definition from [60]. Definition 3.5.1. The amplification factor for a semi—implicit numerical method, Am.(/\), with A = a + ifi, is interpreted as the numerical solution to y'(t) = OW) + 213W). 11(0) = 1, (328) where the explicit component of the numerical method is applied to the imaginary part ifly(t), and the implicit component of the numerical method is applied to the real part (xi/(t), after a time step of size 1 for a E IR and fl 6 IR, i.e., Am()\) = y(l). Definition 3.5.2. The stability region, S, for a semi—implicit numerical method, is the subset of the complex plane C, consisting of all A such that Am()\) S 1 for ODE 3.28, i.e., S = {A : Am()\) 3 1}. 99 In Figure 3.1a, the stability regions for 3rd, 6th, 9th, and 12th order IDC- FEBE (IDC constructed using forward and backward Euler) with 2, 5, 8, and 11 correction 100ps, respectively, are computed numerically and plotted. Figure 3.1b shows 3rd, 6th, 9th, and 12th order IDC-ARK3KC (IDC constructed with 3rd order ARK3KC) with 0, 1, 2, and 3 correction loops, respectively. The shaded regions represent the stable regions. Note that for the same order method, IDC-ARK3KC shows a marked improvement in stability over IDC-FEBE. This correlates with the results in [18], where it was shown that IDC-RK methods were found to possess better stability properties than SDC constructed with forward Euler integrators [18]. Such a stability measure provides insight on the stability property of a scaled IDC-FEBE ’ scaled IDC-ARK3KC 0.6* 0 6 g 0.4% = _ g 0.4» E e "' x 3rd ' " x 3rd 0.2- 0 6th 0.2» 0 6th 0 9th 0 9th 0 . . A 12th 1 O L A 12th -0.05 0 0.05 0.1 -0.05 0 0.05 0.1 Rem R90») (a) IDC-FEBE (b) IDC-ARK3KC Figure 3.1: Stability regions for 3rd, 6th, 9th, and 12th order IDC constructed using (a) forward and backward Euler with 2, 5, 8, and 11 correction loops (and 3, 6, 9, 12 uniform quadrature nodes), respectively, and (b) 3rd order ARK3KC with 0, 1, 2, and 3 correction loops (and 3, 6, 9, and 12 uniform quadrature nodes), respectively. The crosses denote the 3rd order stability regions, diamonds for 6th order, circles for 9th order, and triangles for 12th order. The stability regions are scaled by the number of implicit function evaluations ((#stages)(M + 1)M, where FEBE has 1 stage, ARK3KC has 4 stages). semi—implicit method (such as an IDC-ARK method), especially when it is applied to a convection-diffusion problem (e.g., the problem in Section 3.6.2), where the 100 convection term is treated explicitly and diffusion term implicitly. 3.5.2 Stability for Other Types of Problems We remark that the stability property of the semi—implicit scheme greatly depends on the splitting of equation (3.1), therefore, the stability analysis must be adapted for specific problems. Hence we also present a definition from [56] that may apply to a broader class of problems, although our results are inconclusive. Rather than the Dahlquist equation (2.45), we consider a test problem of the form y'(t) = A5.?! + ANy. y(0) = 1- (3.29) /\ S and AN depict the eigenvalues of the Jacobian of f S and f N, respectively, where Re(AS) S O, Re()\N) S 0, and IRe(AS)| >> |Re(AN)|. The semi-implicit amplification factor, denoted Am(AS,AN), is determined by applying the semi- implicit numerical method (e.g., (3.4)) to the test problem (3.29) for one time step, obtaining 3’1 = Am(AS,AN) (3.30) Definition 3.5.3. [56] An A(a)-stability region for a semi—implicit numerical method is m 35:5 W = {AN 6 C s.t. [Am()\S,/\N)[ $1 VAS 6 55,5}, (3.31) 101 where SCfYS = {)‘S E C s.t. [arg(—/\5)[ S 0:, )‘S 75 0} U {0} (3.32) An L—stable semi-implicit method is an A-stable (i.e., A(%)-stable) which also satisfies lim [Am(AS, )‘Nll = 0. (3.33) AS—>0 Stability regions were numerically and analytically determined for some ARK methods to confirm the (A-stable type) stability regions presented in [56]. Definition 3.5.3 with a = 7r/2 was used to calculate the regions. ARK2A1, ARK3A3, and ARK4A2 were chosen due to the obvious difference in size and / or shape among the methods and between their explicit stability regions (i.e., the stability region arising when the explicit part of the ARK method is calculated alone as an RK method) and their semi-implicit stability regions. The general algorithm for computing the semi-implicit stability regions, using (3.29) and (3.30), follows. 1. Grid the left half complex plane for )‘S' 2. For each fixed AS, find the stability region (according to the usual definition, Definition 2.2.2) in terms of )‘N' 3. Take the intersection of the stability regions from step 2 to obtain the ARK stability region. Remark 3.5.4. In step 1, it is necessary to refine the grid for )‘S close to the imaginary axis. These values of /\ S appear to be the troublesome ones that shrink the ARK stability region to smaller than the explicit stability region. It may be unnecessary to include A S with large negative real parts. These observations are supported by [55, 62/ and by the structure of a typical amplification factor’s equation. 102 E.g., in (1 — AS — 9g.) + (1 — ASMN + §A§V AmA ,A = (S N) 1-2/\s+/\% (3.34) the equation for the amplification factor of ARKZAI, we see that large negative real values of ’\S make Am(AS, AN) approximately constant in modulus, while values of A S closer to the imaginary axis may increase the modulus on their own or through interaction with )‘N' Furthermore, [62/ presents a particular case where it is rigor- ously shown that the maximum of the amplification factor occurs on the imaginary axis. At first glance, ARK methods appear to reduce the stability region; however, when one considers that this restriction only applies to the nonstiff eigenvalue, while the stiff eigenvalue may be anywhere in the left-half complex plane, a smaller stability region may not be unsatisfactory [56]. Nevertheless, we are still limited by the lack of I-stability (imaginary axis stability, [36]) for the nonstiff component. The stability regions for IDC-ARK schemes were ascertained numerically as for ARK schemes. A fourth order IDC-ARK method using ARK2A1 has a semi- semi-impl A-stability: IDC4 w/ARK2A1 6"“ ——-..— —— .— s _-fi ~ 7 -e -s -2: x2 0 (a) IDC-ARK2A1 Figure 3.2: Stability region for fourth order IDC with ARK2A1 in prediction and correction loops. A N is plotted, consistent with Definition 3.5.3, and a = 7r /2. 103 implicit stability region (Figure 3.2a) that is much larger than the semi-implicit stability regions for both fourth order ARK methods shown in [56]. Unfortunately, the 8th order IDC-ARK (with ARK2A1) method’s stability region is not A-stable at all (in the semi-implicit sense of Definition 3.5.3); however, the numerical solution for (3.35) and (3.36) computed via the same 8th order IDC-ARK method seems to be stable, perhaps even for c = 0.01. To investigate the loss of (semi-implicit) A-stability, we first fixed AS, and plotted some explicit stability regions only. We achieved very reasonable stability regions that increased with the order of the IDC method, as expected from the results in [18]. Then we fixed AN, and plotted some implicit stability regions only. We also plotted some stability regions for implicit IDC-RK methods constructed with a sec- ond order DIRK method from [2] (Figures 3.3a and 3.3b). Here we discovered the problem. Clearly the stability regions decrease in size as the number of correction loops increases and frequently the higher order methods lose A-stability. In partic- ular, we notice that in Figure 3.3b, the 10th order method loses A-stability; even A(a)—stability is lost. Hence, although the explicit region increases, as seen in [18], the implicit region loses A(a)-stability as number of IDC corrections increases, and the semi-implicit method struggles between improving and worsening its stability. Thus far the results of measuring semi-implicit stability via Definition 3.5.3 do not look clear, although the main difficulty appears to arise in the implicit part of the method, where increasing the order of an implicit method decreases its stability re- gion. This decrease is well-known for BDF methods, but generally is unpredictable for implicit RK methods. Also, [56]’s method given in Definition 3.5.3 appears to be overly restrictive since some numerical solutions still appeared to be stable despite the numerical method’s loss of semi-implicit stability. Therefore, we chose not to pursue this means of measuring semi-implicit stability further. 104 Stabilityi IDCBDIRK2p106, 0. 1. 2, 3 corrections Stability: |DC10DIRK2p106 0, 1, 2, 3, 4 corrections 400 -200 o 200 (a) IDC-DIRK2, 8 nodes (b) IDC-DIRK2, 10 nodes Figure 3.3: Stability regions for IDC with 2nd order DIRK (from [2]) in prediction and correction loops. The labels are the number of correction loops, and are placed outside of the stability regions. From 0, 1, . . . ,4 correction loops, the IDC methods have corresponding order of accuracy 2,4, . . . , 10. Stability regions are standard, according to Definition 2.2.2. 3.6 Numerical Results In this section we describe some numerical experiments with various IVPs and an advection diflusion equation. Numerical experiments support Theorem 3.4.1, which states that the semi-implicit IDC-ARK methods have similar order of accuracy rules as for IDC with explicit RK methods [18]. Additionally, in most cases, higher order IDC-ARK presents an efficiency advantage over existing IMEX ARK schemes. Sev- eral different ARK methods of various orders were tested Within the IDC prediction and correction loops: second order ARK2A1 [56] and ARK2ARS [3], third order ARK3KC [45] and ARK3BHR1 (Appendix 1. in [8]), and fourth order ARK4A2 [56] and ARK4KC [45]. For efficiency, these IDC-ARK methods were compared with their constituent ARK methods. This section first lists three ODE test problems, followed by the combined results for order of accuracy, comments on order reduction, and efficiency results. Then an advection-diffusion equation is described, followed by its order of accuracy results and comments on an improved CFL condition. 105 We further comment that the focus of our study is not the choice of f N, f S, and thus we make naive choices of f N and f S in each example below, acknowledging that better choices are likely to exist (in particular, there is much study on the choice of f N, f S for PDES, such as in [43, 44, 34]). Note also, that, in most of the literature, the examples of semi-implicit methods used to solve problems with stiff regions only test convergence and efficiency up to the time just before the stiff layer, at which point they claim that adaptive steps should be taken [45, 8, 3]. In this work, we have tested through times that extend beyond the stiff layer, and hence our results may appear to be worse than they should. Testing our methods in the standard way is likely to produce much better results than shown here. 3.6.1 ODE Test Problems Test Problems The following test problems, an initial layer problem and Van der Pol’s oscillator, were used in numerical determination of order of accuracy and efficiency. Each IVP was first computed in a completely nonstiff situation (6 = 1). Then varying levels of stiffness, which caused each IVP to have a stiff and nonstiff component, were considered. 1. Initial layer problem Pareschi and Russo’s example [45], with nonequilibrium, or perturbed, initial conditions y1(0) = 7r/2 and y2(0) = 1/2 is tested: 31’1“) = W20), t E [0, 4]. (3.353.) nan) = 1110) + 1(sins/1(0) — 312(0), (3.35m 106 where we choose fN = 3 f5 = 1 . y1 2(Slny1 — 112) . Van der Pol’s equation A standard nonlinear oscillatory test problem, Van der Pol’s equation [60], is examined: yi(l) = 1120). t e [0, 4], (3.36a) its) 1<—y1+ (1 - ween/2(0), (33%) E with initial conditions y1(0) = 2 and y2(0) = 2/3, and we choose Order of Accuracy In the following convergence study, error is plotted versus the stepsize H, where the error at the final time is calculated by comparing successive solutions. For the IVPs, the error is the absolute error at the final time. In agreement with Theorem 3.4.1, the order of accuracy of the IDC method increases by the order of the ARK method used for each prediction and correction loop, as seen for the initial layer problem in Figure 3.4a and for Van der Pol’s oscillator in Figures 3.4b and 3.5a. Here we choose 6. = 1, pertaining to the nonstiff scenario. The circles in Figures 3.4a and 3.4b correspond to the error of the 3rd order ARK3KC method (equivalent to one prediction loop of IDC). It has slope z 3, as 107 expected, since a third order ARK method is used and there are M +1 2 3 nodes in each subinterval [tn, tn+1]. After one correction loop computed using a third order ARK method (and M + 1 = 6), the crosses show sixth order convergence — an increase in the order by three. Similarly, a second correction increases the order to 9, as seen with the triangles. Thus three loops (1 prediction + 2 corrections) of a 3rd order ARK method produces a 9th order IDC method, when there are M + 1 = 9 nodes in each subinterval [tn, tn+1l- Similarly, in Figure 3.5a, we see a fourth order increase with each loop of the IDC method when the fourth order ARK4KC method is used (giving an eighth order method after two loops). We also solved these test problems with various orders of ARK methods (2, 3, 4) and IDC-ARK (4, 6, 8, 9) constructed using second order ARK2A1 [56] and ARK2ARS [3], third order ARK3KC [45] and ARKBBHRI (Appendix 1. in [8]), and fourth order ARK4A2 [56] and ARK4KC [45] (not all shown). In each implementation, we anticipate that # loops * ARK order = expected IDC order. In general, the expected order of accuracy is seen for nonstiff problems. Order Reduction For stiff problems, a phenomenon known as order reduction is observed. The ef— fects of increasing the stiffness of the problems (i.e. decreasing c) on the order of convergence can be seen in Figures 3.6a and 3.6b. The plots show the error for computing the initial layer problem, IVP 3.35 and Van der Pol’s oscillator, (3.36) in a stiff case (6 = 10-3) using third order ARK3KC in O, 1, and 2 correction loops of IDC. Order reduction that is c-dependent is apparent for both the initial layer problem, IVP (3.35), and the oscillator, IVP (3.36). Once the stepsize is below 5, however, it appears that the expected IDC order is achieved. For the initial layer IVP (3.35), the order reduction clearly still allows for convergence of the method within the order reduction regime (Figure 3.6a). Van der Pol’s oscillator exhibits 108 1O +|DCB, y1 >3“ '9'", y2 ;_ -5 -—|DCG, y1 .>.~ 10 y2 g +|o<329, y1 h _V_", y g 10-10. (D -5 a 10 H (a) IVP 3.35, e = 1 (nonstiff) 10° +|DC3, y1 >‘~\' "'9'", y2 ;_ _5 --|DC6, y1 f 10 ~..-».. y2 3 +|DC9, y1 a; _v_"’ y2 £3 10"“)- (U 10'5 (b) VdP (3.36), e = 1 (nonstiff) Figure 3.4: Convergence study of absolute error at T = 4 vs H, for 3rd, 6th, and 9th order IDC constructed using 3rd order ARK3KC with 0, 1, and 2 correction loops (and 3, 6, and 9 quadrature nodes), respectively. The order of accuracy is clearly seen, as the dotted reference lines (with slopes of 3, 6, 9) indicate. 109 10 +IDC4, y1 N '6‘", y2 >{_ _5 +|DCB, y1 a 10 -.u yz e e 8 2 10‘”- (U 10’4 10'2 H (a) VdP (3.36), e = 1 (nonstiff) .~ 9‘] >5" 100 >1- 9 e 5 ,. a) a -10 10 " -A-", y2 10‘4 10‘3 10'2 H (b) VdP (3.36), e = 10‘3 (mildly stiff) Figure 3.5: Convergence study of absolute error at T = 4 vs H, for 4th and 8th order IDC constructed using 4th order ARK4KC with O and 1 correction loops (and 4 and 8 quadrature nodes), respectively. The order of accuracy is clearly seen in (a), as the dotted reference lines (with slopes of 4, 8) indicate. 110 10 . +1003, y1 N _e_u’y2 >§_ -5 --Ioce,y1 3‘ 10 i---" y2 g +IDCQ,y1 t 'V‘",y2 2 10‘”. (U ’9’ ,0 or 10‘6 (a) IVP 3.35, e = 10”3 (mildly stiff) abs errors: y1, y2 (b) VdP (3.36), e = 10‘3 (mildly stiff) Figure 3.6: Convergence study of absolute error at T = 4 vs H, for 3rd, 6th, and 9th order IDC constructed using 3rd order ARK3KC with O, 1, and 2 correction loops (and 3, 6, and 9 quadrature nodes), respectively. The dotted reference lines (with slopes of 3, 6, 9) indicate the expected order of the methods, but note order reduction. 111 instability for larger timesteps; however, once the stepsize is small enough, order reduction is apparent, and then as the stepsize is further decreased, it appears that the expected IDC order is achieved (Figures 3.6b and 3.5b). The better behavior of solutions to the IVP (3.35) is not surprising, since many of the ARK methods were designed to solve problems similar to (3.35) (see, e.g. [56, 45, 21]), but the oscillator problem (3.36) has recurring stiff peaks, which are difficult to handle even with most popular integrators. Efficiency To measure efficiency, we plot the number of implicit function evaluations versus the error. As in the convergence plots, the error is found by comparing successive solutions (absolute error at the final time). The number of implicit function evalu- ations is counted numerically by tracking the number of times the “stiff” function and the Jacobian are called, including the Newton iterations. In general, the higher order IDC-ARK methods are more efficient than popular ARK methods. We see two general “rules of thumb.” First, if the order of the method is higher, then fewer function evaluations are required to reach a low error tolerance. This trend can be observed for ARK methods (equivalent to a prediction loop) (Figure 3.7a) and when more correction loops are taken (Figure 3.8a). Second, the efficiency of an IDC-ARK method generally depends upon its con- stituent ARK method; i.e., if ARK method a is more efficient than ARK method b, then IDC constructed using ARK a is more efficient than IDC constructed using ARK b. For example, an IDC method using the 3rd order ARK3KC frequently is more efficient than the other methods we tested (see, e.g., Figure 3.7b), but IDC using ARK2A1 performs relatively inefficiently (as expected since ARK2A1 alone requires the highest number of function evaluations among the ARK methods we tested in Figure 3.7a). We note that the behavior for the nonstiff (c = 1, Figure 3.8a) 112 # implicit function evaluations I“, D D 0 E 10- 10' L... of y1, y2 abs errors (a) VdP (3.36), e = 1 (nonstiff) 105. . ; +Ioce-3Kc i ' +IDCG-3BHR13 -e- IDCQ—3KC j -v- lDCQ-3BHR1 - -‘- -- -s Q ‘~ -. # implicit function evaluations 8 1 L... of y1, y2 abs errors (b) IVP 3.35, c = 1 (nonstiff) Figure 3.7: Efficiency studies of # implicit function evaluations vs absolute error at T = 4, using (a) ARK methods of second (2A1 [56], 2ARS [3]), third (3BHR1 [8], 3KC [45]), and fourth (4A2 [56], 4KC [45]) orders (no IDC). (b) for 6th and 9th order IDC constructed using 3rd order ARK3KC and 3rd order ARK3BHR1 with 1 and 2 correction 100ps (and 6 and 9 quadrature nodes), respectively. 113 -v-|DC9 "Q. 2 --<>-IDC3 . ‘° 10“" 10' L... of y1, y2 abs errors # implicit function evaluations (a) VdP (3.36), e = 1 (nonstiff) +|DC6 -v-|DC9 4 "0"IDC3 I; 10-10 -5 L... of y1, y2 abs errors # implicit function evaluations (b) VdP (3.36), e = 10‘3 (mildly stiff) Figure 3.8: Efficiency study of # implicit function evaluations vs absolute error at T = 4, for 3rd, 6th, and 9th order IDC constructed using 3rd order ARK3KC with 0, 1, and 2 correction loops (and 3, 6, and 9 quadrature nodes), respectively). 114 and mildly stiff (c = 10_3, Figure 3.8b) cases is similar. 3.6.2 Advection-Diffusion Example We now consider an advection-diffusion equation Ht = -—’U.1f + 11.55113, t E [0, 0.1], .73 6 [0,1], (3.373.) u(x, 0) = 2 + sin(47rx), (3.37b) with periodic boundary conditions. We solve the advection-diffusion BVP (3.37) via method of lines with fast Fourier transform (FFT) for the spatial derivatives and time-stepping with semi-implicit IDC-ARK, where the advection term is treated explicitly ( f N = —u;1;) and the diffusion term is treated implicitly (f3 = umy). This choice of methods ensures that the error is dominated by the time-stepping. We expect that treating the diffusion term implicitly will remove the requirement that the CFL condition must satisfy At g ch2, and that the CFL condition will be controlled only by the advection term so that At S ch. To verify this hypothesis, we plot the error for IVP 3.37 at T = 0.1 versus H 2 At = 0.5Ax, timestepping with IDC constructed with 3rd order ARK3KC, in Figure 3.9a. The error is the L1 (spatial) norm of the absolute errors at the final time and was computed by comparing the numerical solution to a reference solution computed with small Ax and At. Shown are IDC methods with zero and one correction loops (3rd and 6th order convergence, respectively). Since the expected convergence orders are seen when At is proportional to Ax, we conclude that it is reasonable to say that treating the diffusion term implicitly improves the CFL condition to At g ch. This improvement is exactly what we expect for a semi-implicit method. 115 _ +1003 10 -->-IDCB -10_ 10 abs error 1 I" ,n .- . 10-15r 10 10'2 H (a) Advection-diffusion BVP (3.37) Figure 3.9: CFL study showing absolute error at T = 0.1 vs H 2 At = 0.5Ax. FFT is used for the spatial discretization, and time-stepping is done with 3rd and 6th order IDC constructed with 3rd order ARK3KC with 0 and 1 correction loops of IDC, respectively, where the advection term is treated explicitly and the diffusion term is treated implicitly. This choice of methods ensures that the error is dominated by the time-stepping. Since the expected order of accuracy in time is seen, as the dotted reference lines (with slopes of 3, 6) indicate, it appears that the CFL condition is At 3 ch. 116 3.7 Concluding Remarks We have provided a general framework for the simple construction of arbitrary order IMEX methods. These semi-implicit IDC-ARK methods use ARK integrators as base schemes inside the IDC framework and can thus be constructed to higher order easily, with no need to consider the complicated order conditions that are typically required for the construction of ARK methods. High order IDC-ARK methods can be used to solve multiscale differential equations that involve disparate time scales more efficiently than popular ARK schemes when a low error tolerance is desired. We have performed an analysis of the local truncation error of IDC-ARK methods, shown improved stability over IDC-FEBE methods, and conducted numerical results showing order of convergence, efficiency, and the potential for an improved CFL condition. We plan to examine the use of asymptotic preserving ARK schemes (as in [64], and Section 6.1) in the prediction and correction loops of IDC as a potential way to mitigate the issue of order reduction that arises with increased stiffness. We also have begun an investigation of embedded IDC methods, including IDC-ARK, and their implementation in an adaptive setting (see Chapter 4 for preliminary results). We anticipate that adaptivity will provide additional efficiency gains over IDC—ARK, and in fact, adaptive timestepping is the natural way to handle many multiscale problems. 117 Chapter 4 Adaptive IDC Methods 4. 1 Introduction In this chapter, we consider initial value problems of the form y' = f5(t»y) + fN(t.y). y(0) = yo, te [01'1"], (4-1) where y 6 IR”, and the function containing stiff terms, f S(t,y) : IR x IR" —> IR", and the function containing non-stiff terms, f N(t, y) : IR x IR" ——> IR”, are Lipschitz continuous in the second arguments. When f S = 0, i.e., equation (4.1) contains only non-stiff terms, explicit integrators such as popular explicit Runge—Kutta (RK) methods or Adams—Bashforth methods can be applied. When f N = 0, i.e., equa- tion (4.1) contains only stiff terms, implicit RK methods or Adams—Moulton meth- ods can be applied [2]. When equation (4.1) contains both stiff and non-stiff terms, implicit—explicit (IMEX) methods [4] or additive Runge—Kutta (ARK) methods can be applied [21]. Although in previous chapters we only consider fixed step sizes, to guarantee 118 accuracy in the solution of equation (4.1), estimation of the error and adaptive step size control is often desirable. For example, in Van der Pol’s oscillator (3.36), there are stretches of time where the solution changes only mildly, punctuated by oscillatory stiff peaks in the derivative. The stiff oscillations cause methods imple— mented with a fixed step size to exhibit order reduction, as explained in section 3.6.1. Adaptive time-stepping may be needed for other types of problems as well. For problems containing only non-stiff terms, a widely accepted approach is to use a matched pair of explicit RK methods [37], also known in the literature as embedded RK methods. The premise is to select a pair of RK methods (i) whose Taylor series approximations agree with the corresponding Taylor terms of the solution up to orders p and p + 1, and (ii) whose stage weights, aij, and stage nodes, cj, overlap as much as possible for efficiency. The matched pairs are often expressed compactly in a Butcher Tableau, where A is an s x 3 lower triangular matrix, b,b and c are s-vectors. An exam— ple of a pair of matched explicit RK method is the Runge-Kutta—Fehlberg Method (RKF 45), which generates a fourth and fifth order method using six stages. This ap- proach of matched pairs can be extended to implicit RK and ARK methods. In [33], matched diagonally implicit RK pairs are discussed, and [45] discusses embedded ARK methods and error control. In recent papers [18, 17, 15], the present authors have illustrated that integral deferred correction (IDC) methods are competitive with popular RK and ARK methods. In brief, IDC methods are a class of deferred (or defect) correction methods 119 that generate a sequence of approximations that successively reduce a residual. First introduced by Dutt, Greengard and Rokhlin in [25], high order IDC methods can be constructed systematically without solving tedious order conditions. This work exploits the framework of an IDC method, using computed approximations of the numerical error for step size control. In Section 4.2, IDC methods are reviewed, and numerical comparisons between IDC and existing embedded methods are given in Section 4.3. Conclusions are given. in Section 4.5. 4.2 Integral Deferred Correction IDC methods are essentially predictor-corrector schemes. In this section, we first discuss how the error for an approximate solution is computed. Then, the IDC framework and algorithm is given. Finally, an adaptive time step control is discussed. 4.2.1 Error Equation Given an approximate solution, 17(t), to the exact solution, y(t), of (4.1), the error of the approximate solution is 60) = y(t) — 77(t)- (4-2) Defining the residual, c(t), as ea) = n’ 1: parameter for coarsening time step (in Figures 4.2a, 4.2b, 4.3a, and 4.1 below, we use a = 0.5, 7 = 0.1, fl = 2) Output: At: new time step; accept: flag whether to accept new update; 1 if err > TOL then /* reject step */ 2 [ At = aAt, accept = 0; 3 else if err < 7 >1: TOL then /* coarsen step */ 4 [ At = fiAt, accept = 1; 5 else At = At, accept = 1 ; /* accept and leave unchanged */ Algorithm 2: Determination of step size For embedded RK (or embedded ARK) methods, two solutions of order p and p + 1, resp. are calculated. Then, for example, the error can be approximated by comparing the difference between those two solutions. Each correction iteration of an k—l]) IDC method includes an error estimator (the correction 6i , making it a natural choice for adaptive implementation. For IDC methods, the error approximation is 123 already built into the IDC method, since each correction 100p solves (4.4) for the error. If we consider that IDC-RK (IDC-ARK) methods are RK (ARK) methods, then we notice that embedded IDC-RK (IDC-ARK) methods are a natural analogue to embedded RK (ARK) methods. In particular, if the final correction of an IDC method is order q, and the solution just before the final correction is order p, then the embedded IDC method gives two solutions of order p and p+ q, resp. An obvious choice for the final correction is g = 1, to avoid extra computational expense when approximating the error. Note for IDC methods, At is adapted, while 6t 2 At /]W remains uniform for each choice of At. 4.3 Numerical Comparisons In this section, a numerical study is performed on such as Van der Pol’s oscilla- tor problem, a stiff ODE, comparing the performance of IDC-ARK methods and other popular embedded methods. In the figures and table in this section (with the exception of Figures 4.1a and 4.1b), we use the notation in [59] to compare various characteristics of the numerical methods and to facilitate future tests with additional IV P3 in [59]. The notation is as follows: 0 solver The numerical method with which the ODE was solved. a atol Absolute error tolerance provided by the user. a scd Minimum number of significant correct digits the numerical solution has at the end of the integration interval, i.e. scd = —log10||( relative error at the end of the integration interval )Iloo (4.6) 124 0 steps Total number of steps the solver takes, including steps that are rejected because of error test failures and / or convergence test failures. 0 accept Total number of accepted steps. 0 # imp f Number of evaluations of the implicitly evaluated function, f S (mod- ified from [59] to count only implicit function evaluations). Implicit function evaluations are counted by keeping track of how many times the Matlab script that evaluates the implicit part ( f S) is called. We use only implicit func- tion evaluations based on the assumption (as in [60]) that implicit calls will dominate due to Newton iterations, although this assumption is not always valid. 0 #Jac Number of evaluations of the Jacobian of f5. Note: The Jacobian is used in Newton iterations as part of the implicit method and is calculated analytically at every Newton iteration in every step, which is likely not the best way to handle the Jacobian. A smarter choice of Jacobian calculation should give far fewer evaluations. The tolerance for the Newton iterations is fixed at = 10—13, and the maximum number of iterations is 10. We also include the following additional characteristics with notation: o # refine Number of times the timestep is refined, over the entire calculation. 0 # coarsen Number of times the timestep is coarsened, over the entire calcu- Iation. o minstep The smallest timestep taken over the integration interval. 0 maxstep The largest timestep taken over the integration interval. 125 4.3.1 Van der Pol Oscillator We consider the application of adaptive semi—implicit methods to Van der Pol’s os- cillator, as described by equation (3.36), with varying levels of stiffness (stiffness is increased by decreasing c). As a preliminary test, we considered a recursive im- plementation of some adaptive IDC-ARK and ARK methods where the stepsize is only reduced when the error is too large and then reset to the original At for the next step. Algorithm 5 is far more practical, but we present the results from the recursive implementation as an example of the improvement that can be possible by using IDC methods. In comparison with embedded ARK methods, comparable or improved efficiency seems likely with appropriate choice of base scheme integrators. For example, I constructed adaptive IDC-ARK methods and found that fewer time gridpoints are required when using a 7(6) order embedded IDC—ARK method (Fig- ure 4.1b) rather than a 4(3) order ARK method (from [45]) to solve Van der Pol’s oscillator problem (Figure 4.1a) to the same error tolerance. The 7(6) IDC-ARK method solves the prediction and the first correction loops via a 3rd order ARK method from [45] and the final correction approximates the error via semi-implicit first order forward and backward Euler. Now we examine the results from adapting the time step via Algorithm 5. In the figures, tol represents the quantity explained below as atol, and is the tolerance for the error when choosing whether or not to refine the timestep. The number of implicit function evaluations is described by # imp f below and is the number of times the Matlab script calls the stiff function f S (which is evaluated via the implicit part of the numerical method). From Figures 4.2a, 4.2b, and 4.3a, we see that, for any choice of tolerance, ARK4(3) is more efficient than the same order IDC method constructed via a semi-implicit forward and backward Euler combination, IDC4(3)FEBE. This conclusion holds for the nonstiff (c = 10_1), mildly stiff (e = 126 0 1 f 3 Z (a) ARK4(3)KC 4 o y1 o y2 O . . . . 0. 00000 -2. '40 1 é 3 it t, including IDC subnodes (b) IDC7(6)-ARK3KC/FBE Figure 4.1: Preliminary results comparing (4.1a) an adaptive ARK method (embed- ded order 4(3) from [45]) and (4.1b) an adaptive IDC method constructed with a 3rd order ARK method from [45] (both recursively implemented). Both plots show solution to Van der Pol’s oscillator (c = 1) for the same recursion error tolerance, 10—5. 127 VdP, |C=[2,0], Tf=2, NN=1,£ = 0.1 +ARK4(3)KC , -- lDC7(6)-ARK3 : --|DC7(6)-ARK2, + lDC4(3)-FEBE§ +lDC7(6)-FEBE1 _L O 0’ —L # implicit func evals 10 _ ‘- 10 1° 10 5 10° tol (a) c = 0 l 7 VdP, |C=[2,0], Tf=2, NN=1,e = 0.001 10 . +ARK4(3)KC , -- IDC7(6)-ARK3 : -I-lDC7(6)-ARK2, + lDC4(3)-FEBE§ +lDC7(6)—FEBE1 # implicit func evals 1130'” 10‘5 10° tol (b) c = 0.001 Figure 4.2: Efficiency study comparing semi-implicit ARK and IDC-ARK methods used to solve Van der Pol’s oscillator. Plots show # imp f vs atol, as described in the text. 128 VdP, |C=[2,0], Tf=2, NN=1,e = 1e-06 +ARK4(3)KC , -- lDC7(6)-ARK3 : +|DC7(6)-ARK2, +IDC4(3)-FEBE§ +IDC7(6)—FEBEI # implicit func evals ES Figure 4.3: Efficiency study comparing semi-implicit ARK and IDC-ARK methods used to solve Van der Pol’s oscillator. Plots show # imp f vs atol, as described in the text. 129 8.83 AAtmws «an». AAEN Ewan.“ sewn” ENE $38 2.2 3: 9:; 3:55 22$ 8 32 Emma EAAmm 82A SSA 8.2 8.2 Exam? egos 8133s As was A88 Sam ABA AEA A83 A3 oxAmVA. 3:2 82 was”? a. 2: Es Anew man was. AAV 8A oAm V $3. 383 AAts..A 28 ES @389. 5me $82 Eta oAtA 33A mwAmaAmVAooA 38d 223 As a: 83% $33 Ame: 83A 8.2 Ase mmmaAm V59 Sec 8.3”? 2 3o 83% $82 8% 88 88¢ 83. mmmaAmVAooA 83 $th 3 SA 32: AwmAm an. new A. AV men mmmaAmVAooA 8.3; 2.23 AAA 32 gnaw 839: mg News 2.2 was mmmas V59 82 moses 8 as 882 $22 as «2: 8.2 we... mmmaa V59 83 8138A 8 «AAA 55% Emma a: an 88.9 was mamas V59 82. SANA. on Am 83m 33. N: 8m AAA m3” mmmas V89 Asses 2.83 mam Atom: 388: 8E5 38A 88A 22 3: mmfixmfio V59 AwBAVAV moses Am an. £83 88mm 3: on: 813A New mmfixmioroe 83 Sam: 8 8A can. 82w man as” 82; Es mm<§m (5-12) Haj—H) if f (u) < 0 5.3. 1 Reconstruction For higher than first order schemes one needs a process called reconstruction. The re- construction problem is: given the cell averages fii, reconstruct the solution u(xz+1/2) at the cell boundary 1241/2 with higher order. For the first order stencil, we use only the immediately neighboring cell, either 1,; or Ii+1' More neighboring cells are used to achieve a higher order stencil. In general, let Sr = {Ii-“7”""Ii""’Ii—T+k—1} (513) 144 denote a stencil, where k cells are used for the stencil, and r E {k—l, . . . , 1,0, —1}. In the hyperbolic setting, which we delineate here, one must always use an odd number k of neighboring cells because the characteristics have a direction. One chooses more neighbors from the direction from which the characteristic comes. In a hyperbolic setting, choosing an even number of neighboring cells results in an unstable method (e.g., centered stencil). Under the right conditions, the order of accuracy of the reconstruction equals the number of neighboring cells used in the stencil. This process is very similar to interpolation, e.g., for a 3rd order reconstruction where k = 3, we find the polynomial p(:r.) = a0 + (11:1: + a2x2 such that p satisfies 1 (1 _ _ A: Ii-1p(.r)..;r — "(1.2-_1, (5.14) ~1— .‘ 1:17 : '17- (5.15) A2: 12' MW 2, i as: = 2:.- . (5.16) A9” [2+1 12(2) 2+1 Then, since a0, a1, and a2 are known, we can estimate 2125+”? from the left as ui_+l/2 = p(:ri+1/2) (5.17) when p is constructed from a left stencil, S1 = [1,-_1, 1,, [241}. Similarly, 11;: 1/2 = p(x'i+l/2)’ when p is constructed from a right stencil, 50 = {I.i, 1,41, 1242}. t ' ' . 7L ' Two me beds of reconstructing solutions uz+1/2, uz+1/2 at the cell boundaries for a higher order approximation scheme are EN 0 and WENO methods. 5.3.2 ENO When the solution It contains discontinuities, using all In cells in a stencil 5,» gives oscillatory behavior near the discontinuities. Consider the case when k = 3. Away 145 from a discontinuity, we use three cells for the stencil, e.g., 51 = {I j_1, Ii, Ii +1}, and achieve a third order approximation. Near a discontinuity, e.g., at :13,- +1 /2, where the characteristics move from left to right, instead of using three cells as in 5'1, consider two smaller overlapping stencils 3'1 = {[,_1. It}: 350 = {1751,41}. Note that stencil :51 will not cause oscillations in a reconstructed solution because it contains cells only on one side of the discontinuity, whereas stencil §0 contains the discontinuity and thus will result in oscillatory behavior. Hence in the case |fi,,j_1 - 11,] << Ia,- — '17., +1], ENO will choose stencil 51 as the “better” stencil. Note that choosing these stencils near the discontinuity will result in a lower order approximation, but only near the discontinuity. Away from the discontinuity, the approximation remains fully third order. In general for FV methods, ENO uses a linear combination of the cell averages '11,; to reconstruct the solutions “241/2 at the cell boundaries. ENO with FD methods is similar, except the numerical flux (rather than the solution) is reconstructed at the cell boundaries. The details are explained below for WENO methods. For an ENO reconstruction method that uses It cells, the method is order k away from the discontinuities, and lower order at the discontinuities. 5.3.3 WENO WEN O is similar to ENO, except WEN O uses weighted combinations of the sten- cils to achieve higher order approximations. A WEN O method where each stencil contains k cells has order 21: —- 1 away from discontinuities and order k near dis- continuities. We present the procedures for WEN O methods in detail for the FD case involving one—dimensional scalar flux splitting, followed by an extension to 1-D systems. 146 l-D Scalar Case First, we present the following outline, given by Procedure 2.5 and 2.2 of [70], of the spatial approximation of equation (1.12). 1. Since the numerical flux in equation (5.9) should increase in its first variable and decrease in its second variable, the first step is to form a smooth flux splitting M) = Na) + flu), (5.18) where (If+ (1f— For example, the Lax-Friedrichs flux given by Wu) = (f(u) + an), (5.20) f-(u) = (f(u) — an), (5.21) [\DIr—INIr—I ax where a =mu, [f’(u)|. 2. Identify the cell average of some function v(;r) at 1:3,; with one part of the split flux, i.e., set 27,- : f+(u,-). (5.22) Then reconstruct the cell boundary values from the left, vi+1/2, for all 2', e.g., via a WENO method as in the subsequent procedure. 3. Set the positive numerical flux equal to the reconstructed value from the pre- 147 C51 vious step: [111/2 = 57.11 /2‘ (5.23) Now set = flu.) (5.24) and reconstruct the cell boundary values from the right, val/2, for all i, e.g., via a WENO method as in the subsequent procedure. Set the negative numerical flux equal to the reconstructed value from the previous step: A ,— _ + f5+1/2 — vi+1/2' (5‘25) Form the numerical flux from the positive and negative fluxes at the cell bound- aries: f”Pd/2 Z f211/2 +f511/2: (5'26) i.e., fi+1/2 = Dill/2 + pal/2' (5.27) Form the ODE dug—ft) = “Al—$(fi+1/2 — f3141/2): (5'28) 148 to be solved by the desired numerical integrator. Now we present the 1-D WENO reconstruction procedure from [70]. Given the cell averages i),- of some function v(:r) for each cell 127) the goal is to obtain upwind- based (2k — 1)th order approximations, 17+ and v.— to v(:r) at the cell 2—1/2 2+1/2’ boundaries. 1. First find the k reconstructed values (7‘) k—l Ui-l-I/Q = aCTjEi—T‘l'j’ (5.29) J: j come from reconstructing the polynomial whose cell averages coincide with the cell averages of U (see, e. g., Section 5.3.1 for more de— whcre the coefficients (3,. tails on the case k = 3), and the stencil Sr = {Ii—r, li_,.+1,. . . , li—r+k—1} containing k cells is used, for 7‘ = 0, 1, . . . ,k — 1. Also obtain the k reconstructed values k—l ,(7') _ ~ .—,. ti—l/2 — Z 07'] L,_,.+J- (5.30) 1:0 k—l : Z Cr,j'l7i_(r+1)+j, (5.31) i=0 where Erj = Cr—1,jv over the stencil Sr, for -r = —1,0, 1, . . . ,k — 1. 149 2. Find the constants dr and d7. such that ' — k_ld , (7') 5 32 7"2'+1/2 — Z ‘7‘”’.'+1/2 ( ° ) r=0 = v(a:z-+1/2) + 0(Azr2k—l), (5.33) - — k—lci , (7”) 5 34 Ill-_l/Q — Z TU‘li—l/2 ( . ) r=0 = v(.2:z-_1/2) + 0(Ax2k—l). (5.35) Note that (fr 2 dk_1_,.. 3. Find the smoothness indicators k4 21 1 011207) 2 Hr: E An: “ ——"‘ dz, (5.36) for upwinding in the left to right direction. For upwinding in the opposite direction, modify the procedure symmetrically with respect to 2241/2 to find the smoothness indicators Br- 4. Form the weights wr and LOT: .... — 25:; as, a. (f +61%? (5.37) (Dr = fi, 517- : (“de (5.38) forr=0,l,...,k.—1and0 to + At) to solve: .- 8tf+ A317]: = O, f(.’l?, v, to) = f0(:l3,v), (5.73) Since the exact solution is known, one may simply find the solution as A. f(:L", v, to + At) = f0(x — AAt, 1)). (5-74) 157 0 Take another full timestep (to ——> to + At) to solve: ~ ~ ~ atf+ BO-Uf: = 0, [(131), to) = f(:1:,v,t0 + At). (5.75) Similarly to above, one may find the solution as A. .- f(.r, v, to + At) = fix, '0 — BAt, to + At). (5.76) Then fz is a first order solution to (5.72), if the operators A83; and Bay do not commute, and if (5.73) and (5.75) are solved with a numerical method that is at least first order in time. I.e., the error in the solution is due to both the splitting and the choice of numerical method(s) used at each step of the splitting If A83; and Ba.) do commute, then there is no splitting error. [73, 24]. However, for this splitting, note that each subproblem is a 1D linear advection equation, which can be solved exactly via a semi-Lagrangian method, so the only error is the splitting error. The popular example of a 2nd order splitting method is Strang splitting (other names include Marchuk splitting, fractional step method) [13, 73, 58]. We describe below one timestep of size At, from to to to + At. 0 Take a half timestep (to -—> to + At/2) to solve: am + Aaxfi = 0, [*(x, v, to) = f0(:c, v). (5.77) Since the exact solution is known, one may simply find the solution as f*(.r, v, to + At/2) = f0(:1: — AAt/2, v). (5.78) 158 0 Then take a full timestep (to —> to + At) to solve: 8tf** + Bavf“ = 0, f**(:1:,v, t0) = f*(:r, v, to + At/2). (5.79) Similarly to above, one may find the solution as f**(a:, v, to + At) = f*(:1:, v — BAt, to + At/2). (5.80) 0 Finally, take another half timestep (to + At /2 —+ to + At) to solve: Utfspm + 43wa778 = 0: fsplirlxs 11, t0 + AU?) = f ”(17: v: t0 + At)- (581) Similarly to above, one may find the solution as fame, “U: t0 + At) = f *(r — AAlt/2, 22. to + At). (5.82) fsplz't is a second order solution to (5.72), if the operators A83; and Bay do not commute, and if (5.73) and (5.75) are solved with a numerical method that is at least second order in time. However, for this splitting, just as for the first order splitting, note that each subproblem is a 1D linear advection equation, which can be solved exactly via a semi-Lagrangian method, so the only error is the splitting error. As mentioned above, commuting operators will not give a splitting error [73, 13]. 159 5.6 IDC with Splitting Methods In this section, we expand the results from [67] to examine the same spatial dis- cretization coupled with higher than 2nd order time splitting, using IDC methods to increase the order of low order splitting methods, and in following sections, we examine the effectiveness of our novel construction in solving equations of the form (5.72). Some aspects of our work are similar to constructions in [46, 49, 50]. Note we introduce a special change for the IDC method’s error equation for the Vlasov system, but for constant advection, rotation, and drifting rotation problems, the error equation is the same as in previous versions of IDC. 5.6.1 Overview of IDC methods The basic construction of IDC methods is the same as for SDC methods [25, 18, 15]. Suppose we wish to solve the 1D version of (5.1) (given below by (5.88)) to the final time T. The time interval, [0, T], is discretized into intervals [tn, tn+1], n = 0, 1. . . . , N, such that 0=t0 to + (it) to solve: t até = —7‘ = —8t (77 — f0 '1‘] (148177 + 88117))(17'), IC -'= C(LO), (5.116) 0 Take a full timestep (to —> to + 6t) to solve: a5+8a5=o,1c25, (5mm 0 Take a full timestep (to ——> to + 6t) to solve: 8t68p11+ Aaxespll = 0, IC = (:3. (5.118) Alternatively, a second order splitting of the error equation (5.115) could be given as [29]: 0 Take a half timestep (to —> to + cit/2) to solve: (7)731 +./l(’);17(:1 = 0, IC = C(10), (5.119) 0 Take a full (to —> to + (it) timestep to solve: - Take a half timestep (to ——> to + cit/2) to solve: 167 — Take a full (to ——> to + (5t) timestep to solve: 0,53 = —r, 10 = 52(10 + Lit/2), (5.121) — Take a half (to + (it/2 —+ to + 6t) timestep to solve: 8te4 + 381L164 = 0, IC = e3(t0 + (it), (5.122) 0 Take a half (to + (it/2 —> to + 6t) timestep to solve: atesplg + A8336 2 = 0, IC = e4(t0 + 5t). (5.123) spl 5.7 Efficiency Analysis The following sketch compares the number of split equations that must be solved by high order splitting methods constructed via IDC methods and those methods constructed via Yoshida’s methods in [77]. Note that splitting methods in [39] scale similarly to methods in [77]. Although the metric we choose below is an imperfect metric for comparison, it suggests that there is merit in considering high order split IDC methods. In the comparison we assume the timestep At of Yoshida is the same as (it 2 At / M (M = number of substeps in each IDC timestep) of an IDC method; i.e., for Yoshida, we let T = At, and for IDC methods, we let T = 6t. Consider a pth order method solving 8,5 = (A + B)u, (5.124) with exact solution exp(T(A + B)), after one timestep. In keeping with Yoshida’s notation, a splitting method’s solution S(r) can be written as a composition of 168 exponentials k 3(7) = H exp(e,~rA)exp(d,rB) (5.125) i=1 = exp(T(A + B) + 0(7—P+1)). (5.126) \Ve introduce the idea of counting the number of split solves as a measure of the effi- ciency of the method, with the assumption that each split solve is roughly equivalent in cost (understanding that each split solve could be significantly different and thus nullifying our comparison). For example, a first order splitting could have (:1 = 1 and d1 = 1, to give 51(T) = 674.373, (5.127) which would require two split solves. Similarly, a second order splitting (equivalent to Strang splitting) could be written with c1 = 1/2, d1 = 1, c2 = 1/2, and d2 = 0, to give 17A B 17A 32(7) = e? 67 e? , (5.128) which would require three split solves. For each of Yoshida’s splitting methods that are higher than 2nd order, one must solve a set of equations (not shown here) to find the coefficients c,, 11,-. To find the exact coefficients, one must solve the coefficient equations analytically. In this situation, if the order of the method is p (and p is even), then the number of nonzero coefficients is 2k — 1, where k = 313/271 + 1. These coefficients will be referred to as the analytic coefficients. To find a simpler higher order splitting method, one may solve an approximate set of coefficient equa- tions, giving fewer coefficients, but they are approximate and not exact coefficients. 169 These coefficients will be referred to as the nonanalytic coefficients. In practice, the nonanalytic coefficients are used since they make the computation much more cost effective (see Table 5.1). The number of nonanalytic coefficients is also given by 2k — 1, but 16 is different: I: = 21 + 2, where t = 2p/2-1 — 1, and the (even) order is p > 4. Some of these calculations can also be clarified by examining [77, 49, 50]. The number of nonzero coefficients (analytic or nonanalytic) can be thought of as the number of split solves required for one timestep of that method. order Yoshida analytic Yoshida nonanalytic IDC-82 IDC—Sl 2 3 - 3 3 4 7 - 7:1,- 9 6 19 15 12§ 15 8 55 31 17% 21 Table 5.1: Number of split solves required for one step r of each splitting method. Yoshida’s methods are from [77], and r 2 At. IDC-S2 refers to split IDC methods constructed with second order splitting, and IDC-SI refers to split IDC methods constructed with first order splitting. 'r 2: 6t for IDC methods. Unlike Yoshida’s methods, it is not entirely correct to write split IDC methods as compositions of exponentials. However, we can still count the number of split solves in one timestep in a similar manner. Recall we are considering the solution of (5.124), so there will be no partial derivatives in the following analysis. One substep r = 6t of the split IDC prediction loop is simply a basic splitting method, which can be written as (5.127) or (5.128), for example. Thus the cost for one substep of the prediction loop is 2 solves for a first order splitting or 3 solves for a 2nd order prediction. One substep of the correction loop, solved according to the first order 170 splitting given by (5.116)-(5.118), can be written as (11(7) = 5746731117), (5.129) where R(r) is the solution to (5.116). Thus the cost of this correction substep is 3 solves. However, if we look carefully at Algorithm 4, line 6, we see that the very first substep of an IDC interval, from t = to to t1, only requires the solution of (5.116), but (5.117) and (5.118) are unnecessary. Thus the first substep of the correction loop only requires I solve, while the last M — 1 substeps require the full 3 solves seen in (5.129). Similarly, a second order splitting given by (5.119)-(5.123) can be expressed as 1 1 1 1 02(7) = 62”.;273 12(552736274, (5.130) where R(r) is the solution to (5.121). Thus the cost of this correction substep is 5 solves. Again, since the initial error at each correction is zero, we may adjust the second order splitting to be more efficient (similarly to Algorithm 4, line 6). Thus we do not need to solve (5.119) or (5.120) over the first substep, but only (5.131) Therefore the first substep of the correction loop only requires 3 solves, while the last M - 1 substeps require the full 5 solves seen in (5.130). To coalesce the number of solves required for an IDC prediction and the number required for a correction into the total number of solves required for one substep 6t, we must first consider one step of size At. For a pth order IDC method, we subdivide At into M = p - 1 substeps (see (5.83), (5.86) in Section 5.6.1). We solve for the prediction solution over each of those M substeps. Then we solve for 171 the error correction J times over each of those M substeps. Letting P denote the number of solves required for a prediction substep, 30 the number of solves required for the first substep of a correction, and s the number of solves required for the subsequent substeps of a correction, the resulting number of solves over At is MP + J(50 + (M —1)s). (5.132) Then we divide by M to obtain the number of solves required for one substep T 2 6t, solves = P + J(sO + (M — l)s)/M. (5.133) For example, if we consider a fourth order split IDC method constructed with first order splittings, then one must solve one prediction and 3 correction loops, requiring 9 solves per substep: solves=2+3(l+2*3)/3=9. (5.134) For a fourth order split IDC method constructed with second order splittings, one must. solve one prediction and one correction loop, requiring 7:1; solves per substep: 1 solves =3+1(3+2*5)/3=7§. (5135) See the Table 5.1 for other higher order results comparing Yoshida’s split methods and split IDC methods. For 2nd, 4th, and 6th order methods, Yoshida and IDC methods appear roughly equivalent, whereas for 8th order (and higher, though not shown), IDC methods appear to have an advantage in requiring fewer split solves. 172 5.8 Conservation of Mass In this section, we prove that the split IDC methods constructed in this work pre- serve total mass in the case of periodic boundary conditions; i.e., for ftnj+1 the numerical solution to the Vlasov-Poisson system at t = Un+11 :1; = 1:,- = x0 + iAx, v = vj = ”0 + j Av, we desire the discrete version of :r.v,t d2; dv =/ f f a:,v,t dxdv 5.136 /X /V f( n+1) X V ( n) ( l to hold: N1 —1Ng:— 1 1 qu-‘lNgy—l 71+ ,. , _ n :0 :6 1,3. ALAv— Z 2 mm. (5.137) .7 i=0 For simplicity. below we drop ArAv and the limits on the sums. Theorem 5.8.1. A split IDC method that is constructed via first order splitting methods in the prediction and correction loops together with conservative semi- Lagrangian WEN 0 methods, and whose spatial and velocity derivatives in the resid- ual are approximated via WENO reconstruction, conserves total mass if periodic boundary conditions are imposed. Proof. By Proposition 5.8.3 and repeated application (according to the number of correction loops) of Proposition 5.8.5, both given below. E] Proposition 5.8.2. 1;?“ = f," — 100,711 flee) — 13:1 flee» (5.138) conserves total mass if periodic boundary conditions are imposed and f," 1/2 is the numerical flux given by a 2k + 1 order conservative semi-Lagrangian WENO 173 reconstruction. For example, when interpolation and reconstruction is from the left, the flux: is defined as “n _ n n L 2k. for C2Lk+1 a matrix of coefficients and 5 E [0,1/2]. Further details are in [67]. Proof. The sums are over i = 0, 1, . . . , N — 1 (i = N is excluded due to periodicity). Z r?“ = Ev," — €00?“ flee) — 1,711 /2(€0))) = Zfzn — {0 Z(f:::_1/2(€0) _fA,‘n_1/2(€Ol) = 2 f," - {of-f31/2(€0) +f)7(;_1/2(€0)) =21", i where the last equality holds from periodicity, which can be seen by writing out the fluxes as defined in the proposition. D Proposition 5.8.2 represents the solution of only one of the split equations in (5.94). However, conservation also carries through for consecutive splittings. Proposition 5.8.3. The first order split prediction of a split IDC method solved via conservative semi-Lagrangian WENO conserves total mass if periodic boundary conditions are imposed. Proof. Let f 3‘- denote the solution updating 7‘. in the x direction, given by 2] U] 1;; = 13’;- —€0(f;:1/2,j(€0) — 131/Ma». (5.140) 174 We know by Proposition 5.8.2 that :1 __ n 2],,- _ 2ij (5.141) 2 i for all j. Then let fig-+1 denote the solution updating f;- in the v direction, given by +1 5 * f‘n' 2 f2? — €0(f:j+1/2(€0) “ fitj—l/QKOD' (5142) ”(J Again, by Proposition 5.8.2, we have n+1 _ =1: 2:ij —Zfz'j (5.143) j i for all i. Therefore F‘- (5.144) =2 2);} =2 2f};- . (5.145) 1 j 2' Lemma 5.8.4. If the derivatives are found via WENO reconstruction, then 2211-8 ~+EU+€a ~—0 (5146) : ' J 23% ,- 21772,] — - - Proof. From [67], we see that the WEN O flux coefficients can be written as a vector, which for example, when using a 2k + 1th order reconstruction from the left, we 175 denote as 1021A: +1. Then the numerical flux P,” can be written —1/2 A n _ n .n ,L Fi—1/2—(fi—k—1"”’ 1+k—ll'u2k+1’ (5147) and a derivative can be approximated as l A A r. n ~ _ ,n _ 11 Now consider the sum in the lemma, with the derivatives approximated as in (5.148). 77+e ZZUJUI’Ivt +Ei Uv’h‘j 2 =2% 2: A1x( Fi+1/2,J' Ft-1/2,j) n+6 jn _‘n +ZE W2 A1U(Fi,j+1/2 Fi,j—1/2) _ _ p77,, “n An ‘ Z A, (‘EL,—1/2 + ELM—U2) where the last inequality holds by periodicity, which can be seen by writing out the fluxes as defined above. [:1 Proposition 5.8.5. The first order split correction of a split IDC method solved via conservative semi-Lagrangian WENO methods (and whose spatial and velocity derivatives in the residual are approximated via WENO reconstruction) conserves total mass if the prediction conserves total mass and periodic boundary conditions are imposed. 176 ‘_ Proof. Since (5.117) and (5.118) are solved in the same way as the prediction, we only need to prove that the numerical solution for (5.116) conserves mass. Denoting e"+1 as the error at the updated time Un+11 we have, as the solution to (5.116), t ,n+1_ ,n , n+1 n n+1 __ 477+e. .. (ll-j — (lj — 17.. + 7],,‘j — £7, ’1‘] (17:71,] (T) + El: (2117],? (T) (17'. Then 224'.“=zzei—Zzn:;+l+:zni i j i J i j i j t ._ Z Z]; n+1 i1ji);rr),:j(r) + Eg+€fh1r),-j(r)dr i j n t. = 2:62, -[ n+1 221583377000 + E27+68vn,-j(r) dT 2' 1' U 2‘ 1' =23), , J where the second equality holds by the hypothesis, and the last equality holds by Lemma 5.8.4. Cl Although the results in this section are for split IDC methods constructed via first order splittings, we anticipate that conservation of mass for split IDC methods constructed via second order splittings can be similarly shown. 5.9 Numerical Results We apply the split IDC methods to constant advection, rotating problems, and Vlasov-Poisson problems, where we study classic plasma problems, such as the warm two stream instability and Landau damping. Although some of the following notation used may be obvious, we provide it 177 here for complete clarity: At is the timestep as given by (5.84), N33, NU are the number of space and velocity (or second space, where relevant) steps. Tf is the final time to which the solution was calculated. When explaining which part of the operator is used in the splitting, we use the notation of A, B as in Section 5.5. For spatial differentiations and interpolations, we use WEN O reconstruction [19] and conservative semi-Lagrangian W EN O [67]. For error plots, the error is given as the sum of the l1 and loo norms (l1+loo). Comparison to exact, reference, or successive solutions is clarified for each test problem. 5.9.1 Constant Advection Here we apply split IDC methods to a constant advection equation, Ht + U1; + U1) 2 0, (IE,L‘) 6 [—1,].] X [—1,1], t) 0, (5.149) u(O, x, v) = sin (47r(x + v)), u(t,1,v) = u(t, —1.v), u(t,:r, l) = u(t,x, —1), using simple numerical integrators in each split step rather than exact (semi-Lagrangian) solves. We verify that the split IDC framework attains the expected increase in or- der of accuracy with each correction loop. Here we use A = B = 1. For the spatial derivatives, we use fifth or ninth order WENO reconstruction. In Figure 5.1a, we show error plots from solving (5.149) via split IDC that incorporates first order splitting where each split step is solved with a forward Euler (FE) integration. In Figure 5.1b, we show error plots from solving (5.149) via split IDC that incorporates second order splitting where each split step is solved with a second order Runge—Kutta integration. Clearly the order of accuracy increases and the error decreases with each successive correction loop as expected. 178 adv2D: split IDC, w5, Nx=Nv=100 0 10 E __-.r—:"'fi x N -_*_-—l'"'* V d r""r- ,.-‘ '. II "'7' v" -'V‘ v’0 H— ,' ’0 I— V---'V ,3 , '0-0 5 ,v” ‘0," ’v' (u —, .s 810 ‘, 7V 0" ” at ., FE g ,x" v IDCZ “63 ,x" O IDC3 '- _,0 x’ x IDC4 10 -4 . ...-...._3 . 1....-.-2 . 10 1O 10 At (a) lst order split (0, 1, 2, 3 corrections) 0 adv2D: split IDC, Nx=NV=4OO N o' __--..7 V II ---- ‘7' -5 V ————— V “I: 10 * v ‘ a x. B .—""‘ 510-10, v RK2 ”2°" *5 x IDC4 "’ ax” “ 1: IDCG ,n’ -4 -3 ' A A H A —2 10 1O 10 At (b) 2nd order split (0, 1, 2 corrections) Figure 5.1: Convergence study for (5.149) using (5.1a) IDC constructed with lst order split methods, each split equation is solved via forward Euler, combined with 5th order WENO. The order of accuracy is clear, as reference lines (with slopes of 1, 2, 3, 4) indicate. (5.1b) IDC constructed with 2nd order split methods, each split equation is solved via 2nd order Runge-Kutta, combined with 9th order WENO. The order of accuracy is clear, as reference lines (with slopes of 2, 4, 6) indicate. In both figures, the error is in the norm ll + ’00. compared to a reference solution computed on a finer time mesh. 179 Although we also used semi-Lagrangian split IDC methods with the conserva- tive WENO construction, the plots are not shown here because the error obtained was effectively machine precision. There is no splitting error for (5.149) since the operators commute, so those results are not unexpected. However, since there is no splitting error, the test of constant advection is not sufficient to verify split IDC methods’ increase in order of accuracy as the number of correction loops increases. Thus we consider the rotating problem next. 5.9.2 Rotation Using split IDC methods with conservative semi-Lagrangian 5th order WENO, we solved the problem of rotating a nonsymmetric Gaussian initial condition “t — vug; + xuv = 0, (x. v) E [—1, 1] X [—-1,1], t > 0, (5.150) u(0. x, v) = exp (—20(x2 + 5v2)), u(t,1,v) = u(t, —1, v). u(t,x,1) = u(t, x. —1). Here we used A = —v, B = x for our splittings. Time convergence results can be seen in Figure 5.2a for solving (5.150) via split IDC that incorporates a first order splitting in prediction and correction loops. We see that the temporal order of accuracy increases by 0(At) with each additional correction. Although the figure shows results when N3; = NU = 40, which may raise questions since the Gaussian is sharp enough that more spatial resolution is needed to capture its peak, we found identical results for N1; = ND 2 400. Figure 5.2b shows convergence for split IDC that incorporates second order splitting in prediction and correction loops with N1: = NU = 400. A prediction loop only (or simply a standard implementation of Strang splitting) gives second order. For a prediction and one 180 correction loop, an increase in accuracy by (2(At2) is expected, and indeed clean fourth order is acheived. For a prediction and two correction loops, sixth order is expected, but it is hard to verify due to acheiving Matlab precision quickly. However, the magnitude of its error is less than the fourth order error. Since the split operators used here for (5.150) do not commute, the splitting error will contribute to the error that the IDC methods must correct, and in light of our results, we conclude that IDC methods are able to correct splitting errors to achieve higher order split methods. Unfortunately, however, the IDC corrections appear to introduce a CF L that is not present in the prediction loop. The reason remains to be discovered. Now we move on to more interesting examples. 5.9.3 Drifting Rotation In this section, we consider a problem that is similar to the rotating problem, except there is a drift in time. We refer to the problem as a drifting rotating problem, but it can also be considered a linear Vlasov equation [69]: Htf+1’v$f+E(t'T)v'Uf:01 f(01miv) =f0(.77,U), (5151) where E (t, x) : [0, 00) x RN -+ RN is given and smooth. As in [69], we consider :1: and v E R with E(t,x) = —x — sin(2l), .7: : R —-> ’R. defined by .7:(0) = 0 and f’( ) 3&7: sin7(rrs) if 0 < s < 1 S = 0 otherwise 181 rot: SLw5, ichLsplit1, Nx=Nv=4O 0 10 10 V IDCZ II ..I: 0 IDC3 :7 ...—o -5 x .V"' p‘ (010 - ”304 W,»- ,5 8 37—“ ,'0 ,x L- V" ,6 x a) '6' ’I‘“ U ” g’ S -10 o x” a) 10 " " ‘ _4 ...._3 ...._2 10 1O 10 At (a) lst order split (1, 2, 3 corrections) _5 rot: SLw5, ichLstr, Nx=Nv=4OO 1o - . . ~ - V d ,,V'"v D ,’ II ”.V' x ’ ’ ‘0': V D ’1 .... ...-10 I ‘5 10 1 , " lDC4 l— X I e ’x , V str 5 ,x c1 IDCB o x’ . g X ’I’ I” a) _15 " X D 10 ...---LL’ 1:1 1 -3 -2 -1 10 1O 10 At (b) 2nd order split (0, 1, 2 corrections) Figure 5.2: Convergence study for (5.150). (5.2a) IDC constructed with lst order split methods, each split equation is solved in a semi-Lagrangian setting with con- servative 5th order WENO. The order of accuracy is clear, as reference lines (with slopes of 2, 3, 4) indicate. (5.2b) IDC constructed with 2nd order split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th or- der WENO. The order of accuracy is clear for 2nd and 4th order, but unclear for 6th order, as reference lines (with slopes of 2, 4, 6) indicate. In both figures, the error is in the norm ll + 100. comparing successive solutions. 182 so that if F(0, x, v) = f/(x cosl + v sinl — 2sin1+ sin 2) - f,(—x sin 1 + v cos 1 - 2cos1+ 2cos 2), then the exact solution [68] is f(t, x, v) 2: F’((x — sin(2t))cos(t—1)—(v — 2cos(2t)) sin(t — 1) + sin 2) -.7-'/((x — sin(2t)) sin(t — 1) — (v — 2cos(2t)) cos(t — 1) + 2cos2), which for t = 1 gives f(l, x, v) = .77’(x).7'"(v). Here we use A = v, B = E for our splittings. Our convergence results can be seen in Figures 5.3a, 5.3b, 5.4a, and 5.4b, where the error is found by comparing the numerical and exact solutions in Figures 5.3a and 5.3b and successive solutions in Figures 5.4a and 5.4b. We see that split IDC that incorporates a first order splitting in prediction and correction loops has the expected convergence, where the order of accuracy increases by 0(At) with each additional correction loop. As with the solution of the rotating problem, we observed a CFL limitation for split IDC methods. Thus a future step for us is to determine error bounds that clarified the CFL condition for split IDC methods. However, other higher order split methods also exhibit a CFL limitation. For example, Schaeffer tested the second order splitting method from [13] and his own fourth order splitting method [69] combined with a third order interpolation scheme on (5.151). He tested the 183 driftrot202: SLw5, idCSLspIit1, Nx=200, Nv=150 10° : v, ‘17 . r: E -1 §10r 5'1 *5 (0 X a) -2 1o_ - 1o3 102 101 At (a) 1 prediction, 0 corrections driftrotZDZI SLw5, idCSLspIit1, Nx=800, Nv=800 —2 10 : II : u— -3 :10 ‘ 0 (U h 8 : - a, 10’4. , |1+|°° error "('5 i o g —slope2 CD . -5 10 A A . . . -1 A m n . . .. -4 -3 -2 10 10 10 (b) 1 prediction, 1 correction l | At Figure 5.3: The order of accuracy is clear, as reference lines (with slopes of 1, 2) indicate. In all figures, the error is in the norm l1 +100. 184 driftrotZDZ: SLw5, idCSLspIit1, Nx=200, Nv=200 0 10 : : .- , l1+|°° error ll ~I_—_ —slope2 '26 o 5 o 5 . 8 . 3 (D -5 1O . 5* “Ar : . -4 -3 -2 10 10 10 At (a) 1 prediction, 1 correction driftrotZD3: SLw5, idCSLspIit1, Nx=800, Nv=800 -5 10 I +I' ~ .— . 1 ”error ll ‘ ..l: -—slope3 «a o L g . (D 0 l O 3 w l -10 10 _4 ....1-3 . .....-.-2 1O 10 1O ' | At (b) 1 prediction, 2 corrections Figure 5.4: The order of accuracy is clear, as reference lines (with slopes of 2, 3) indicate. In all figures, the error is in the norm l1 + loo- 185 4 methods numerically with Ax 2 Av and At = c(Ax)5 for a fixed constant c for his method and At = C(AIE)§ for a (possibly different) fixed constant c for [13]’s method. The values of At were chosen based on a CFL condition found through a rigorous error bound, and so that the error bounds in each case should be optimized [69]. Although Schaeffer’s fourth order method showed an improvement over Cheng and Knorr’s method, we note that Schaeffer’s method was only implemented for a linear problem. He displays no results for the nonlinear Vlasov system, and it is not clear how his method extends to the nonlinear problem. In the following section, we show that we may apply our split IDC methods to the nonlinear Vlasov—Poisson system. 5.9.4 Vlasov—Poisson Now we consider the implications of using split IDC methods with conservative semi- Lagrangian WENO to solve the Vlasov-Poisson system (5.88) with initial conditions of a two stream instability and Landau damping. In addition to temporal conver- gence results, we also include plots of solutions and physical quantities that should be conserved theoretically. Since the spatial resolution and methods are largely what influence the accurate numerical representation of the physical quantities, we do not expect that high temporal order via split IDC methods will improve them; however, we hope to assert that the physical properties are not unduly degraded by the split IDC methods’ error correction process. For both the two stream instability and Landau damping, we impose periodic boundary conditions in the x-direction and Neumann boundary conditions in the v-direction. Periodicity in space allows the use of a fast Fourier transform (F FT) to solve the 1D Poisson equation. The density p(x, t) is computed by the rectan- gular rule, p(x,t) = ff(x,v,t)dv z 29:1 f(x,vj,t)6v and then is normalized as 186 described in Section 5.6.3. For the splittings, we use A = v, B = E. Convergence studies are performed using N3; = NU = 400. Classical theoretical results for the Vlasov-Poisson system include conservation of the following quantities (i.e., the time derivative is zero): 0 The LP norm, for 1 S p < 00, f x.v,t pdxdv, 5.152 ./V /X| ( )l ( ) 1' l o Entropy, — f(x,v,t)ln(f(x,v,t) dx dv, 5.153 it. > < > a Total energy, 1//f(r1l)2(l(l+l/F(t)2d (5154) — .2, .1, , v x .v - 1 x, . x. . 2 v X 2 X In our numerical experiments, we checked the time evolution of the discrete versions of these theoretically preserved quantities. Entropy, energy, and the L1 and L2 norms are approximated by the rectangular rule. Two Stream Instability For the two stream instability, we use (5.88) with the initial condition [26] f(0. x, v) = %(1+ 502) (1 + a((cos(2kx) + cos(3kx))/1.2 + cos(kx))) - exp(—i,12/2). (5.155) 187 with a = 0.1, k = 0.5, x E [0,27r/k], v E [—7r/k,7r/k], periodic BCs in the x direction, and Neumann BCs in the 2) direction. Figures 5.5a and 5.5b show time convergence for split IDC methods with three and four, respectively, nodes in each subinterval, and in both figures, the error is shown for the prediction loop alone, followed by the prediction plus an increasing number of correction loops. As with the constant advection, rotation, and drifting rotation examples, we see a decreased error and very clear improved order of accuracy with each additional IDC correction loop. A first order split solution and a fourth order split IDC solution are shown in Figures 5.6a and 5.6b, respectively. In both plots, N1; = NU = 400, At = 1/300, and the final time is Tf = 25. The physical quantities are shown in Figures 5.7a, 5.7b, 5.7c, 5.7d, and 5.7e. In all plots, N]; = 64, NU = 128, At = 1/40, and M+l = 4. The energy, entropy, and L2 norms are virtually indistinguishable among the numerical methods, so split IDC methods conserve these quantities equally as well as the first order splitting method. Although the L1 norm (Figure 5.7a) has one peak for both the first order splitting (lower peak) and all the split IDC methods shown (higher peak), note that the peak is very small (< 0(10—5), though the figure shows only 3 0(10’4». We note that it should not be surprising that split IDC methods do not preserve positivity because the structure of IDC methods closely resembles (and when RK methods are used in the IDC construction, actually is identical to, see [17]) the structure of RK methods, which are known to be non-positivity preserving for orders greater than one without special assumptions on the type of problem solved or the size of the timestep used [7, 61]. The integral f f f (x,v,t)dx dv is clearly conserved (Figure 5.7b), so these split. IDC methods maintain conservation of mass, consistent with the proof in Section 5.8. 188 =0.1 succ error at Tf Figure 5.5: Convergence study for (5.88) with ICs as in (5.155). IDC constructed with lst order split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th order WENO. (5.5a) Here three IDC nodes are used, with 0, 1, and 2 correction loops. The order of accuracy is clear, as reference lines (with slopes of 1, 2, 3) indicate. (5.5b) Here four IDC nodes are used, with 0, 1, 2, and 3 correction loops. The order of accuracy is clear, as reference lines (with slopes of 1, 2, 3, 4) indicate. In both figures, the error is in the norm ll +100, comparing = 0.1 succ error at Tf 10- 10 5 r ”u. -----¢- ------ o ------ ....-- in." ' IDC4J0 "7 """ V V lDC4J1 ....... V ° l004.12 ”Ob—v ...... v ‘0 x IDC4J3 "7 """ data5 "0"." data6 —°"" u" data7 'D'flfl ”x“ dataB -15' y XJ'L';:: ‘ A A A A A L A 10'4 12:, 1O- 5 10 ...... ‘ ----r ------ «I ------ «Ir-~— r .--- ...... v ————— -v ...... V’ ...... v 10-10 V """" v «0 """ 0 4 o"‘0 """ Me " IDC3J0 QMMV VI003.11 - o 1O 15 . ‘ 1-4 + L . 1-3 IDC3J2 10 10 two stream; w5; IDC3, split1, 0,1,2 corr; Nx=Nv=400 At (a) 3 IDC nodes twostream; w5; IDC4 split1; 0,1,2,3 corr; Nx=Nv=400 successive solutions. (b) 4 IDC nodes 189 re- w5. |DCm4J0. Tf=25. Nx=Nv=400. A t=0.0033 w5, IDCm4J3. Tf=25, Nx=NV=400, A l=0.0033 (a) lst order (b) 4th order Figure 5.6: Solution for (5.88) with ICs as in (5.155). IDC constructed with lst order split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th order WENO, NI = NU = 400, At = 1/300. Landau Damping Landau damping is given by (5.88) with the initial condition [26] f(0, x, v) = (1+ (1 cos(kx)) exp(—0.5122) NR (5.156) with k = 1, x E [0, 27r], v E [-27r, 27r], periodic BCs in the x direction, and Neumann BCs in the 1) direction, a = 0.5 for strong Landau damping, and a = 0.01 for weak Landau damping. Figures 5.8a and 5.8b show time convergence for split IDC methods applied to weak and strong Landau damping problems, respectively, and in both figures, the error is shown for the prediction loop alone, followed by the prediction plus an increasing number of correction loops. As with the previous examples, we see a decreased error and very clear improved order of accuracy with each additional IDC correction loop. For weak damping, a first order split solution and a fourth order split IDC solution are shown in Figures 5.9a and 5.9b, respectively. In both plots, Nag = 190 twostream twostream 30.4655- 30.4655 ---2nd . ---2nd 30. 4655 _ - 3rd 30.4655 . _ 3rd """ 4th 13 30.4655 "“‘4th 3 30.4655 A 1’; ‘1 3 5 30.4655» l I 30-4655’ M 30.4655» \ 30.46550 50 1 00 30.46550 go 1 00 t 1 (a) L (b) fff(x,v,t)dx dv twostream twostream 14.6” 28' —1st 14.4’ "'an - - 3rd 14 2_ ..... 4th a 27.5 N ' 9 -' ‘E 14. LU 27‘ —1st --2nd 13.8‘ ‘- - 3rd 13 6 26 5 . . ‘ """ 4th . ' 0 20 40 t 60 80 100 ' 0 20 40 t 60 80 100 (c) L2 (d) Entropy twostream A 1 A ..... 4th I ‘o 20 40 t 60 so 100 (e) Energy Figure 5.7: Physical quantities for (5.88) with ICs as in (5.155). IDC constructed with lst order split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th order WENO, N1; = 64, N0 = 128, At = 1/40, and M + 1 = 4. The labels 1st, 2nd, 3rd, and 4th in the legends denote first order splitting (prediction only), 2nd order split IDC (prediction, 1 correction), 3rd order split IDC (prediction, 2 corrections), and 4th order split IDC (prediction, 3 corrections) methods, resp. 191 weakLandau; w5; idc4 split1; 0,1,,2 3 corr; Nx=Nv- 4—00 10‘5 ". o lllllllllll . . ......... I" ll . , ........ I: “I: x 1- IDC4J0 ‘ ..... V a 10 ‘°- V ........... .v 1111111 v IDC4J1 . § v ....... V lllllllllllll 0 IDC4J2 ‘3 .......... 0 x IDC4J3 8 .......... 0‘“ , 3 -15 ........ ~0” ...... m 10 0 """ 8 x P" 4 _3 2 10 10 10 At (a) weak Landau 1strgngLandau; w5; idc4 split1, 0,1,,2 3 corr; Nx= NV: 400 .__ . . ........... ." - IDC4J0 T.) v v IDC4J1 t ,_ NV v o IDC4J2 a 10-105 v .V A ‘ __o x IDC4J3_ o v, 0 _ O " a 3 10'15 . , x Io" 10‘3 10‘2 A t (b) strong Landau Figure 5.8: Convergence study for (5.88) with ICs as in (5.156). In Figure 5.8a, a = 0.01, and in Figure 5.8b, a = 0.5. IDC constructed with lst order split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th order WENO. Here four IDC nodes are used, with 0, 1, 2, and 3 correction loops. The order of accuracy is clear, as reference lines (with slopes of 1, 2, 3, 4) indicate. The error is in the norm 11 + loo, comparing successive solutions. 192 w5. lDCm4JO. Tf=16. Nx=NV=400. A t=0.0033 w5. lDCm4J3. Tf=16. Nx=NV=400. A t=0.0033 5 (a) lst order (b) 4th order Figure 5.9: Solution for (5.88) with ICs as in (5.156), a = 0.01. IDC constructed with lst order split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th order WENO, N3; = NU = 400. (5.9a) At = 1/300. (5%) At = 1/300. NU = 400 and At = 1/300, and the final time is Tf = 16. The physical quantities for weak damping are shown in Figures 5.10a, 5.10c, 5.10d, and 5.10e. We can see that the quantities are essentially conserved for all methods shown, and IDC methods may have a slight improvement over the first order splitting alone. Second order split IDC methods were not included in most graphs since there appear to be stability issues due to the choice of time step. The expected damping of the electric field is seen qualitatively in Figure 5.11. For strong Landau damping, a first order split solution and a fourth order split IDC solution are shown in Figures 5.12a (At = 1/80) and 5.12b (At = 1/80), respectively. In both plots, N3; = NU = 400, and the final time is Tf = 16. The physical quantities for strong damping are shown in Figures 5.13a, 5.13c, 5.13d, and 5.13e. We can see that the quantities are essentially conserved for all methods shown, and IDC methods produce nearly identical results to first order splitting alone. Second order split IDC methods were not included in most graphs since there appear to be stability issues due to the choice of time step. The expected damping of the electric field is seen here qualitatively in Figures 5.14. 193 weakLandau weakLandau 6.2832 6.2832; —1st —1St ---2nd 1" 3rd > _ .-.- 3rd 6.2832 ----- 4th 3 62832 ----- 4th 1- 'U "' 3‘; 628320 50 Too 6"28320 50 160 t (a) L1 (b) fff(x,v,t)dx dv weakLandau weakLandau 1.7725' 8.9153( -1st 13.9153~""3rd 1.7725» 3. ----- 4th S g 3.91531 . 1.772 W 5 8.9153,——/__/__./ 17725 ..... 4th . 4 89153 ' 1 . ' 0 5t0 100 ' 0 5t0 100 (c) L2 (d) Entropy weakLandau 3.1418' 3.1418; > E’ L“ —1st 3.1417» .-.-3,d 31416 . """ 4th. ' 0 5t0 100 (e) Energy Figure 5.10: Physical quantities for (5.88) with ICs as in (5.156), a = 0.01. IDC constructed with 1st order split methods, each split equation is solved in a semi- Lagrangian setting with conservative 5th order WENO, N1; = 64, NU = 128, and M + 1 = 4. The labels lst, 2nd, 3rd, and 4th in the legends denote first order splitting (prediction only, At = 1 / 40), 2nd order split IDC (prediction, 1 correction, At = 1/40), 3rd order split IDC (prediction, 2 corrections, At = 1/80), and 4th order split IDC (prediction, 3 corrections, At = 1/80) methods, resp. 194 ...—l. .....— weakLandau 0. —1st A-10' ---2nd 5 ----'3rd 0 "20' ----- 4th 3 5-30- 0 _l _40. '500 10 2'0 3'0 t Figure 5.11: Electric field for (5.88) with ICs as in (5.156) and a = 0.01. IDC constructed with 1st order split methods, each split equation is solved in a semi— Lagrangian setting with conservative 5th order WENO, NI = 64, N1; = 128, M + 1 = 4. For 1st order: At = 1/80, 2nd order: At = 1/80, 3rd order: At = 1/80, 4th order: At = 1/80. w5, |DCm4J0, Tf=16, Nx=Nv=400, A t=0.0125 w5, |DCm4J3. T1=16, Nx=NV=400. A i=0.0033 (a) lst order (b) 4th order Figure 5.12: Solution for (5.88) with ICs as in (5.156), a = 0.5. IDC constructed with lst order split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th order WENO, Nx = NU = 400, M + 1 = 4. (5.12a) At = 1/80. (5.12b) At = 1/300. 195 strongLandau strongLandau L1 6.2832' 6.2832' ---2nd 6.2832» 3rd 5 6.2832 ----3rd """ 4t“ g; m--4th 6.2832- 2 6.2832 628320 50 160 628320 5? 160 t (a) L1 (b)fff(x,v,t)dxdv strongLandau strongLandau 2. 95' —1st ’---3rd A/ 1.9» ----- 4th a 9» N 2 _J E —1st 1.8l LU 8.5- 3rd \ ----- 4th 1'7 . I A I 8 l . + r . 0 20 40 t 60 80 100 0 20 40 t 60 80 100 (c) L2 (d) Entropy strongLandau 3.6' >.3.55' F.” “J 3.5» --3rd ----- 4th 3'450 5;o 100 (e) Energy Figure 5.13: Physical quantities for (5.88) with ICs as in (5.156), a = 0.5. IDC constructed with 1st order split methods, each split equation is solved in a semi- Lagrangian setting with conservative 5th order WENO, N3; = 64, N0 = 128, and M + 1 = 4. The labels 1st, 2nd, 3rd, and 4th in the legends denote first order splitting (prediction only, At = 1/40), 2nd order split IDC (prediction, 1 correction, At = 1/40), 3rd order split IDC (prediction, 2 corrections, At = 1/80), and 4th order split IDC (prediction, 3 corrections, At = 1/80) methods, resp. 196 strongLandau 0 —1st ---2nd fi_10_ - "3rd “5 ----- 4th N :l §’—20- "300 10 210 3‘0 - t Figure 5.14: Electric field for (5.88) with ICs as in (5.156), a = 0.5. IDC constructed with lst order split methods, each split equation is solved in a semi-Lagrangian setting with conservative 5th order WENO, Nx = 64, N1) = 128, M + 1 = 4. For lst order: At = 1/80, 2nd order: At = 1/80, 3rd order: At = 1/80, 4th order: At = l / 80. 5. 10 Conclusions Here we have investigated high order splitting methods using IDC methods con- structed with low order splitting methods. We have found that these split IDC methods have potential for improved efficiency over other general order splitting methods, requiring a number of solves that increases linearly rather than expo— nentially with the order of the method. We also proved conservation of mass, and applied some split IDC methods to constant advection, rotating, and Vlasov-Poisson equations to verify the order of accuracy of split IDC methods. We found that split IDC methods maintain conservation properties at least as well as low order splitting methods. Ongoing work is to determine why a CFL restriction is introduced in the IDC correction step, whether that restriction may be mitigated, and to investigate high order splitting for other PDEs. The principles for split IDC methods could 197 be extended easily from the Vlasov—Poisson Operators to other operators suitable for splitting methods, such as the Vlasov—Maxwell, BGK—Poisson, or BGK—Maxwell equations. Future work should also consider how well the high order split IDC meth— ods handle Vlasov equations with certain collision operators (such as BGK), and how boundary conditions may be handled in such a way that high order accuracy is still maintained. 198 Chapter 6 Future Work 6.1 Asymptotic Preserving Methods As mentioned in Chapter 3, semi-implicit IDC methods seem ideal to investigate higher order asymptotic preserving (AP) methods by utilizing popular low order AP ARK methods. An AP method is a numerical scheme designed for specific models containing a small scale 6 (such as shallow water equations with strong diffusive terms), such that numerical solutions behave similarly to analytic solutions in the asymptotic limit, and the timestep is independent of e [43, 44, 64, 63]. Both the numerical method and the form of the model are interdependent and must be chosen carefully to achieve AP results. However, AP schemes are not yet higher than third order in time, and it is the existing first and second order schemes that work best thus far [64, 34]. Using IDC methods, we expect to be able to extend these methods to higher order in time. We anticipate that incorporating AP schemes into the IDC framework will make it possible to maintain stability and mitigate order reduction that occurs in many semi-implicit methods as e shrinks. Additionally, since AP ARK methods may also be described as splitting methods which solve the stiff part of the system first [64, 63], we expect some relationship between the split IDC methods 199 II of Chapter 5 and IDC methods that are constructed via AP ARK methods. In the subsequent sections, we outline details of our current plans for investigating whether IDC methods may possess AP properties. 6.2 IDC with AP-ARK Theory 6.2.1 Moving Towards an AP IDC Proof The algorithm for IDC methods constructed with AP ARK integrators is identical to the algorithm given in Chapter 3 for IDC—ARK methods, being careful to note that the CZN aé C? for AP ARK integrators. A main concern is whether the resulting IDC—ARK method will preserve the AP property of its constituent ARK method. Applying the ideas in [64, 63], we present a sketch that indicates that, under certain assumptions, and provided we may answer some questions, it is likely that we may obtain IDC-ARK methods that are AP methods. From [64, 63], we wish solve 1 0,0 = —(‘);I;F(U) + 2aw), (6.1) (6.2) where U 6 RD, F : RD x RD, and R : RD x RD. If R is a relaxation operator, then R(U) = 0 in the limit as e —> 0, and there exists a constant d x D matrix W with rank(W’) = d < D and WR(U) = 0 VU 6 RD. (6.3) Then also there are d conserved quantities u such that u = WU and R(U) = 0 can be solved uniquely for U as a function of u, that is, U = 8(a) and [2(5) (11.)) = 0. Then we may obtain a system of d conservation laws by applying W and (6.3) to 200 a,(WU) + 0;,;(WF(U)) = 0, (6.4) which is holds for every solution U of (6.1). Since 6 ——> 0 gives R(U) = 0, then the equilibrium system (9t(u) + Elli-flu) = 0, (6.5) where 6(21.) = W F (8 (11.)), is a good approximation to (6.1). 0 Prediction Consider the following prediction step to solve (6.1) via an AP ARK method for the predicted solution 771 at the next timestep y(t) = L70 H N (1 VS 81 H + Al 2 (lij(—5)IF(U J )lt=t0+CNAt) + Z aijERa] J )lt=t0+CSAt ’ 't=1,...,US, (6.6) ”1 = U0 + Al VN VS 1 N (2') S - (2') / _ F — 1 Z aI/‘Vi( 81: (U )ltzt0+cNAt) + Z Gust 6 RH] )|t=t0-I~C‘-S At 2:1 ’1. 1:1 I When 5. —r 0, S S 5.1 (1') _ -_ S Zlal/SiaULR(U )—0, 2—1,2,...,V . J: 201 From a lemma in [63], if AS (the matrix from the Butcher table coefficients for at?) is nonsingular, then R(U(J j—)) — 0. If the implicit part of the ARK method has stiff decay, i.e., if aVS]. = b‘JS, then 771 = U(VS ), and R071) = 0 also. By a theorem in [63], then the AP ARK scheme in (6.6) becomes the RK scheme characterized by the explicit part of the ARK scheme applied to the equilibrium limit system (6.5). The proof for the prediction step only requires that all R(U(j)) = 0, but the proof for the correction step likely requires that R(771) = 0 also. However, if R(771) = 0 also, then the AP ARK scheme is also asymptotically accurate for the prediction step; i.e., the order of accuracy of not only the d conserved quantities but also of the D — d nonconserved quantities should be preserved. 0 Correction Now we consider the error equation in order to analyze a correction loop. The error equation is given by (976 = 8,0 — 87,7] 1 1 = —81.F(U) + ERW) + 511F107) — 23(0) - where 'r is the residual, 1 r = 077; + 0:7F(77) - 23(7)) /tr(7')d’r— (t)—U +/t8 F(7)—1R( )dr .tO "‘77 O .to 113 7 C -77 . Then the error equation becomes ())(e+/r)=—()FF<77+1((+/7‘) —7"-:-/)+—R(77+(c+/7') —/) + 81, F( ) — —R() 202 or, letting Q = e + f r, . .. 1 . 1 (W = 4711‘ (77+ Q — fr) + 211(7) + Q — fr) + 0.217(7)) — 213(7)). 01' t 870 = 49,51? (Q + U0 —/ 8315(7)) —1R(n)dr) t0 ‘5 1 t 1 —R — aL-Fw ——R7d'r +. (ewe fto. <1) . <1) ) + 8:1:F(77)— %R(n), 0 Q0=co+/0 T=0. (6.7) Now we solve equation (6.7) via an AP ARK method whose implicit part has stiff decay, so Q(VS) is also the approximate solution to (6.7) at the next timestep. . VN . t0+cNAt 1 Q“) = Q0 + At]; 11,1] ($5.77 ((20) + U0 — [t0 9 613mm — ERMW) +2),F(17(70 + cyan) ) VS 31 . t0+CSAt 1 + At aij; R 62(3) + U0 —/ J 833F(77) — ;R(77)d’r . t —R(n(tO+cfAt)) ), i=1,...,I/S. 203 Then multiply by c to obtain . ”N N ( .) t0+CNAt 1 EQU) : EQO + At 2 atj —58;,;F Q J + U0 —/t J 3317(77) — ;R(77)d'r j=1 0 +6033F(77(t0 + cj-VAtD ) S s V . t0+c,- At + At 2 655']. R (.20) + U0 —/ J 633F(77) — 17(7))dT j=1 t0 —R(n(t0 + c3980) ) , (6.8) and consider the limit as c approaches zero. The left hand side of (6.8) goes to zero. The first term on the right hand side also goes to zero. We wish to prove that the entire explicit sum and the second part of the implicit sum go to zero, which would leave us with (j) 10+CJS At . S R Q + U0 — t 81,1707) — R(77)d7' = 0, ] = 1,. . . ,I/ , (6.9) 0 since AS is nonsingular. If (6.9) holds, then, by a similar process as in the rest of the proof of the theorem in [63], the correction step possesses the AP property. It seems likely that the explicit sum goes to zero, but it is important to take care with the fact that the quantity inside F in the first part of the explicit sum contains a 1/ 6. term. It also seems likely that the second part of the implicit sum goes to zero, but the fact that the values in the second part of the implicit sum are interpolated values may affect that conclusion. However, we know that R(77(t0)) goes to zero, and if the AP ARK scheme has stiff decay, then we know that R(771) go to zero as well. By the same reasoning as for the prediction, the correction step is also asymptotically accurate if the AP ARK method’s implicit part has stiff decay. 204 6.3 AP IDC Methods for P1 Equations Here we outline a plan for understanding how to implemen AP IC methods applied to the P1 equations in [34]. 6.3.1 The P1 Problem and the Splitting We begin iwth a presentation of the P1 system and its splitting that is used in AP methods. The original P1 system that we desire to solve is (97p + aaxm = 0, (6.10) b a (97771 + :26” = -—€—2m, 10(13 to) = 100(1), m(x, t0) = m0(x), where a, b are constants, 0 = 0(x), and t E (to, to + At]. The first order splitting of (6.10) is given by (97,00) 2 0, (6.11) atm(1)+—b§(1— (2)3320“) = —%m(1), e c p<1> 0, the system (6.10) becomes b 207 For an asymptotic preserving method, the numerical solution of (6.10) should match the numerical solution of (6.17) in the C ——+ 0 limit, for At > 62. 6.3.3 Splitting Error Now we analyze the splitting error by Taylor expanding the exact solution 11 of (6.14), A12 u.(;1.',t0 + At) = (I + Am, + Ti)“ + . . . )'u,0(x), where I is the identity matrix. We replace the t derivatives with their equivalent expressions from (6.14), 1 8 — 1F G 1F G 11 - ‘65 + 3 + 1 = —4F2+—12-(FG+GF)+GQ, e c and obtain 1 u(x,t0 + At) = (I + At (7F + C) (6.18) 6 At? 1 1 +— —F2+—-(FG+GF)+02 +... 110(5). 2 54 £2 208 Similarly, we expand the solutions to (6.15) and (6.16) to obtain 1 At21 u(1)(x,t0+At)—-—(I +At 2F+— 2 4F2+...) u(1)(x,t0), (6.19) 1 At21 I + At— 217+ —2-—F2+. ) 110(x), {4 u(2)(x, (.0 + At) 2 1 + AtG + 91—024mm) u(2)(x, 10) < ( 2 ( At < I+AtG+ —02+ ...)u(1)(x,t0+At) +A__t_2 2+ 2 , ,2 Since F and G do not commute, then the exact (6.18) and split (6.19) solutions I+At( differ in the At2 term, giving a first order splitting. Alternatively, the solution to the split equations is given by u(2)(x, to + At) = exp ((to + At)ti2F) exp ((tO + At)G) u0(x), which indicates that, if F is negative semi-definite, then it seems more likely that the method will approximate (6.17) as c —-> 0. Also, for the split solution not to blow up as c —> 0, we require 6.3.4 Split Scheme with one IDC Correction Now we present the form that the error equation and the split error equation should take in one correction loop of an IDC method. Let u (as defined in (6.13)) be the exact solution to (6.14), 7) be the split solution after solving (6.15) and (6.16) exactly, 209 i.e. 77 = 11(2), and define the error, e, and residual, 7‘, as €=u—77, 1 7:8777— (7F+G) 77. c Differentiating 6 with respect to t, we obtain the error equation, ate = 87,11 — 6,77 (6.20) = (—1§F+G)u—((—12—F+G)77+r) e e 1 = (—2F+G)e—r, 6 e(x,t0) = 0. Then we approximate the solution to the error equation (6.20) by the three-part first order splitting ate“) = 32-186“), (6.21) 6(1)($, to) = 60(4); (976(2) 2 06(2), (6.22) 6(2)(x, t0) = 6(1)(.’E,t0 + At); 876(3) = —7‘ 2 "‘8t (77 - 110 — ft; (ciQF + G) ndr) , (6.23) 6(3)(x, t0) = 8(2)(.’L‘, to + At). 210 Now for comparison with the common form of (6.10), we change from the vector notation back to the component-wise notation. The error equation (6.20) becomes t (976,0 + adzem = —Bt (77p — p0 +./t (183377771 d7) , 0 b _0 t b 1 — 62 0 876m. + —8;L-ep- — —em- 8,: 77m — mg + (b + —(——))8;L-77p + —77m d7 , 62 £2 to £2 62 ep(x, t0) = ep’0(x), em(x, to) = 8771,0(5’3lv and the split form (6.21), (6.22), (6.23) of the error equation becomes 6,421) = 0, (6.24) 64.? 6326* 185-J )— ”2 5.1.). 2% to) = e 0(4). 4336 to) — em,6< ) 6,4,2) + 6638),?) = 0 (6.25) 876m (2 l + 56,8(2)=0 0, (6.26) 6)?)(A10) = 621%,10 + At), (6.27) (4,2,)(5, t0) = 6);)(5, 10 + At); (6.28) 211 3 t ate), ) = —8L (77p — p0 +/ aaxnm d7) , (6.29) 3 t b 1 — c2 ‘ 0 atein) = ‘at (77m - m0 +/t (5 + i17—2Wx77p + —2'77m dT) , (6-30) 6 0 e£3)(x,10) = e£2)(1‘,t0 + At), (6.31) 12,2)(1‘, t0) = 65.3)(610 + At) (6.32) 6.3.5 Analytic Solutions of Splitting Parts Here we give the analytic solutions to the split equations (6.11), (6.12), and the split error equations (6.24), (6.25), and (6.29). To simplify the presentation, let t = to + At. The solution to (6.11) is exactly 6%) = Mao), (6.33) mum) = _exp (22.)) [5) 32‘3"" (£2) awn) d. where the second equality for ma) holds if there is enough smoothness such that p“) constant in t gives Egg/2(1) also constant in t. Now we solve (6.12). Left multiplying (6.12) by 1 1 95 2 b VL: 1 ia— ’ ‘95 NEE one obtains 0,111+ + Marv = 0, (6.34) 0,11" — Hagar = 0, <2) (2) (2) (2) + _ p m _ 2 _p m ‘1’ ‘2a+2\/25’ q’ 2113.471)“ Then the solution for (6.34) is given by 111+(1;,1) = 111+(1— M160), \I'—(:1:,t)= \II—(x + \/Zl_bt,t0). Then some simple calculations give the solution to (6.12) as p(2)(.1 1)=,%(() ()(2 (1: — V7110) +p(2)($+ $610)) (6.35) + —\/\7__- ("1(2) (:1: — at, 1.0) — 7n,(2)(:17 + \/(Tbl,, 1.0)) , 2 , x/5 111.( )(1.1)= 375(p(2)( (.1—x/a—bt, 10)— [)(2)(:I:+\/a—bt,t0)) +%(m( 2)(a:—\/.(16t,t0)+m(2)(:1:+\/o—5t.t0)). The solutions for the split error equations (6.24) and (6.25) are analogous to the solutions (6.33) and (6.35) of (6.11) and (6.12), respectively, and the solution for the part of the error splitting (6.29) is given by 1 15.311. 1) = «22%,10) — (111111») — 120(1) + / 1161111116, 1) d1) , (6.36) t0 12911.1.) = 125%, to) t — 62 a — (7777203, t) — m0(x) +/t (b + M—lcz—lwxanJ) + 277m(a:,7) dT) . 6.3.6 Anticipated Numerical Tests Here we discuss the plans for numerical tests to determine the effectiveness of IDC in the AP setting, in particular, to determine whether the IDC framework is capable of maintaining the AP property. To start with, we plan to use periodic boundary con- ditions and smooth initial conditions. We will solve each splitting analytically, which is essentially a semi-Lagrangian method. It may be necessary to use WENO recon- struction for some places where 83; appears, but in other places, setting Ax 2 At may mean WENO reconstruction is unnecessary. Note that within the WENO method, the Lax-Friedrichs flux may cause problems for these P1 equations. Per- haps instead a cubic spline may be used, especially if the nonoscillatory features of WENO are not required. We intend to test various values of c, for At > 62. Note, however, that the CF L limitation found numerically for the split IDC meth- ods in Chapter 5.6 will need to be understood much better before these tests can be performed reasonably. 214 [ll [2] [3] [4] [5] l6] [7] [8] [9] BIBLIOGRAPHY A. L. Araujo, A. Murua, and J. M. Sanz—Serna. Symplectic methods based on decompositions. SIAM J. Numer. Anal., 34(5):1926—1947, 1997. Uri M. Ascher and Linda R. Petzold. Computer methods for ordinary difler- ential equations and diflerential—algebraic equations. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1998. Uri M. Ascher, Steven J. Ruuth, and Raymond J. Spiteri. Implicit-explicit Runge—Kutta methods for time—dependent partial differential equations. Appl. Numer. Math., 25(2—3):151—167, 1997. Special issue on time integration (Ams- terdam, 1996). Uri M. Ascher, Steven J. Ruuth, and Brian T. R. Wetton. Implicit—explicit methods for time-dependent partial differential equations. SIAM J. Numer. Anal., 32(3):797—823, 1995. W. Auzinger, H. Hofstatter, W. Kreuzer, and E. Weinmiiller. Modified de- fect correction algorithms for ODEs. I. General theory. Numer. Algorithms, 36(2):135—155, 2004. W. Auzinger, H. Hofstatter, W. Kreuzer, and E. Weinmiiller. Modified de- fect correction algorithms for ODEs. II. Stiff initial value problems. Numer. Algorithms, 40(3):285—303, 2005. C. Bolley and M. Crouzeix. Conservation de la positivité lors de la discrétisation des problemes dévolution paraboliques. RAIRO Anal. Nume’r, 12(3):237—245, 1978. Sebastiano Boscarino. On an accurate third order implicit-explicit Runge— Kutta method for stiff problems. Appl. Numer. Math., 59(7):1515—1528, 2009. Anne Bourlioux, Anita T. Layton, and Michael L. Minion. High-order multi- implicit spectral deferred correction methods for problems of reactive flow. J. Comput. Phys, 189(2):651—675, 2003. 215 [10] T. J. M. Boyd and J. J. Sanderson. The physics of plasmas. Cambridge Uni- versity Press, Cambridge, 2003. [11] M. P. Calvo, J. de Frutos, and J. Novo. Linearly implicit Runge—Kutta methods for advection-reaction-diffusion equations. Appl. Numer. Math., 37(4):535—549, 2001. [12] J. A. Carrillo and F. Vecil. Nonoscillatory interpolation methods applied to Vlasov-based models. SIAM J. Sci. Comput., 29(3):1179—1206 (electronic), 2007. [13] C.Z. Cheng and G. Knorr. The integration of the Vlasov equation in configu- ration space. J. Comput. Phys, 22(3):330—351, 1976. [14] A. Christlieb, R. Krasny, and B. Ong. Particle simulation of Weak Solutions of the Vlasov—Poisson equations. In preparation. [15] A. Christlieb, M. Morton, B. Ong, and J.-M. Qiu. Semi-implicit integral de- ferred correction constructed with high order additive Runge—Kutta methods. Submitted. [16] Andrew Christlieb, Colin Macdonald, and Benjamin Ong. Parallel high-order integrators. SIAM J. Sci. Comput, 32(2):818—835, 2010. [17] Andrew Christlieb, Benjamin Ong, and Jing-Mei Qiu. Comments on high or- der integrators embedded within integral deferred correction methods. Comm. Appl. Math. Comput. Sci, 4(1):27—56, 2009. [18] Andrew Christlieb, Benjamin Ong, and Jing—Mei Qiu. Integral deferred cor- rection methods constructed with high order runge—kutta integrators. Math. Comput., 79:761—783, 2010. [19] B. Cockburn, C. Johnson, C.-W. Shu, and E. Tadmor. Advanced numerical approximation of nonlinear hyperbolic equations. l697:vi+433, 1998. Papers from the C.I.M.E. Summer School held in Cetraro, June 23—28, 1997, Edited by Alfio Quarteroni, Fondazione C.I.M.E.. [C.I.M.E. Foundation]. [20] G. J. Cooper and A. Sayfy. Additive methods for the numerical solution of ordinary differential equations. Math. Comp, 35(152):1159—1172, 1980. [21] G. J. Cooper and A. Sayfy. Additive Rungo—Kutta methods for stiff ordinary differential equations. Math. Comp, 40(161):207——218, 1983. [22] G.-H. Cottet and P.-A. Raviart. Particle methods for the one-dimensional Vlasov-Poisson equations. SIAM J. Numer. Anal, 21(1):52—76, 1984. [23] N. Crouseilles, M. Mehrenberger, and E. Sonnendriicker. Conservative semi- Lagrangian schemes for Vlasov equations. J. Comput. Phys, 229(6):1927 — 1953, 2010. 216 [24] P. Csomos and I. Farago. Error analysis of the numerical solution of split differential equations. Math. Comput. Modelling, 48(7—8):1090—1106, 2008. [25] Alok Dutt, Leslie Greengard, and Vladimir Rokhlin. Spectral deferred correc- tion methods for ordinary differential equations. BIT, 40(2):241—266, 2000. [26] F. Filbet and E. Sonnendriicker. Comparison of Eulerian Vlasov solvers. Com- put. Phys. Comm., 150(3):247—266, 2003. [27] R. Frank and C. W. Ueberhuber. Iterated defect correction for differential equations. I. Theoretical results. Computing, 20(3):207—228, 1978. [28] Reinhard Frank and Christoph W. Ueberhuber. Iterated defect correction for the efficient solution of stiff systems of ordinary differential equations. Nordisk Tidskr. Informationsbehandling (BIT), 17(2):146—159, 1977. [29] D. Gottlieb. Strang—type difference schemes for multidimensional problems. SIAM J. Numer. Anal, 9:650—661, 1972. [30] Sigal Gottlieb, David I. Ketcheson, and Chi-Wang Shu. High order strong stability preserving time discretizations. J. Sci. Comput., 38(3):251—289, 2009. [31] Sigal Gottlieb, Chi-Wang Shu, and Eitan Tadmor. Strong stability-preserving high-order time discretization methods. SIAM Rea, 43(1):89—112 (electronic), 2001. [32] L. Greengard. Spectral integration and two—point boundary value problems. SIAM J. Numer. Anal, 28(4):1071—1080, 1991. [33] Kjell Gustafsson. Control-theoretic techniques for stepsize selection in implicit Runge—Kutta methods. ACM Trans. Math. Software, 20(4):496—517, 1994. [34] JR. Haack and CD. Hauck. Oscillatory behavior of asymptotic-preserving splitting methods for a linear model of diffusive relaxation. Los Alamos Report LA-UR-08-0571. [35] Wolfgang Hackbusch. Bemerkungen zur iterierten Defektkorrektur und zu ihrer Kombination mit Mehrgitterverfahren. Rev. Roumaine Math. Pures Appl., 26(10):1319—1329, 1981. [36] E. Hairer, S. P. Norsctt, and G. Wanner. Solving ordinary diflercntial equations. 1, volume 8 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1987. Nonstiff problems. [37] E. Hairer, S. P. Norsett, and G. Wanner. Solving ordinary diflerential equations. I, volume 8 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 1993. Nonstiff problems. 217 [38] Ernst Hairer, Christian Lubich, and Gerhard Wanner. Geometric numeri- cal integration, volume 31 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 2006. Structure-preserving algorithms for ordinary differential equations. [39] E. Hansen and A. Ostermann. High order splitting methods for analytic semi- groups exist. BIT Numerical Mathematics, 49(3):527—542, 2009. [40] P. W. Hemker. Mixed defect correction iteration for the solution of a singu- lar perturbation problem. In Defect correction methods (Oberwolfach, 1983), volume 5 of Comput. Suppl., pages 123—145. Springer, Vienna, 1984. [41] Jingfang Huang, Jun Jia, and Michael Minion. Accelerating the convergence of spectral deferred correction methods. J. Comput. Phys, 214(2):633—656, 2006. [42] Jingfang Huang, Jun Jia, and Michael Minion. Arbitrary order Krylov de— ferred correction methods for differential algebraic equations. J. Comput. Phys, 221(2):739—760, 2007. [43] Shi Jin. Efficient asymptotic-preserving (AP) schemes for some multiscale ki- netic equations. SIAM J. Sci. Comput., 21(2):441—454 (electronic), 1999. [44] Shi Jin, Lorenzo Pareschi, and Giuseppe Toscani. Uniformly accurate diffusive relaxation schemes for multiscale transport equations. SIAM J. Numer. Anal, 38(3):913—936 (electronic), 2000. [45] Christopher A. Kennedy and Mark H. Carpenter. Additive Runge—Kutta schemes for convection-diffusion-reaction equations. Appl. Numer. Math., 44(1- 2)2139—181, 2003. [46] Wendy Kress and Bertil Gustafsson. Deferred correction methods for initial boundary value problems. In Proceedings of the Fifth International Conference on Spectral and High Order Methods (ICOSAHOM-OI) (Uppsala), volume 17, pages 241—251, 2002. [47] Anita T. Layton. On the choice of correctors for semi-implicit Picard deferred correction methods. Appl. Numer. Math., 58(6):845—858, 2008. [48] Anita T. Layton and Michael L. Minion. Implications of the choice of predictors for semi-implicit Picard integral deferred correction methods. Commun. Appl. Math. Comput. Sci, 221—34 (electronic), 2007. [49] Jongwoo Lee and Bengt Fornberg. A split step approach for the 3-D Maxwell’s equations. J. Comput. Appl. Math., 158(2):485—505, 2003. [50] Jongwoo Lee and Bengt Fornberg. Some unconditionally stable time stepping methods for the 3D Maxwell’s equations. J. Comput. Appl. Math., 166(2):497— 523, 2004. 218 [51] Randall J. LeVeque. Numerical methods for conservation laws. Lectures in Mathematics ETH Ziirich. Birkhauser Verlag, Basel, second edition, 1992. [52] R.J. LeVeque. Intermediate boundary conditions for LOD, ADI and approxi- mate factorization methods. 1985. [53] R.J. LeVeque. Intermediate boundary conditions for time-split methods applied to hyperbolic partial differential equations. Math. Comp, 47(175):37—54, 1986. [54] M.A. Lieberman and A.J. Lichtenberg. Principles of plasma discharges and materials processing. Wiley-Blackwell, 2005. [55] Hongyu Liu. Personal correspondence, 2008. [56] Hongyu Liu and Jun Zou. Some new additive Runge—Kutta methods and their applications. J. Comput. Appl. Math., 190(1—2):74—98, 2006. [57] Yuan Liu, Chi-Wang Shu, and Mengping Zhang. Strong stability preserving property of the deferred correction time discretization. J. Comput. Math., 26(5):633—656, 2008. [58] G. I. Marcuk. Some application of splitting-up methods to the solution of mathematical physics problems. Apl. Mat., 13:103—132, 1968. [59] Francesca Mazzia and Cecilia Magherini. Test set for initial value problem solvers, release 2.4. Technical Report 4, Department of Mathematics, University of Bari, Italy, February 2008. Available at http://pitagora.dm.uniba. it/ ~testset. [60] Michael L. Minion. Semi-implicit spectral deferred correction methods for or- dinary differential equations. Commun. Math. Sci, 1(3):471-—500, 2003. [61] A. Ostermann and M. van Daele. Positivity of exponential Runge—Kutta meth- ods. BIT Numerical Mathematics, 47(2):419—426, 2007. [62] Lorenzo Pareschi and Giovanni Russo. Implicit-explicit Runge—Kutta schemes for stiff systems of differential equations. In Recent trends in numerical analy- sis, volume 3 of Adv. Theory Comput. Math., pages 269—288. Nova Sci. Publ., Huntington, NY, 2001. [63] Lorenzo Pareschi and Giovanni Russo. High order asymptotically strong- stability—preserving methods for hyperbolic systems with stiff relaxation. In Hyperbolic problems: theory, numerics, applications, pages 241—251. Springer, Berlin, 2003. [64] Lorenzo Pareschi and Giovanni Russo. Implicit-Explicit Runge—Kutta schemes and applications to hyperbolic systems with relaxation. J. Sci. Comput., 25(1- 2)129—155,2005. 219 [65] Victor Pereyra. Iterated deferred corrections for nonlinear operator equations. Numer. Math., 10:316-323, 1967. [66] J.-M. Qiu. Lecture, IPAM, UCLA, 2009. [67] J.-M. Qiu and A. Christlieb. A conservative high order semi-Lagrangian WENO method for the Vlasov equation. J. Comput. Phys, 229(4):1130—1149, 2010. [68] J. Schaeffer. Personal correspondence, 2010. [69] Jack Schaeffer. Higher order time splitting for the linear Vlasov equation. SIAM J. Numer. Anal, 47(3):2203—2223, 2009. [70] Chi-Wang Shu. Essentially non—oscillatory and weighted essentially non- oscillatory schemes for hyperbolic conservation laws. In Advanced numerical approximation of nonlinear hyperbolic equations (Cetraro, 1997), volume 1697 of Lecture Notes in Math., pages 325—432. Springer, Berlin, 1998. [71] E. Sonnendriicker, J. Roche, P. Bertrand, and A. Ghizzo. The semi-Lagrangian method for the numerical resolution of the Vlasov equation. J. Comput. Phys, 149(2):201—220, 1999. [72] Hans J. Stetter. The defect correction principle and discretization methods. Numer. Math., 29(4):425—443, 1977/78. [73] Gilbert Strang. On the construction and comparison of difference schemes. SIAM J. Numer. Anal, 5:506—517, 1968. [74] J. van Dijk, G.M.W. Kroesen, and A. Bogaerts. Plasma modelling and numer- ical simulation. J. Phys. D: Appl. Phys, 42(19). [75] JP. Verboncoeur. Particle simulation of plasmas: review and advances. Plasma physics and controlled fusion, 47:A231—A260, 2005. [76] Yinhua Xia, Yan Xu, and Chi-Wang Shu. Efficient time discretization for local discontinuous Galerkin methods. Discrete Contin. Dyn. Syst. Ser. B, 8(3):677— 693, 2007. [77] H. Yoshida. Construction of higher order symplectic integrators. Physics Letters A, 150(5-7):262—268, 1990. [78] Pedro E. Zadunaisky. On the estimation of errors propagated in the numeri- cal integration of ordinary differential equations. Numer. Math., 27(1):21—39, 1976/77. 220