HHigh-Order Multiderivative IMEX Schemes
Alexander J. Dittmann ∗ Department of Astronomy, University of Maryland, College Park, MD 20742-2421
Abstract
Recently, a 4th-order asymptotic preserving multiderivative implicit-explicit(IMEX) scheme was developed [1]. This scheme is based on a 4th-order Her-mite interpolation in time, and uses an approach based on operator splittingthat converges to the underlying quadrature if iterated sufficiently. Hermiteschemes have been used in astrophysics for decades, particularly for N-bodycalculations, but not in a form suitable for solving stiff equations. In thiswork, we extend the scheme presented in [1] to higher orders. Such high-orderschemes offer advantages when one aims to find high-precision solutions tosystems of differential equations containing stiff terms, which occur through-out the physical sciences. We begin by deriving Hermite schemes of arbitraryorder and discussing the stability of these formulas. Afterwards, we demon-strate how the method of [1] generalises in a straightforward manner to anyof these schemes, and prove convergence properties of the resulting IMEXschemes. We then present results for methods ranging from 6th to 12th orderand explore a selection of test problems, including both linear and nonlin-ear ordinary differential equations and Burgers’ equation. To our knowledgethis is also the first time that Hermite time-stepping methods have been ap-plied to partial differential equations. We then discuss some benefits of theseschemes, such as their potential for parallelism and low memory usage, aswell as limitations and potential drawbacks.
Keywords:
High-order accuracy, Multiderivative, IMEX scheme, Hermiteinterpolation ∗ Corresponding author
Email address: [email protected] (Alexander J. Dittmann)
Preprint submitted to Journal of Computational Physics August 12, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] A ug . Introduction Numerous physical systems occur in nature that are difficult to modelnumerically because of stiff terms in the differential equations that governtheir evolution. For example, when simulating the evolution of fluid flowsthat include chemical reactions, cooling, or viscosity, the time step requiredfor stability can be reduced by orders of magnitude when using explicit meth-ods alone [e.g. 2, 3, 4]. Stiff differential equations also appear in combustionmodels [5], and circuitry [6], for example. If one wishes to efficiently carry outhigh-accuracy simulations of such systems, especially in multiple dimensions,high-order numerical methods are required. This work explores high-orderschemes constructed using Hermite interpolation in time, adapting them tosystems where one may use operator splitting to separate stiff and non-stiffterms.Hermite schemes were introduced in celestial mechanics as an alterna-tive to Adams-Bashforth-Moulton predictor-corrector schemes [7] and can bethought of as a multistep generalisation of earlier schemes based on Taylorseries [8]. As the name implies, Hermite schemes are based on constructinga Hermite interpolation polynomial in time. Thus, it follows that by calcu-lating up to the ( q − th -order derivative of a function f at r points in time,one can construct a Hermite interpolation polynomial of order ( qr − f when compared to Adams-Bashforth-Moulton schemes.Various high-order time-stepping algorithms exist for solving differentialequations. Runge-Kutta methods of orders higher than four typically involveeither very complex Butcher tableaus or tens of stages [9, 10]. Another familyof schemes, based on deferred corrections in a Picard integral formalism,requires large numbers of stages and multiple iterations to reach high orders[11, 12]. Thus, high-order algorithms using values only at the beginning andend of an interval offer potential improvements.The schemes in this work are largely based on previous investigationsfocused on N-body algorithms [13, 14, 15], and readers are encouraged toreference those works for further details concerning deriving the schemespresented herein. These schemes are usually applied using an explicit Taylorseries predictor, utilising derivatives of the interpolation polynomial at theend of the previous step, followed by a correction step (see, for example, theappendix of [15]). In general, this style of application is not suitable for solv-2ng stiff equations. However, Hermite quadrature rules can also be applied ina fully implicit manner, as well as the operator-split manner described in [1].To elucidate this, we follow [1] and consider the coupled set of differentialequations y (cid:48) ( t ) = z ( t ) , z (cid:48) ( t ) = g ( y ( t ) , z ( t )) (cid:15) , ≤ t ≤ T, (1)where g : R → R is a smooth function and we choose 0 < (cid:15) <
1. Theseequations are supplemented by initial conditions at t = 0 y (0) = y , z (0) = z . (2)The stiff relaxation parameter (cid:15) causes Equation (1) to be a singularly per-turbed equation, which exhibits multiscale structure that is challenging formany explicit methods. These problems are well documented in the liter-ature, and a number of books discuss the analysis and application of theseequations [5, 16, 17].To develop more general schemes, we follow [1] and also consider thefollowing generic system of differential equations with an additive right-handside dudt = φ ( u ) = φ E ( u ) + φ I ( u ) , (3)where φ E and φ I contain the non-stiff and stiff terms respectively.When problems can be split in the manner of Equation (3), the stiffcomponents may be treated implicitly and non-stiff parts explicitly, greatlysimplifying the implicit part of the calculation. This procedure leads toimplicit-explicit (IMEX) methods. Popular families include those based onRunge-Kutta and multistep methods. Recently, [1] introduced a multideriva-tive IMEX method, which effectively implements a Hermite scheme for stiffproblems. Because of its connection to Hermite schemes, it is straightfor-ward to generalise the method presented in [1] to higher orders by drawingon previous advances in N-body methods [15, 14]. Developing these high-order methods is the primary goal of this paper, but before doing so it isnecessary to introduce Hermite quadrature schemes of arbitrary order. Af-terwards, we show how these can be used to generalise the multiderivativeIMEX scheme of [1]. We provide a few theoretical results on their conver-gence, and apply these Hermite IMEX schemes to example stiff ordinary andpartial differential equations. 3 . Hermite Schemes Instead of repeating the entirety of derivations found elsewhere, we statea number of results [13, 14] and proceed to discussing stability properties.Let us consider the equation dudt = f ( u, t ) (4)and a time step ∆ t , defining h ≡ ∆ t/
2. Then, let us define a fitting polyno-mial at the centre of the interval [ t − h, t + h ] using the notation f ( n ) ≡ d n dt n fF n = hn ! f ( n ) ( u ( t ) , t ) , (5)which will be determined by linear combinations of values at the beginningend of the interval F ± n = 12 h n n ! (cid:0) f ( n ) ( u ( t + h ) , t + h ) ± f ( n ) ( u ( t − h ) , t − h ) (cid:1) . (6)We can then calculate an update to u , ∆ u ≡ u ( t + h ) − u ( t − h ), as∆ u = (cid:90) h − h f ( t ) dt ≈ (cid:18) F + 13 F + ... + 12 p + 1 F p (cid:19) ∆ t, (7)where the series is truncated at n = 2 p and the order of the scheme is 2( p +1).Then, we can solve for F n in terms of linear combinations of F ± n using F +0 F − F +2 F − ... = ( ) ( ) ( ) ( ) · · · ) ( ) ( ) · · · ) ( ) ( ) · · · ) ( ) · · · ... ... ... ... . . . F F F F ... = ¯ A F F F F ... (8)where ¯ A is a matrix constructed from the even columns of an upper triangularPascal matrix, and parentheses indicate binomial coefficients. Thus, for a2( p + 1)th-order scheme with p + 1 terms,∆ u = (cid:0) c p F +0 + c p F − + ... + c pp F ± p (cid:1) ∆ t, (9)4here c pk is defined for 0 ≤ k ≤ p as c pk = 1( − k (2 k )!!(2 k + 1)!! (cid:0) p − k + mp − k (cid:1) (2 k + 1)!!(2 k + 1 − m )!! (2 p + 1 − m )!!(2 p + 1)!! , (10)where m ≡ (cid:98) ( k + 1) / (cid:99) [13]. As an example, let us consider the resulting4th-order scheme ∆ u = ∆ t f + f ) − ∆ t (cid:16) f (1)1 − f (1)0 (cid:17) (11)where the subscript 0 indicates a value at the beginning of a time step andthe subscript 1 indicates a value at the end of a time step. Explicit formsup to 12th order are provided in Appendix A. Typically these schemes areapplied in a predictor-corrector fashion using an explicit Taylor series predic-tor, constructed using derivatives of the interpolating polynomial evaluatedat the end or the previous interval. However, this procedure is not suitablefor stiff equations, and we will instead follow [1], using a mix of lower-orderforward and backward Taylor series as a predictor. But first, we consider thestability of the underlying Hermite schemes.Let us apply Equation (11) to the test function dudt = f ( u ) = λu, f ( n ) ( u ) = λ n u : ∆ u = u − u = ∆ t λu + λu ) − ∆ t (cid:0) λ u − λ u (cid:1) . (12)Rearranging, and defining Φ( λ ∆ t ) ≡ u /u we see thatΦ( K ∆ t ) = 1 + λ ∆ t + λ ∆ t − λ ∆ t + λ ∆ t . (13)Analogously to the trapezoidal rule, we see that for any λ < t > | Φ( λ ∆ t ) | <
1. Thus, the 4th-order method is A-stable, with a linear stabilityregion consisting of the left half of the complex plane. Now we considerΦ( λ ∆ t ) for higher-order methods. Based on the sign of the coefficients c pk inEquation (10), we can see that any odd power of ∆ t will be accompanied bya negative coefficient in the denominator of Φ( λ ∆ t ) but a positive coefficientin the numerator. Thus, all Hermite schemes are A-stable. In this section, we show how Hermite schemes can be used to constructhigh-order IMEX schemes in the manner of [1]. Concerning notation, we rep-resent the m -th time derivative of g ( y ( t n ) , z ( t n )) as g ( m ) n . As an example, let5s first consider the 6th-order Hermite scheme as a base for solving Equation(1), ∆ y = ∆ t z + z ) − ∆ t
10 ( g − g ) + ∆ t
120 ( g (1)1 + g (1)0 )∆ z = ∆ t (cid:15) ( g + g ) − ∆ t (cid:15) ( g (1)1 − g (1)0 ) + ∆ t (cid:15) ( g (2)1 + g (2)0 ) . (14)6he corresponding Hermite IMEX method is as follows. Algorithm 1:
For the solution of Equation (1), we propose the follow-ing 6th-order Hermite IMEX method to advance the solution from t n to t n +1 Predict.
Given the solution ( y n , z n ), we compute a 3rd-order IMEXTaylor approximation y [0] = y n + ∆ tz n + ∆ t (cid:15) g n + ∆ t (cid:15) g (1) n (15) z [0] = z n + ∆ t(cid:15) g [0] − ∆ t (cid:15) g (1)[0] + ∆ t (cid:15) g (2)[0] (16)for the unknowns y [0] and z [0] that will be initial guesses for ourapproximation to y n +1 and z n +1 . For the 6th-order Hermite IMEXscheme, we use a 3rd-order forward Taylor series in y and a 3rd-orderbackwards Taylor series in z . For a (2 n )th-order Hermite IMEXscheme we would use an n th-order Taylor series in this step.2. Correct.
Based on the initial step, for 0 ≤ k ≤ k max − y [ k +1] = y n + ∆ t z [ k ] + z n ) − ∆ t (cid:15) ( g [ k ] − g n ) + ∆ t (cid:15) ( g (1)[ k ] + g (1) n ) (17) z [ k +1] = z n + ∆ t (cid:15) ( g [ k ] + g n ) − ∆ t (cid:15) ( g (1)[ k ] − g (1) n ) + ∆ t (cid:15) ( g (2)[ k ] + g (2) n )+ ∆ t(cid:15) ( g [ k +1] − g [ k ] ) − ∆ t (cid:15) ( g (1)[ k +1] − g (1)[ k ] ) + ∆ t (cid:15) ( g (2)[ k +1] − g (2)[ k ] ) (18)for y [ k max ] and z [ k max ] .3. Update.
The update for the solution is then carried out as y n +1 = y [ k max ] , z n +1 = z [ k max ] ..Note that the corrector for y is simply the Hermite corrector, and thatas y [ k ] and z [ k ] converge, g ( m )[ k +1] = g ( m )[ k ] and the corrector for z converges tothe Hermite scheme. Thus, it is simple to implement any of the methods inAppendix A or the natural extensions thereof in IMEX form. Subsequently,we will extend this algorithm to arbitrary splittings, following [1]. We shall7o so for the 8th-order Hermite IMEX scheme in order to demonstrate howeasily one can translate a given Hermite scheme to a Hermite IMEX scheme.One may note that higher-order derivatives of g are used to calculatethe update to z than for y , and that one has enough information to use ahigher-order quadrature when updating y , although this would not changethe overall order of the scheme. However, when modelling physical systems,there are often physical quantities that ought to be conserved, but are notconserved by a given algorithm. One such system is a Keplerian orbit, wherethe Laplace-Runge-Lenz vector should be conserved, but is not when thesystem is treated by many time-stepping algorithms such as the one outlinedin Algorithm 1, even when those algorithms conserve energy and the scalareccentricity of the orbit. In such cases one may tune the truncation error of y using higher-order derivatives of g in order to improve the conservation ofsuch quantities, as outlined in [14].We now consider an 8th-order Hermite IMEX scheme, using Equation (3).8 lgorithm 2: For general problems, we propose the following 8th orderHermite IMEX method to advance the solution of Equation (3) from t n to t n +1 , emphasising that one may construct a Hermite IMEX schemeof any even order in this manner.1. Predict.
Given the solution u n , we compute a 4th-order IMEXTaylor approximation to u [0] , u [0] = u n + ∆ t (cid:0) φ E ( u n ) + φ I ( u [0] ) (cid:1) + ∆ t (cid:16) φ (1) E ( u n ) − φ (1) I ( u [0] ) (cid:17) + ∆ t (cid:16) φ (2) E ( u n ) + φ (2) I ( u [0] ) (cid:17) + ∆ t (cid:16) φ (3) E ( u n ) − φ (3) I ( u [0] ) (cid:17) , (19)performing a forward Taylor expansion for φ E and a backward Taylorexpansion for φ I .2. Correct.
Based on the initial step, for 0 ≤ k ≤ k max − u [ k +1] = u n + ∆ t (cid:0) φ I ( u [ k +1] ) − φ I ( u [ k ] ) (cid:1) − ∆ t (cid:16) φ (1) I ( u [ k +1] ) − φ (1) I ( u [ k ] ) (cid:17) + ∆ t (cid:16) φ (2) I ( u [ k +1] ) − φ (2) I ( u [ k ] ) (cid:17) − ∆ t (cid:16) φ (3) I ( u [ k +1] ) − φ (3) I ( u [ k ] ) (cid:17) + ∆ t (cid:0) φ ( u [ k ] + φ ( u n ) (cid:1) − t (cid:0) φ (1) ( u [ k ] ) − φ (1) ( u n ) (cid:1) + ∆ t (cid:0) φ (2) ( u [ k ] ) + φ (2) ( u n ) (cid:1) − ∆ t (cid:0) φ (3) ( u [ k ] − φ (3) ( u n ) (cid:1) (20)for u [ k max ] .3. Update.
The update for the solution is then defined as u n +1 = u [ k max ] .Note that in order to construct this scheme, one simply includes onemore term in the Taylor series during prediction and correction, and uses ahigher-order Hermite scheme during the correction step. We emphasise that,similarly to the 4th-order method presented in [1], intermediate steps do notneed to be stored, and that the algorithm only requires storage of values at t n and t n +1 . This is advantageous over multistep methods as well as mostRunge-Kutta methods. 9 .2. Accuracy each iteration In this section, we review how the accuracy of the approximate solution z [ k max ] improves each iteration. It is clear that Algorithms 1 and 2 have thesame order of accuracy as their predictor when k max = 0, and approach theorder of accuracy of the underlying Hermite quadrature rule as k max → ∞ .Here, we show that with each iteration k > z as in Algorithm 1, although the same result holds for y and forproblems cast in the form of Algorithm 2. First, we assume that g ( y, z ) andits derivatives are Lipschitz continuous, such that (cid:107) g ( y , z ) − g ( y , z ) (cid:107) ≤ L g (cid:107) y − y , z − z (cid:107) , (cid:107) g (1) ( y , z ) − g (1) ( y , z ) (cid:107) ≤ L g (1) (cid:107) y − y , z − z (cid:107) , (cid:107) g (2) ( y , z ) − g (2) ( y , z ) (cid:107) ≤ L g (2) (cid:107) y − y , z − z (cid:107) , (21)where L g , L g (1) , and L g (2) are constants that subsume the respective Lipschitzconstants and factors of (cid:15) , and Euclidean norms are denoted by (cid:107) x (cid:107) .Then, consider a 6th-order Hermite quadrature rule for a generic function f , such that I [ f n , f [ k ] ] = ∆ t (cid:0) f [ k ] + f n (cid:1) − ∆ t (cid:16) f (1)[ k ] − f (1) n (cid:17) + ∆ t (cid:16) f (2)[ k ] + f (2) n (cid:17) (22)and such that I [ g n , g n +1 ] = (cid:90) t n +1 t n g ( y, z ) dt + O (∆ t ) . (23)In this notation, one updates z according to z [ k +1] = z n + ∆ t(cid:15) (cid:0) g [ k +1] − g [ k ] (cid:1) − ∆ t (cid:15) (cid:16) g (1)[ k +1] − g (1)[ k ] (cid:17) + ∆ t (cid:15) (cid:16) g (2)[ k +1] − g (2)[ k ] (cid:17) + 1 (cid:15) I [ g n , g [ k ] ] . (24)Then, assuming that z n +1 is the true solution at time t n +1 , the error of theapproximate solution z [ k ] is δ z : k = z [ k ] − z n +1 and the error of y [ k ] is δ y : k =10 [ k ] − y n +1 . Then, the error of the approximate solution is δ k = (cid:107) ( δ y : k , δ z : k ) (cid:107) .Thus, the error of the quadrature rule given by Equation (22) at iteration k is given by (cid:12)(cid:12) I [ z n , z n +1 ] − I [ z n , z [ k ] ] (cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) ∆ t (cid:0) z n +1 − z [ k ] (cid:1) − ∆ t (cid:16) z (1) n +1 − z (1)[ k ] (cid:17) + ∆ t (cid:16) z (2) n +1 − z (2)[ k ] (cid:17)(cid:12)(cid:12)(cid:12)(cid:12) ≤ ∆ t (cid:12)(cid:12) z n +1 − z [ k ] (cid:12)(cid:12) + ∆ t (cid:12)(cid:12)(cid:12) z (1) n +1 − z (1)[ k ] (cid:12)(cid:12)(cid:12) + ∆ t (cid:12)(cid:12)(cid:12) z (2) n +1 − z (2)[ k ] (cid:12)(cid:12)(cid:12) ≤ ∆ t δ k + ∆ t L g δ k + ∆ t L g (1) δ k (25)and (cid:12)(cid:12) I [ g n , g n +1 ] − I [ g n , g [ k ] ] (cid:12)(cid:12) ≤ ∆ t L g δ k + ∆ t L g [1] δ k + ∆ t L g (2) δ k (26)for integrating z and g respectively. Then, since the true solution satisfies z n +1 = z n + 1 (cid:15) (cid:90) t n +1 t n g ( y, z ) dt, (27)the error of the approximate solution is given by | δ z : k +1 | = (cid:12)(cid:12)(cid:12)(cid:12) ∆ t(cid:15) (cid:0) g [ k +1] − g [ k ] (cid:1) + ∆ t (cid:15) (cid:16) g (1)[ k +1] − g (1)[ k ] (cid:17) + ∆ t (cid:15) (cid:16) g (2)[ k +1] − g (2)[ k ] (cid:17) + 1 (cid:15) I [ g n , g n +1 ] − (cid:15) (cid:90) t n +1 t n g ( y, z ) dt (cid:12)(cid:12)(cid:12)(cid:12) ≤ ∆ t(cid:15) (cid:12)(cid:12) g [ k +1] − g [ k ] (cid:12)(cid:12)(cid:124) (cid:123)(cid:122) (cid:125) I + ∆ t (cid:15) (cid:12)(cid:12)(cid:12) g (1)[ k +1] − g (1)[ k ] (cid:12)(cid:12)(cid:12)(cid:124) (cid:123)(cid:122) (cid:125) II + ∆ t (cid:15) (cid:12)(cid:12)(cid:12) g (2)[ k +1] − g (2)[ k ] (cid:12)(cid:12)(cid:12)(cid:124) (cid:123)(cid:122) (cid:125) III + 1 (cid:15) (cid:12)(cid:12)(cid:12)(cid:12) I [ g n , g [ k ] ] − (cid:90) t n +1 t n g dt (cid:12)(cid:12)(cid:12)(cid:12)(cid:124) (cid:123)(cid:122) (cid:125) IV . (28)We treat each part separately,I ≤ | g n +1 − g [ k +1] | + | g n +1 − g [ k ] | ≤ L g δ k + L g δ k +1 , (29)II ≤ | g (1) n +1 − g (1)[ k +1] | + | g (1) n +1 − g (1)[ k ] | ≤ L g (1) δ k + L g (1) δ k +1 , (30)11II ≤ | g (2) n +1 − g (2)[ k +1] | + | g (2) n +1 − g (2)[ k ] | ≤ L g (2) δ k + L g (2) δ k +1 , (31)IV ≤ (cid:12)(cid:12) I [ g n , g [ k ] ] − I [ g n , g n +1 ] (cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) I [ g n , g n +1 ] − (cid:90) t n +1 t n g dt (cid:12)(cid:12)(cid:12)(cid:12) ≤ ∆ t L g δ k + ∆ t L g (1) δ k + ∆ t L g (2) δ k + O (∆ t ) . (32)Thus, | δ z : k +1 | ≤ ∆ t(cid:15) L g ( δ k + δ k +1 ) + ∆ t (cid:15) L g (1) ( δ k + δ k +1 )+ ∆ t (cid:15) L g (2) ( δ k + δ k +1 ) + ∆ t (cid:15) L g ( δ k ) + ∆ t (cid:15) L g (1) ( δ k )+ ∆ t (cid:15) L g (2) ( δ k ) + O (∆ t ) . (33)Note that regardless of the order of the Hermite quadrature, the limitingterm determining the order of δ z : k each iteration only changes by one factorof ∆ t until it reaches the underlying order of the quadrature rule. That is tosay, no matter what order the overall scheme, the order of the algorithm atiteration k + 1 is given by | δ k +1 | = O (∆ tδ k +1 ) + O (∆ tδ k ) + O (∆ t ) . (34)Thus, for a (2 n ) th -order scheme with an ( n ) th -order accurate predictor, asdescribed above, k max = n is sufficient to achieve the maximal order of thealgorithm. However, for some problems where these algorithms experienceorder reduction, or for problems where the underlying system has symmetryin time that is preferable to maintain numerically, more iterations may beuseful [1, 14].
3. Numerical Experiments
In this section, we apply Hermite IMEX schemes of order 6–12 to a selec-tion of test problems. We begin in Section 3.1 by examining the van der Poloscillator and making comparisons to the 4th-order algorithm with respect toefficiency. We find that higher-order methods are more susceptible to order12eduction, particularly for very stiff problems. For very stiff problems, the4th-order method may be preferable to higher-order algorithms because itappears less susceptible to order reduction. In Section 3.2, we apply HermiteIMEX methods to a stiff linear ODE and find less-significant order reduction.When applying the 6th-order scheme to Burgers’ equation in Section 3.3, weobserver order reduction in cases when the time scale for the diffusive termis more than ∼
30 times smaller than that for the convective term.
To directly compare our methods to that presented in [1], we consider thevan der Pol oscillator, which has the form of Equation (1) with g ( y, z ) = (1 − y ) z − y. (35)We use the same initial conditions as [1] when testing the 6th-order method, y = 2 , z = −
23 + 1081 (cid:15) − (cid:15) , (36)which are a common choice in the literature [e.g. 16, 18] and facilitate directcomparison with the results of [1]. When testing the 8th-order method, weuse the initial conditions y = 2 , z = −
23 + 1081 (cid:15) − (cid:15) + 1526659049 (cid:15) , (37)which ensures that g (3) ( y , z ) → (cid:15) →
0, such that the system recoversthe correct behavior in the asymptotic limit. We integrate these equationsuntil t end = 0 .
5. Errors are calculated by comparing y ( t end ) and z ( t end ) toreference values y ref and z ref , computed using an additional integration usinga time step half as small as that for the shortest-timestep result presented.We compare the values at t end to the reference values using the Euclideannorm δ ≡ (cid:107) y ( t end ) − y ref , z ( t end ) − z ref (cid:107) .We plot results for this test in Figures 1 and 2. The schemes, whenapplied using the minimum k max to achieve their nominal order, suffer fromsignificant order reduction for some values of (cid:15) . However, by increasing k max to 20 or 100, order reduction can be fully or partially ameliorated. Whencomparing the efficiency of the the 6th- and 8th-order schemes to the 4th-order version, we primarily consider the cost of calculating higher-order timederivatives. Thus, we approximate that the 6th-order scheme requires about13wice the total floating point operations of the 4th-order scheme, and the8th-order scheme requires about ∼ . k max required to reach the nominal order of thescheme differ between schemes, we compare the other choices of k max . Con-cerning the predictor alone, k max = 0, the 4th-order algorithm presented in[1] reaches an error of ∼ − using about 17000 time steps. The 6th-orderalgorithm reaches similar error using about 500 time steps, while the 8th-order algorithm reaches similar errors using about 150 time steps. Thus,in this scenario, the 6th-order scheme is roughly ∼
17 times more efficientthan the 4th-order scheme, and the 8th-order scheme is roughly ∼
33 timesmore efficient than the 4th-order scheme. Naturally, the higher-order schemeswould be even more beneficial if small enough time steps were used to achieverounding-dominated errors.Concerning k max = 20, the 4th-order scheme experiences little order re-duction and reachs rounding-limited errors using around 1000 steps. The6th-order schemes suffer order reduction for some values of (cid:15) , but still achieveorders better than the predictor alone, reaching rounding-limited errors us-ing ∼ − −
500 steps depending on the value of (cid:15) . Thus, 6th-order schemeranges from roughly ∼ (cid:15) . For somevalues of (cid:15) , the 8th-order scheme reaches errors around ∼ − using just32 time steps, and is more efficient than the 6th- or 4th-order schemes, whilefor other values of (cid:15) the 8th-order scheme is less efficient. Using k max = 100,the 6th-order scheme suffers virtually no order reduction and is more efficientthan the 4th-order scheme. The 8th-order scheme still suffers order reductionfor many values of (cid:15) even when using k max = 100, and is generally about asefficient as the 6th-order scheme in this case, although this depends on (cid:15) . Although we observed order reduction down to the order of the predictorwhile applying the 8th-order Hermite IMEX method for both k max = 20 and k max = 100 for various values of (cid:15) in the van der Pol oscillator problem, weshow that this does not preclude the application of higher-order methods toall problems. We demonstrate this by examining a stiff linear system, findingthat the 10th- and 12th- order schemes performed at or around their nominalorder accuracy for k max = 20. 14 steps (N)10 e rr o r () k max = 0 = 10 = 10 = 10 = 10 = 10 steps (N)10 e rr o r () k max = 3 = 10 = 10 = 10 = 10 = 10 steps (N)10 e rr o r () k max = 20 = 10 = 10 = 10 = 10 = 10 steps (N)10 e rr o r () k max = 100 = 10 = 10 = 10 = 10 = 10 Figure 1: The convergence of 6th-order Hermite IMEX schemes applied to the van der Polequation. Dashed black lines indicate a 3rd-order scaling, and dotted black lines indicate a6th-order scaling. Results for (cid:15) = 10 − are plotted using blue circles, results for (cid:15) = 10 − are plotted using vermilion squares, results for (cid:15) = 10 − are plotted using turquoise left-pointing triangles, results for (cid:15) = 10 − are plotted using orange stars, and results for (cid:15) = 10 − are plotted using purple right-facing triangles. steps (N)10 e rr o r () k max = 0 = 10 = 10 = 10 = 10 = 10 steps (N)10 e rr o r () k max = 4 = 10 = 10 = 10 = 10 = 10 steps (N)10 e rr o r () k max = 20 = 10 = 10 = 10 = 10 = 10 steps (N)10 e rr o r () k max = 100 = 10 = 10 = 10 = 10 = 10 Figure 2: The convergence of 8th-order Hermite IMEX schemes applied to the van derPol equation. Dashed black lines indicate a 4th-order scaling, and dotted black linesindicate an 8th-order scaling. Results for (cid:15) = 10 − are plotted using blue circles, resultsfor (cid:15) = 10 − are plotted using vermilion squares, results for (cid:15) = 10 − are plotted usingturquoise left-pointing triangles, results for (cid:15) = 10 − are plotted using orange stars, andresults for (cid:15) = 10 − are plotted using purple right-facing triangles. dydt = − Ky, (38)with initial value y (0) = 1, which is trivial to solve analytically but can bedifficult to solve numerically for large values of K . We solve for the value of y at t = 0 .
5, and compare to the true value of exp( − K/
2) by calculating thefractional error. In Figure 3 we present the results for the 10th- and 12th-order algorithms using k max = 20 in all cases. We follow the methodologypresented in in Algorithm 2, neglecting the explicit components. We find thatalthough there is some scatter, the order of convergence tracks the nominalvalue until rounding errors begin to dominate. Notably, for this problem thecost of evaluating higher-order derivatives is trivial. For both values of K tested, the 12th-order algorithm is preferable for reaching small errors, butfor lenient tolerances the 10th-order algorithm may be preferable. N10 e rr o r K=10 N10 e rr o r K=50
Figure 3: Convergence of the 10th- and 12th-order versions of Algorithm 2 applied toEquation (38), using k max = 20. Results for the 12th-order algorithm are plotted usingorange circles, and results for the 10th-order algorithm are plotted using blue squares. .3. Burgers’ equation In this section, we test the efficacy of the 6th-order Hermite IMEX schemefor solving nonlinear PDEs. We solve Burgers’ equation, ∂u∂t = − u ∂u∂x + ν ∂ u∂x . (39)This equation is in the form of Algorithm 2, where the second term on theright-hand side of Equation (39) can lead to stringent time-step constraints.Specifically, for an explicit scheme, the convective term typically requiresa time step ∆ t (cid:46) ∆ x/u , while the parabolic term typically requires ∆ t (cid:46) ∆ x /ν .We solve Equation (39) according to Algorithm 2 subject to periodicboundary conditions, taking derivatives using a Fourier pseudospectral method.We use an exponential filter for de-aliasing [19], and discretize the domaininto N x = 64 and N x = 256 nodes. It is useful to investigate different spa-tial resolutions because evolving the diffusive term becomes more challeng-ing compared to the convective term as resolution increases. We calculatehigher time derivatives of Equation (39) using the Cauchy-Kowalevski, orLax-Wendroff, procedure, converting time derivatives to spatial derivativesbut limiting ourselves to problems without discontinuities in u or its deriva-tives. To complete the picture, we specify initial conditions as u ( x ) = 2 . πx/L ) (40)where L is the length of the domain, which we set to 1. We evolve the systemuntil t = 0 .
15 using various values of ν . We quantify the time step in termsof α ≡ ∆ t max( u ) / ∆ x , where ∆ x ≡ L/N x and α determines the step size.We calculate errors by running an additional simulation as a reference,using a time step half that of the smallest used in the presented results. Wethen calculate the L norm of the difference between each solution and thereference solution, plotting the results in Figure 4. Concerning the N x =64 tests, ∆ x/u ∼ / x /ν was 1 / k max = 3 in allcases. Concerning the N x = 256 tests, ∆ x/u ∼ / x /ν was 1 / ∼
128 times smaller. For more extreme timescale differences, the Hermite IMEX method converges at third order using k max = 3, although for smaller values of ν the method still converges at6th order. Even when suffering order reduction, applying a few correctoriterations still reduces the error by more than an order of magnitude.18 L e rr o r N x = 64 = 1.0= 0.5= 0.1= 0.01 10 L e rr o r N x = 256 = 1.0, k max = 0= 1.0= 0.5= 0.1= 0.01 Figure 4: L errors as a function of time step when solving Burgers’ equation using a6th-order Hermite IMEX scheme. The left panel plots results computed using a resolutionof N x = 64, while the right panel plots results using a spatial resolution of N x = 256. Inboth panels, blue squares plot results for ν = 1, orange circles plot results for ν = 0 . ν = 0 .
1, and purple left-facing triangles plotresults for ν = 0 .
01. Solid markers indicate results that were computed using k max = 3,while open symbols indicate results computed using k max = 0. Dotted black lines indicatea 3rd-order scaling, while dashed black lines indicate a 6th-order scaling.
4. Conclusions
In this work, we have presented a class of multiderivative IMEX schemesof arbitrary order, focusing on orders 6–12 following recent work on the 4th-order method of this family [1]. These schemes have a number of advanta-geous qualities. Although calculating higher-order derivatives can be expen-sive, such calculations are highly parallelizable, and may be advantageouson current and future parallel architectures. Additionally, many high-ordermultistep and Runge-Kutta methods require storing a number of interme-diate stages of the solution, while the methods presented here only requirestoring values at the beginning and end of a time step, regardless of the or-der of the algorithm. However, because higher-order schemes appear moresusceptible to order reduction for more stiff problems, lower-order schemes19an be more efficient depending on the problem at hand. Similarly, tuning k max adaptively may be fruitful. We note that the 4th-order version of theseschemes has been proven to be asymptotic preserving [1], although we havenot extended this proof to arbitrary orders in this work.For nonlinear problems, especially systems of nonlinear PDEs, computinghigher-order time derivatives can be very expensive. This is one reason thatthe original manner of implementing ADER schemes, based on the Cauchy-Kowalevski procedure and forward Taylor approximations [20], has fallen outof favour compared to implementations that implicitly solve a weak formu-lation of the PDE in time [21]. However, the methods presented here canhandle stiff terms implicitly, and also require significantly fewer time deriva-tives of the governing PDE to be calculated to achieve an algorithm of agiven order. Thus, the algorithms presented here have a number of benefits,and do not suffer the same drawbacks as previous very-high-order schemes. Acknowledgements
The author is thankful to Cole Miller and Derek Richardson for theirfeedback during the preparation of this manuscript, and to Bill Dorlandand Keigo Nitadori for their comments, all of which improved the qualityof this work. We gratefully acknowledge support from NASA under grantNNX17AF29G. 20 ppendix A. Hermite schemes up to 12th order u = ∆ t f + f ) − ∆ t
10 ( f (1)1 − f (1)0 ) + ∆ t
120 ( f (2)1 + f (2)0 ) (A.1)8th:∆ u = ∆ t f + f ) − t
28 ( f (1)1 − f (1)0 ) + ∆ t
84 ( f (2)1 + f (2)0 ) − ∆ t f (3)1 − f (3)0 ) (A.2)10th:∆ u = ∆ t f + f ) − ∆ t f (1)1 − f (1)0 ) + ∆ t
72 ( f (2)1 + f (2)0 ) − ∆ t f (3)1 − f (3)0 ) + ∆ t f (4)1 + f (4)0 ) (A.3)12th: ∆ u = ∆ t f + f ) − t
44 ( f (1)1 − f (1)0 ) + ∆ t
66 ( f (2)1 + f (2)0 ) − ∆ t
792 ( f (3)1 − f (3)0 ) + ∆ t f (4)1 + f (4)0 ) − ∆ t f (5)1 − f (5)0 ) (A.4) References [1] J. Sch¨utz, D. C. Seal, An asymptotic preserving semi-implicit multi-derivative solver, arXiv e-prints (2020) arXiv:2001.08268 arXiv:2001.08268 .[2] R. J. Leveque, H. C. Yee, A Study of Numerical Methods for HyperbolicConservation Laws with Stiff Source Terms, Journal of ComputationalPhysics 86 (1) (1990) 187–210. doi:10.1016/0021-9991(90)90097-K .[3] H. P. Langtangen, K.-A. Mardal, R. Winther, Numerical methodsfor incompressible viscous flow, Advances in Water Resources 25 (8)(2002) 1125 – 1146. doi:https://doi.org/10.1016/S0309-1708(02)00052-0 . 214] O. Tes , ileanu, A. Mignone, S. Massaglia, Simulating radiative astro-physical flows with the PLUTO code: a non-equilibrium, multi-speciescooling function, Astronomy and Astrophysics 488 (1) (2008) 429–440. doi:10.1051/0004-6361:200809461 .[5] R. E. O. Jr., Singular Perturbation Methods for Ordinary DifferentialEquations, 1st Edition, Applied Mathematical Sciences 89, Springer-Verlag New York, 1991.[6] B. van der Pol, J. van der Mark, Frequency Demultiplication, Nature120 (3019) (1927) 363–364. doi:10.1038/120363a0 .[7] J. Makino, Optimal Order and Time-Step Criterion for Aarseth-TypeN-Body Integrators, The Astrophysical Journal 369 (1991) 200. doi:10.1086/169751 .[8] M. Lecar, R. Loeser, J. R. Cherniack, Numerical integration of gravi-tational N-body systems with the use of explicit Taylor series., in: Nu-merical Solution of Ordinary Differential Equations, 1974, pp. 451–470.[9] S. Gottlieb, C.-W. Shu, E. Tadmor, Strong Stability-Preserving High-Order Time Discretization Methods, SIAM Review 43 (1) (2001) 89–112. doi:10.1137/S003614450036757X .[10] D. I. Ketcheson, Highly efficient strong stability-preserving rungekuttamethods with low-storage implementations, SIAM Journal on ScientificComputing 30 (4) (2008) 2113–2136. doi:10.1137/07070485X .[11] A. Dutt, L. Greengard, V. Rokhlin, Spectral deferred correction methodsfor ordinary differential equations, BIT Numerical Mathematics 40 (2)(2000) 241–266. doi:10.1023/A:1022338906936 .[12] A. Christlieb, B. Ong, J.-M. Qiu, Comments on high-order integratorsembedded within integral deferred correction methods, Commun. Appl.Math. Comput. Sci. 4 (1) (2009) 27–56. doi:10.2140/camcos.2009.4.27 .URL https://doi.org/10.2140/camcos.2009.4.27 [13] K. Nitadori, Hermite integrators and Riordan arrays: When one feelstired in round-off error (Jan 2016).22RL [14] A. J. Dittmann, Modified Hermite integrators of arbitrary order,Monthly Notices of the Royal Astronomical Society 496 (2) (2020) 1217–1223. doi:10.1093/mnras/staa1631 .[15] K. Nitadori, J. Makino, Sixth- and eighth-order Hermite integrator forN-body simulations, New Astronomy 13 (7) (2008) 498–507. arXiv:0708.0738 , doi:10.1016/j.newast.2008.01.010 .[16] G. W. E. Hairer, Solving Ordinary Differential Equations II: Stiffand Differential-Algebraic Problems, Springer Series in ComputationalMathematic, Springer-Verlag, 1991.[17] J. D. C. a. J. Kevorkian, Perturbation methods in applied mathematics,1st Edition, Applied Mathematical Sciences 34, Springer-Verlag NewYork, 1981.[18] S. Boscarino, On an accurate third order implicit-explicit rungekuttamethod for stiff problems, Applied Numerical Mathematics 59 (7) (2009)1515 – 1528. doi:https://doi.org/10.1016/j.apnum.2008.10.003 .[19] A. Majda, J. McDonough, S. Osher, The fourier method for nonsmoothinitial data, Mathematics of Computation 32 (144) (1978) 1041–1081. doi:10.1090/S0025-5718-1978-0501995-4 .[20] V. Titarev, E. Toro, Ader: Arbitrary high order godunov approach, J.Sci. Comput. 17 (2002) 609–618. doi:10.1023/A:1015126814947 .[21] M. Dumbser, C. Enaux, E. F. Toro, Finite volume schemes of very highorder of accuracy for stiff hyperbolic balance laws, Journal of Compu-tational Physics 227 (8) (2008) 3971–4001. doi:10.1016/j.jcp.2007.12.005doi:10.1016/j.jcp.2007.12.005