aa r X i v : . [ m a t h . NA ] M a y Volume-preserving exponential integrators
Bin Wang ∗ Xinyuan Wu † May 31, 2018
Abstract
As is known that various dynamical systems including all Hamiltonian systems preservevolume in phase space. This qualitative geometrical property of the analytical solution shouldbe respected in the sense of Geometric Integration. This paper analyses the volume-preservingproperty of exponential integrators in different vector fields. We derive a necessary and sufficientcondition of volume preservation for exponential integrators, and with this condition, volume-preserving exponential integrators are analysed in detail for four kinds of vector fields. It turnsout that symplectic exponential integrators can be volume preserving for a much larger class ofvector fields than Hamiltonian systems. On the basis of the analysis, novel volume-preservingexponential integrators are derived for solving highly oscillatory second-order systems andextended Runge–Kutta–Nystr¨om (ERKN) integrators of volume preservation are presented forseparable partitioned systems. Moreover, the volume preservation of Runge–Kutta–Nystr¨om(RKN) methods is also discussed. Four illustrative numerical experiments are carried out todemonstrate the notable superiority of volume-preserving exponential integrators in comparisonwith volume-preserving Runge-Kutta methods.Keywords: exponential integrators, volume preservation, geometric integrators, extended RKNintegrators, highly oscillatory systemsMSC (2000): 65L05, 65L06, 65L99, 34C60
Geometric integrators (also called as structure-preserving algorithms) have been an area of greatinterest and active research in the recent decades. The main advantage of such methods for solvingordinary differential equations (ODEs) is that they can exactly preserve some qualitative geometricalproperty of the analytical solution, such as the symplecticity, energy preservation, and symmetryfor long-term integration. Various geometric integrators have been designed and analysed recentlyand the reader is referred to [2, 3, 5, 9, 10, 12, 19, 21, 25, 27] for example on this topic. For agood theoretical foundation of geometric numerical integration for ODEs, we refer the reader to[7, 13]. Surveys of structure-preserving algorithms for oscillatory differential equations are referredto [31, 33]. ∗ School of Mathematical Sciences, Qufu Normal University, Qufu 273165, P.R.China; Mathematisches Institut,University of T¨ubingen, Auf der Morgenstelle 10, 72076 T¨ubingen, Germany. E-mail: [email protected] † School of Mathematical Sciences, Qufu Normal University, Qufu 273165, P.R.China; Department of Mathematics,Nanjing University, Nanjing 210093, P.R. China. E-mail: [email protected]
1t is well known that volume preservation is an important property shared by several dynamicalsystems. By the classical theorem due to Liouville, it is clear that all Hamiltonian systems arealso of volume preservation. Preservation of volume by a numerical method has become a desirableproperty and many methods have been proposed and shown to be (or not to be) of volume preser-vation. We refer the reader to [4, 8, 14, 18, 20, 24, 35] and references therein. From the researchesof this topic, it follows that all symplectic methods are of volume preservation for Hamiltoniansystems. However, this result does not hold for the system beyond Hamiltonian systems. Feng andShang have proved in [8] that no Runge-Kutta (RK) method can be of volume preservation evenfor linear divergence free vector fields. The authors in [4, 18] showed that no B-series method canbe of volume preservation for all possible divergence-free vector fields. Thus, the design of efficientvolume-preserving (VP) methods is still an open problem in the geometric numerical integration(see [22]). Recently, various VP methods have been constructed and analysed for different vec-tor fields, such as splitting methods (see [20, 35]), RK methods (see [1]) and the methods usinggenerating functions (see [24, 36]).On the other hand, exponential integrators have been proposed and researched as an efficientapproach to integrating ODEs/PDEs. The reader is referred to [15, 16, 17, 29] for example. How-ever, it seems that exponential integrators with volume-preserving property for different vectorfields have not been researched yet in the literature, which motives this paper.In this paper, we study the volume preservation of exponential integrators when solving thefollowing first-order ODEs y ′ ( t ) = Ky ( t ) + g ( y ( t )) := f ( y ( t )) , y (0) = y ∈ R n , (1)where K is an n × n matrix which is assumed that (cid:12)(cid:12) e hK (cid:12)(cid:12) = − < h <
1, and g : R n → R n is a differentiable nonlinear function. In this paper, |·| denotes the determinant. The function f is assumed to be divergence free such that this system is of volume preservation. It is well knownthat the exact solution of (1) can be presented by the variation-of-constants formula y ( t ) = e tK y + t Z e (1 − τ ) tK g ( y ( τ t )) dτ. (2)The main contributions of this paper are to derive the volume-preserving condition for exponen-tial integrators and analyse their volume preservations for different classes of vector fields which arelarger than Hamiltonian systems. Furthermore, based on the analysis, volume-preserving adaptedexponential integrators are formulated for highly oscillatory second-order systems and extendedRunge–Kutta–Nystr¨om (ERKN) integrators of volume preservation are derived for separable parti-tioned systems. We also discuss the volume preservation of Runge–Kutta–Nystr¨om (RKN) methodsby considering them as a special class of ERKN integrators. To our knowledge, the results pre-sented in this paper are novel that rigorously study the volume-preserving properties of exponentialintegrators and ERKN/RKN integrators.We organise the remainder of this paper as follows. In Section 2, the scheme of exponentialintegrators is presented and some useful results of these integrators are summarised. Then a neces-sary and sufficient condition for exponential integrators to be of volume preservation is derived inSection 3. On the basis of this condition, we study the volume-preserving properties of exponentialintegrators for four kinds of vector fields in Section 4. Section 5 discusses the application to variousproblems including highly oscillatory second-order systems and separable partitioned systems, andshows the volume preservation of adapted exponential integrators and ERKN/RKN integrators for2hese different problems. Four illustrative numerical experiments are implemented in Section 6 andthe concluding remarks are included in Section 7. In order to solve (1) effectively, one needs to approximate the integral appearing in (2) by a quadra-ture formula with suitable nodes c , c , . . . , c s . This leads to the following exponential integratorsproposed in [16], which have been successfully used for solving different kinds of ODEs/PDEs. Definition 2.1 (See [16]) An s -stage exponential integrator for numerical integration of (1) isdefined by k i = e c i hK y n + h s P j =1 ¯ a ij ( hK ) g ( k j ) , i = 1 , , . . . , s,y n +1 = e hK y n + h s P i =1 ¯ b i ( hK ) g ( k i ) , (3) where h is a stepsize, c i are real constants for i = 1 , · · · , s , ¯ b i ( hK ) and ¯ a ij ( hK ) are matrix-valuedfunctions of hK for i, j = 1 , . . . , s . The coefficients of the integrator can be displayed in a Butcher tableau (omit ( hK ) for brevity): c ¯ A ¯ b ⊺ = c ¯ a . . . ¯ a s ... ... . . . ... c s ¯ a s · · · ¯ a ss ¯ b · · · ¯ b s It is noted that when K = , this integrator reduces to a classical s -stage RK method representedby the Butcher tableau c Ab ⊺ = c a . . . a s ... ... . . . ... c s a s · · · a ss b · · · b s In this paper, a kind of special and important exponential integrators will be considered andanalysed, which was proposed in [23].
Definition 2.2 (See [23]) Define a kind of s -stage exponential integrators by ¯ a ij ( hK ) = a ij e ( c i − c j ) hK , ¯ b i ( hK ) = b i e (1 − c i ) hK , i, j = 1 , . . . , s, (4) where c = ( c , . . . , c s ) ⊺ , b = ( b , . . . , b s ) ⊺ , A = ( a ij ) s × s (5) are the coefficients of an s -stage RK method. With regard to this kind of exponential integrators, two useful properties are shown in [23] andwe summarise them as follows. 3 heorem 2.3 (See [23]) If an RK method with the coefficients (5) is of order p , then the expo-nential integrator given by (4) is also of order p . Theorem 2.4 (See [23]) The exponential integrator defined by (4) is symplectic if the correspondingRK method (5) is symplectic.
In this paper, we supplement an additional requirement to b of (5) and define the followingspecial symplectic exponential integrators (SSEI). Definition 2.5 An s -stage exponential integrator (4) is called as special symplectic exponentialintegrator (SSEI) if the RK method (5) is symplectic and b j = 0 for all j = 1 , . . . , s. Remark 2.6
We note that a kind of special symplectic RK (SSRK) methods was considered in [1]and our SSEI integrators reduce to the SSRK methods when K = . For each stepsize h , denote the s -stage exponential integrator (3) by a map ϕ h : R n → R n , which is ϕ h ( y ) = e hK y + h s P i =1 ¯ b i ( hK ) g ( k i ( y )) ,k i ( y ) = e c i hK y + h s P j =1 ¯ a ij ( hK ) g ( k j ( y )) , i = 1 , , . . . , s. (6)It is noted that k i ( y ) for i = 1 , , . . . are identical to ϕ c i h ( y ) for i = 1 , , . . . in (6).We firstly derive the Jacobian matrix and its determinant for the general exponential integrator(6). Lemma 3.1
The Jacobian matrix of the exponential integrator (6) can be expressed as ϕ ′ h ( y ) = e hK + h ¯ b ⊺ F ( I s ⊗ I − h ¯ AF ) − e chK , where F = diag ( g ′ ( k ) , . . . , g ′ ( k s )) , I s and I are the s × s and n × n identity matrices, respectively,and e chK = ( e c hK , . . . , e c s hK ) ⊺ . Its determinant reads | ϕ ′ h ( y ) | = (cid:12)(cid:12) e hK (cid:12)(cid:12) (cid:12)(cid:12) I s ⊗ I − h ( ¯ A − e ( c − ) hK ¯ b ⊺ ) F (cid:12)(cid:12)(cid:12)(cid:12) I s ⊗ I − h ¯ AF (cid:12)(cid:12) , (7) where is an s × vector of units and e ( c − ) hK = ( e ( c − hK , . . . , e ( c s − hK ) ⊺ . Here we make useof the Kronecker product ⊗ throughout this paper. Proof According to the first formula of (6), we obtain ϕ ′ h ( y ) = e hK + h s X i =1 ¯ b i g ′ ( k i ( y )) k ′ i ( y ) = e hK + h ¯ b ⊺ F ( k ′ , . . . , k ′ s ) ⊺ . (8)4ikewise, it follows from k i ( y ) in (6) that I − h ¯ a g ′ ( k ) − h ¯ a g ′ ( k ) · · · − h ¯ a s g ′ ( k s ) − h ¯ a g ′ ( k ) I − h ¯ a g ′ ( k ) · · · − h ¯ a s g ′ ( k s )... ... ... ... − h ¯ a s g ′ ( k ) − h ¯ a s g ′ ( k ) · · · I − h ¯ a ss g ′ ( k s ) k ′ k ′ ... k ′ s = e chK , which can be rewritten as ( I s ⊗ I − h ¯ AF )( k ′ , . . . , k ′ s ) ⊺ = e chK . (9)Substituting (9) into (8) yields the first statement of this lemma.For the second statement, we will use the following block determinant identity (see [1, 13]): | U | (cid:12)(cid:12) X − W U − V (cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) U VW X (cid:12)(cid:12)(cid:12)(cid:12) = | X | (cid:12)(cid:12) U − V X − W (cid:12)(cid:12) , which is yielded from Gaussian elimination. By letting X = e hK , W = − h ¯ b ⊺ F, U = I s ⊗ I − h ¯ AF, V = e chK , it is obtained that (cid:12)(cid:12) I s ⊗ I − h ¯ AF (cid:12)(cid:12) | ϕ ′ h ( y ) | = (cid:12)(cid:12) e hK (cid:12)(cid:12) (cid:12)(cid:12) I s ⊗ I − h ¯ AF + he chK e − hK ¯ b ⊺ F (cid:12)(cid:12) = (cid:12)(cid:12) e hK (cid:12)(cid:12) (cid:12)(cid:12) I s ⊗ I − h ( ¯ A − e ( c − ) hK ¯ b ⊺ ) F (cid:12)(cid:12) , which leads to the second result (7).For the SSEI methods of Definition 2.5, a necessary and sufficient condition for these integratorsto be of volume preservation is shown by the following lemma. Lemma 3.2 An s -stage SSEI method defined in Definition 2.5 is of volume preservation if andonly if the following VP condition is satisfied | I s ⊗ I − h ( A ⊗ I. ∗ E ( hK )) F | = (cid:12)(cid:12) e hK (cid:12)(cid:12) | I s ⊗ I + h ( A ⊺ ⊗ I. ∗ E ( hK )) F | , (10) where E ( hK ) is a block matrix defined by E ( hK ) = ( E i,j ( hK )) s × s = I e ( c − c ) hK · · · e ( c − c s ) hK e ( c − c ) hK I · · · e ( c − c s ) hK ... ... . . . ... e ( c s − c ) hK e ( c s − c ) hK · · · I , (11) and . ∗ denotes the element-wise multiplication of two matrices. Proof In terms of the choice (4) of the coefficients, it is computed that¯ A − e ( c − ) hK ¯ b ⊺ = ( A ⊗ I ) . ∗ E ( hK ) − ( e ( c − hK , . . . , e ( c s − hK ) ⊺ ( b e (1 − c ) hK , . . . , b s e (1 − c s ) hK )= ( A ⊗ I ) . ∗ E ( hK ) − ( b ⊺ ⊗ I ) . ∗ E ( hK )= ( A − b ⊺ ) ⊗ I. ∗ E ( hK ) . | ϕ ′ h ( y ) | = (cid:12)(cid:12) e hK (cid:12)(cid:12) | I s ⊗ I − h ( A − b ⊺ ) ⊗ I. ∗ E ( hK ) F | (cid:12)(cid:12) I s ⊗ I − h ¯ AF (cid:12)(cid:12) . (12)Moreover, with careful calculations, it can be verified that for B = diag( b , . . . , b s ), the followingresult holds | I s ⊗ I − h ( A − b ⊺ ) ⊗ I. ∗ E ( hK ) F | = (cid:12)(cid:12) I s ⊗ I − h ( B ⊗ I )( A − b ⊺ ) ⊗ I. ∗ E ( hK ) F ( B − ⊗ I ) (cid:12)(cid:12) = (cid:12)(cid:12) I s ⊗ I − h ( B ⊗ I )( A − b ⊺ ) ⊗ I. ∗ E ( hK )( B − ⊗ I ) F (cid:12)(cid:12) = (cid:12)(cid:12) I s ⊗ I − hB ( A − b ⊺ ) B − ⊗ I. ∗ E ( hK ) F (cid:12)(cid:12) . (13)Since the RK method is symplectic, one has that BA + A ⊺ B − bb ⊺ = 0 (see [13]), which leads to B ( A − b ⊺ ) B − = − A ⊺ . Therefore, the result (13) can be simplified as | I s ⊗ I − h ( A − b ⊺ ) ⊗ I. ∗ E ( hK ) F | = | I s ⊗ I + h ( A ⊺ ⊗ I. ∗ E ( hK )) F | . The proof is complete by considering (12).
Remark 3.3
It is noted that when K = , this VP condition (10) reduces to the condition of RKmethods presented in [1]. Consequently, the condition (10) can be regarded as a generalisation ofthat of RK methods. In this section, we study the volume-preserving properties of exponential integrators for four kindsof vector fields, which are defined as follows.
Definition 4.1 (See [1]) Define the following four classes of vector fields on Euclidean space usingvector fields f ( y ) H = { f | there exists P such that for all y, P f ′ ( y ) P − = − f ′ ( y ) ⊺ } , S = { f | there exists P such that for all y, P f ′ ( y ) P − = − f ′ ( y ) } , F ( ∞ ) = { f ( y , y ) = ( u ( y ) , v ( y , y )) ⊺ where u ∈ H ∪ F ( ∞ ) | there exists P such that for all y , y , P ∂ y v ( y , y ) P − = − ∂ y v ( y , y ) ⊺ } , F (2) = { f ( y , y ) = ( u ( y ) , v ( y , y )) ⊺ where u ∈ H ∪ S ∪ F (2) | there exists P such that for all y , y , either P ∂ y v ( y , y ) P − = − ∂ y v ( y , y ) ⊺ , or P ∂ y v ( y , y ) P − = − ∂ y v ( y , y ) } . Remark 4.2
It has been proved in [1] that all these sets are equal to divergence free vector fields.The relationships of these vector fields are also given in [1] by
H ⊂ F ( ∞ ) ⊂ F (2) and S ⊂ F ( ∞ ) ⊂F (2) . According to Lemma 3.2 of [1], we know that the set H contains all Hamiltonian systems.Denote by H the set of Hamiltonian systems. See Figure 1 for the venn diagram illusting therelationships. This figure clearly shows that the sets H , F ( ∞ ) and F (2) are larger classes of vectorfields than Hamiltonian systems. H Theorem 4.3
All SSEI methods for solving (1) are of volume preservation for vector fields f and g in H with the same P . Proof Let P be such that for all y , P f ′ ( y ) P − = − f ′ ( y ) ⊺ and P g ′ ( y ) P − = − g ′ ( y ) ⊺ . Accordingto this condition and the expression f ( y ) = Ky + g ( y ), one has that P KP − = − K ⊺ . Thus it iseasily obtained that P e hK P − = e − hK ⊺ . In the light of this result, one gets (cid:12)(cid:12) P e hK P − (cid:12)(cid:12) = (cid:12)(cid:12) e hK (cid:12)(cid:12) = (cid:12)(cid:12) e − hK ⊺ (cid:12)(cid:12) = (cid:12)(cid:12) e − hK (cid:12)(cid:12) = (cid:12)(cid:12) ( e hK ) − (cid:12)(cid:12) = | e hK | , which yields (cid:12)(cid:12) e hK (cid:12)(cid:12) = 1 (it is assumed that (cid:12)(cid:12) e hK (cid:12)(cid:12) = − | I s ⊗ I − h ( A ⊗ I. ∗ E ( hK )) F | = (cid:12)(cid:12) ( I s ⊗ P )( I s ⊗ P − ) − h ( I s ⊗ P )( A ⊗ I. ∗ E ( hK ))( I s ⊗ P − )( I s ⊗ P ) F ( I s ⊗ P − ) (cid:12)(cid:12) = (cid:12)(cid:12) I s ⊗ I + h ( I s ⊗ P )( A ⊗ I. ∗ E ( hK ))( I s ⊗ P − ) F ⊺ (cid:12)(cid:12) = | I s ⊗ I + h ( A ⊗ I. ∗ E ( − hK ⊺ )) F ⊺ | = | I s ⊗ I + hF ( A ⊺ ⊗ I. ∗ E ( − hK ⊺ ) ⊺ ) | (transpose)= | I s ⊗ I + h ( A ⊺ ⊗ I. ∗ E ( − hK ⊺ ) ⊺ ) F | (Sylvester’s law) . It follows from the definition (11) that E ( − hK ⊺ ) ⊺ = ( E i,j ( − hK ⊺ )) ⊺ s × s = ( E j,i ( − hK )) s × s = ( E i,j ( hK )) s × s = E ( hK ) . (14)Therefore, one arrives at | I s ⊗ I − h ( A ⊗ I. ∗ E ( hK )) F | = | I s ⊗ I + h ( A ⊺ ⊗ I. ∗ E ( hK )) F | , which shows that all SSEI methods are of volume preservation for vector fields in H by consideringLemma 3.2. S Theorem 4.4
All one-stage SSEI methods and all two-stage SSEI methods with c = c (and anycomposition of such methods) are of volume preservation for vector fields f and g in S with thesame P . P be such that for all y , P f ′ ( y ) P − = − f ′ ( y ) and P g ′ ( y ) P − = − g ′ ( y ). Similarly tothe proof of last theorem, we obtain that P KP − = − K . Thus it is true that P e hK P − = e − hK and (cid:12)(cid:12) e hK (cid:12)(cid:12) = 1.For the one-stage SSEI, according to Lemma 3.2, it is of volume preservation if and only if | I − ha g ′ ( k ) | = | I + ha g ′ ( k ) | , which can be verified by considering | I − ha g ′ ( k ) | = (cid:12)(cid:12) P P − − ha P g ′ ( k ) P − (cid:12)(cid:12) = | I + ha g ′ ( k ) | . For a two-stage SSEI method, according to Lemma 3.2 again, this two-stage SSEI method is ofvolume preservation if and only if (cid:12)(cid:12)(cid:12)(cid:12) I − ha g ′ ( k ) − ha e ( c − c ) hK g ′ ( k ) − ha e ( c − c ) hK g ′ ( k ) I − ha g ′ ( k ) (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) I + ha g ′ ( k ) ha e ( c − c ) hK g ′ ( k ) ha e ( c − c ) hK g ′ ( k ) I + ha g ′ ( k ) (cid:12)(cid:12)(cid:12)(cid:12) , which gives the condition | I − ha g ′ ( k ) − ha g ′ ( k ) + h a a g ′ ( k ) g ′ ( k ) − h a a e ( c − c ) hK g ′ ( k ) e ( c − c ) hK g ′ ( k ) | = | I + ha g ′ ( k ) + ha g ′ ( k ) + h a a g ′ ( k ) g ′ ( k ) − h a a e ( c − c ) hK g ′ ( k ) e ( c − c ) hK g ′ ( k ) | . (15)It also can be verified thatthe left hand side of (15)= | P P − − ha P g ′ ( k ) P − − ha P g ′ ( k ) P − + h a a P g ′ ( k ) g ′ ( k ) P − − h a a P e ( c − c ) hK g ′ ( k ) e ( c − c ) hK g ′ ( k ) P − | = | I + ha g ′ ( k ) + ha g ′ ( k ) + h a a P g ′ ( k ) P − P g ′ ( k ) P − − h a a P e ( c − c ) hK P − P g ′ ( k ) P − P e ( c − c ) hK P − P g ′ ( k ) P − | = | I + ha g ′ ( k ) + ha g ′ ( k ) + h a a g ′ ( k ) g ′ ( k ) − h a a e ( c − c ) hK g ′ ( k ) e ( c − c ) hK g ′ ( k ) | . Under the assumption that c = c , this result becomes | I + ha g ′ ( k ) + ha g ′ ( k ) + h a a g ′ ( k ) g ′ ( k ) − h a a e ( c − c ) hK g ′ ( k ) e ( c − c ) hK g ′ ( k ) | . Thus (15) is obtained immediately and then all two-stage SSEI methods with c = c are of volumepreservation. Remark 4.5
It is noted that for the vector field S and two-stage SSEI methods, the condition c = c (16) is supplemented in order to make the following condition be true e ( c − c ) hK g ′ ( k ) e ( c − c ) hK g ′ ( k ) = e ( c − c ) hK g ′ ( k ) e ( c − c ) hK g ′ ( k ) . We remark that condition (16) is not necessary for some special matrix K such as K = 0 or somespecial function g such as scalar functions. The same situation happens in the analysis of Subsection4.4. .3 Vector field F ( ∞ ) For the vector field F ( ∞ ) , if the function f ( y ) := Ky + g ( y ) has the pattern ( u ( y ) , v ( y , y )) ⊺ , thismeans that K and g can be expressed in blocks as K = (cid:18) K K (cid:19) , g ( y ) = (cid:18) g ( y ) g ( y , y ) (cid:19) . (17)Then the following relation is true u ( y ) = K y + g ( y ) , v ( y , y ) = K y + g ( y , y ) . (18) Theorem 4.6
Consider an s -stage SSEI method for solving y ′ = u ( y ) that is of volume preser-vation for the vector field u ( y ) : R m → R m . Let v ( y , y ) : R m + n → R m + n and assume that thereexists an invertible matrix P such that for all y , y ,P ∂ y v ( y , y ) P − = − ∂ y v ( y , y ) ⊺ , P ∂ y g ( y , y ) P − = − ∂ y g ( y , y ) ⊺ . Then the SSEI method is of volume preservation for vector fields f ( y , y ) = ( u ( y ) , v ( y , y )) ⊺ in F ( ∞ ) . Proof From the property of v , we have P K P − = − K ⊺ and (cid:12)(cid:12) e hK (cid:12)(cid:12) = 1. Thus (cid:12)(cid:12) e hK (cid:12)(cid:12) = (cid:12)(cid:12) e hK (cid:12)(cid:12) (cid:12)(cid:12) e hK (cid:12)(cid:12) = (cid:12)(cid:12) e hK (cid:12)(cid:12) . The Jacobian matrix of g ( y ) is block triangular as follows g ′ ( y , y ) = (cid:18) ∂ y g ( y ) 0 ∗ ∂ y g ( y , y ) (cid:19) . In what follows, we prove the condition (10). Using the block transformation, we can bring theleft-hand side of (10) to the block form | I s ⊗ I − h ( A ⊗ I. ∗ E ( hK )) F | = (cid:18) Φ ∗ Φ (cid:19) , where Φ = I − h ¯ a ( hK ) ∂ y g ( k ) · · · − h ¯ a s ( hK ) ∂ y g ( k s )... . . . ... − h ¯ a s ( hK ) ∂ y g ( k ) · · · I − h ¯ a ss ( hK ) ∂ y g ( k s ) , Φ = I − h ¯ a ( hK ) ∂ y g ( k ) · · · − h ¯ a s ( hK ) ∂ y g ( k s )... . . . ... − h ¯ a s ( hK ) ∂ y g ( k ) · · · I − h ¯ a ss ( hK ) ∂ y g ( k s ) . Let F = diag( ∂ y g ( k ) , . . . , ∂ y g ( k s )) and F = diag( ∂ y g ( k ) , . . . , ∂ y g ( k s )) . The above resultcan be simplified as | I s ⊗ I − h ( A ⊗ I. ∗ E ( hK )) F | = | I s ⊗ I − h ( A ⊗ I. ∗ E ( hK )) F | | I s ⊗ I − h ( A ⊗ I. ∗ E ( hK )) F | . Since the SSEI method is of volume preservation for the vector field u ( y ), the following conditionis true | I s ⊗ I − h ( A ⊗ I. ∗ E ( hK )) F | = (cid:12)(cid:12) e hK (cid:12)(cid:12) | I s ⊗ I + h ( A ⊺ ⊗ I. ∗ E ( hK )) F | .
9n the other hand, we compute | I s ⊗ I − h ( A ⊗ I. ∗ E ( hK )) F | = (cid:12)(cid:12) ( I s ⊗ P )( I s ⊗ P − ) − h ( I s ⊗ P )( A ⊗ I. ∗ E ( hK ))( I s ⊗ P − )( I s ⊗ P ) F ( I s ⊗ P − ) (cid:12)(cid:12) = (cid:12)(cid:12) I s ⊗ I + h ( I s ⊗ P )( A ⊗ I. ∗ E ( hK ))( I s ⊗ P − ) F ⊺ (cid:12)(cid:12) = | I s ⊗ I + h ( A ⊗ I. ∗ E ( − hK ⊺ )) F ⊺ | = | I s ⊗ I + hF ( A ⊺ ⊗ I. ∗ E ( − hK ⊺ ) ⊺ ) | (transpose)= | I s ⊗ I + h ( A ⊺ ⊗ I. ∗ E ( − hK ⊺ ) ⊺ ) F | (Sylvester’s law)= | I s ⊗ I + h ( A ⊺ ⊗ I. ∗ E ( hK )) F | (property (14)) . Therefore, the VP condition (10) holds and the SSEI method is of volume preservation for vectorfields in F ( ∞ ) . F (2) Suppose that the function f ( y ) of (1) falls into F (2) . Under this situation, (17) and (18) are stilltrue. We obtain the following result about the VP property of SSEI methods. Theorem 4.7
Consider a one-stage or two-stage SSEI with c = c (or a composition of suchmethod) that is of volume preservation for the vector field u ( y ) : R m → R m . Letting v ( y , y ) : R m + n → R m + n , we assume that there exists an invertible matrix P such that for all y , y ,P ∂ y v ( y , y ) P − = − ∂ y v ( y , y ) , P ∂ y g ( y , y ) P − = − ∂ y g ( y , y ) . Then the SSEI method is of volume preservation for the vector fields f ( y , y ) = ( u ( y ) , v ( y , y )) ⊺ in F (2) . Proof From the conditions of this theorem, it follows that
P K P − = − K and (cid:12)(cid:12) e hK (cid:12)(cid:12) = 1.For the one-stage SSEI, the condition for volume preservation is | I − ha g ′ ( k ) | = | I + ha g ′ ( k ) | , which can be rewritten as | I − ha ∂ y g | | I − ha ∂ y g | = | I + ha ∂ y g | | I + ha ∂ y g | . (19)Since the method is of volume preservation for the vector field u ( y ), we have | I − ha ∂ y g | = | I + ha ∂ y g | . On the other hand, | I − ha ∂ y g | = (cid:12)(cid:12) P P − − ha P ∂ y g P − (cid:12)(cid:12) = | I + ha ∂ y g | . Thus (19) is proved. 10or the two-stage SSEI, it is of volume preservation if and only if (15) is true. According to thespecial result of g ′ , one hasthe left hand side of (15)= | I − ha ∂ y g ( k ) − ha ∂ y g ( k ) + h a a ∂ y g ( k ) ∂ y g ( k ) − h a a ∂ y g ( k ) ∂ y g ( k ) || I − ha ∂ y g ( k ) − ha ∂ y g ( k ) + h a a ∂ y g ( k ) ∂ y g ( k ) − h a a ∂ y g ( k ) ∂ y g ( k ) | = | I + ha ∂ y g ( k ) + ha ∂ y g ( k ) + h a a ∂ y g ( k ) ∂ y g ( k ) − h a a ∂ y g ( k ) ∂ y g ( k ) || I − ha ∂ y g ( k ) − ha ∂ y g ( k ) + h a a ∂ y g ( k ) ∂ y g ( k ) − h a a ∂ y g ( k ) ∂ y g ( k ) | . It then can be verified that | I − ha ∂ y g ( k ) − ha ∂ y g ( k ) + h a a ∂ y g ( k ) ∂ y g ( k ) − h a a ∂ y g ( k ) ∂ y g ( k ) | = | P P − − ha P ∂ y g ( k ) P − − ha P ∂ y g ( k ) P − + h a a P ∂ y g ( k ) P − P ∂ y g ( k ) P − − h a a P ∂ y g ( k ) P − P ∂ y g ( k ) P − | = | I + ha ∂ y g ( k ) + ha ∂ y g ( k ) + h a a ∂ y g ( k ) ∂ y g ( k ) − h a a ∂ y g ( k ) ∂ y g ( k ) | . Consequently, the left hand side of (15)= | I + ha ∂ y g ( k ) + ha ∂ y g ( k ) + h a a ∂ y g ( k ) ∂ y g ( k ) − h a a ∂ y g ( k ) ∂ y g ( k ) || I + ha ∂ y g ( k ) + ha ∂ y g ( k ) + h a a ∂ y g ( k ) ∂ y g ( k ) − h a a ∂ y g ( k ) ∂ y g ( k ) | = the right hand side of (15) . Therefore, all two-stage SSEI methods with c = c are of volume preservation. Remark 4.8
We note that when K = , all the results given in this section reduce to those proposedin [1], which demonstrate the wilder applications of the analysis. In this section, we pay attention to the application of the SSEI methods to various problems andshow the volume preservation of different exponential integrators and ERKN/RKN methods byusing the analysis given in Section 4.
Consider the following first-order systems y ′ ( t ) = J − M y ( t ) + J − ∇ V ( y ( t )) , (20)11here the matrix J is constant and invertible, M is a symmetric matrix and V is a differentiablefunction. This system is the exact pattern (1) with K = J − M, g ( y ( t )) = J − ∇ V ( y ( t )) . (21)It can be verified that Jg ′ ( y ) J − = JJ − ∇ V ( y ) J − = ∇ V ( y ) J − = − g ′ ( y ) ⊺ and J ( K + g ′ ( y )) J − = − ( K + g ′ ( y )) ⊺ . This shows that the set H contains all vector fields of (20) with the same P = J . Thus in the lightof Theorem 4.3, all SSEI methods are of volume preservation for solving the system (20).When J = (cid:18) I − I (cid:19) , the system (20) is a Hamiltonian system y ′ ( t ) = J − ∇ H ( y ( t )) withthe Hamiltonian H ( y ) = 12 y ⊺ M y + V ( y ) . All SSEI methods are of volume preservation for thisHamiltonian system. This is another explanation of the fact that symplectic exponential integratorsare of volume preservation for Hamiltonian systems.Consider another special and important case of (20) by choosing y = (cid:18) qp (cid:19) , J − = (cid:18) I − I N (cid:19) , M = (cid:18) Ω 00 I (cid:19) , V ( y ) = V ( q ) , which gives the following second-order ODE q ′′ − N q ′ + Ω q = −∇ V ( q ) . (22)This system stands for highly oscillatory problems and many problems fall into this kind of equationsuch as the dissipative molecular dynamics, the (damped) Duffing, charged-particle dynamics in aconstant magnetic field and semidiscrete nonlinear wave equations. Applying the SSEI methods to(22) and considering Theorem 4.3, we obtain the following corollary. Corollary 5.1
The following s -stage adapted exponential integrator k i = exp ( c i hK ) q n + exp ( c i hK ) q ′ n − h s P j =1 a ij exp (( c i − c j ) hK ) ∇ V ( k j ) ,i = 1 , , . . . , s,q n +1 = exp ( hK ) q n + exp ( hK ) q ′ n − h s P i =1 b i exp ((1 − c i ) hK ) ∇ V ( k i ) ,q ′ n +1 = exp ( hK ) q n + exp ( hK ) q ′ n − h s P i =1 b i exp ((1 − c i ) hK ) ∇ V ( k i ) (23) are of volume preservation for the second-order highly oscillatory equation (22) , where exp( hK ) ispartitioned into (cid:18) exp ( hK ) exp ( hK )exp ( hK ) exp ( hK ) (cid:19) and the same denotations are used for other matrix-valued functions. Here c = ( c , . . . , c s ) ⊺ , b = ( b , . . . , b s ) ⊺ and A = ( a ij ) s × s are given in Definition .2. If N commutes with Ω , the results of exp ij for i, j = 1 , appearing in (23) can be expressedexplicitly. After some calculations, it is obtained that exp ( hK ) = e h N (cid:16) cosh (cid:0) h √ N − (cid:1) − N sinh (cid:0) h √ N − (cid:1) ( √ N − − (cid:17) , exp ( hK ) = 2 e h N sinh (cid:0) h √ N − (cid:1) ( √ N − − , exp ( hK ) = − Ω exp ( hK ) , exp ( hK ) = e h N (cid:16) cosh (cid:0) h √ N − (cid:1) + N sinh (cid:0) h √ N − (cid:1) ( √ N − − (cid:17) . (24) The results are still true if we replace h by kh with any k ∈ R . If we further assume that Ω = 0, the equation (22) becomes q ′′ = N q ′ − ∇ V ( q ) . (25)One typical example of this system is charged-particle dynamics in a constant magnetic field, whichcan be expressed by (see [11]) x ′′ = x ′ × B + F ( x ) . (26)Here x ( t ) ∈ R describes the position of a particle moving in an electro-magnetic field, F ( x ) = −∇ x U ( x ) is an electric field with the scalar potential U ( x ), and B = ∇ x × A ( x ) is a constantmagnetic field with the vector potential A ( x ) = − x × B .Under the condition that Ω = 0, formulae (24) can be written more succinctly as:exp ( hK ) = I, exp ( hK ) = hϕ ( hN ) , exp ( hK ) = 0 , exp ( hK ) = ϕ ( hN ) , where the ϕ -functions are defined by (see [16, 17]) ϕ ( z ) = e z , ϕ k ( z ) = Z e (1 − σ ) z σ k − ( k − dσ, k = 1 , , . . . . (27)We then get the following volume preserving methods for the special and important second-ordersystem (25). Corollary 5.2
The following s -stage integrator k i = q n + c i hϕ ( c i hN ) q ′ n − h s P j =1 a ij ( c i − c j ) ϕ (( c i − c j ) hN ) ∇ V ( k j ) ,i = 1 , , . . . , s,q n +1 = q n + hϕ ( hN ) q ′ n − h s P i =1 b i (1 − c i ) ϕ ((1 − c i ) hN ) ∇ V ( k i ) ,q ′ n +1 = ϕ ( hN ) q ′ n − h s P i =1 b i ϕ ((1 − c i ) hN ) ∇ V ( k i ) (28) are of volume preservation for the highly oscillatory second-order system (25) , where c = ( c , . . . , c s ) ⊺ , b =( b , . . . , b s ) ⊺ and A = ( a ij ) s × s are given in Definition 2.2. Remark 5.3
It is noted that the above two corollaries are new discoveries which are of great im-portance to Geometric Integration for second-order highly oscillatory problems. .2 Separable partitioned systems It has been proved in [1] that the set S contains all separable partitioned systems. As an example,we consider (cid:18) qp (cid:19) ′ = (cid:18) p − Ω q + ˜ g ( q ) (cid:19) = (cid:18) I − Ω 0 (cid:19) (cid:18) qp (cid:19) + (cid:18) g ( q ) (cid:19) , (29)which is exactly the system (1) with K = (cid:18) I − Ω 0 (cid:19) , g = (cid:18) g ( q ) (cid:19) , f = (cid:18) p − Ω q + ˜ g ( q ) (cid:19) . It can be checked that f and g both fall into S with the same P = diag( I, − I ). For this specialmatrix K , it is clear that e xK = (cid:18) φ ( x Ω) xφ ( x Ω) − x Ω φ ( x Ω) φ ( x Ω) (cid:19) , for x ∈ R (30)with φ i (Ω) := ∞ P l =0 ( − l Ω l (2 l + i )! for i = 0 , . Thus the exponential integrator (3) has a special form,which is called as extended RKN (ERKN) integrators and is represented below.
Definition 5.4 (See [34]) An s -stage ERKN integrator for solving (29) is defined by Q i = φ ( c i V ) q n + hc i φ ( c i V ) p n + h s P j =1 ¯ a ij ( V )˜ g ( Q j ) , i = 1 , . . . , s,q n +1 = φ ( V ) q n + hφ ( V ) p n + h s P i =1 ¯ b i ( V )˜ g ( Q i ) ,p n +1 = − h Ω φ ( V ) q n + φ ( V ) p n + h s P i =1 b i ( V )˜ g ( Q i ) , where c i for i = 1 , . . . , s are real constants, b i ( V ) , ¯ b i ( V ) for i = 1 , . . . , s , and ¯ a ij ( V ) for i, j =1 , . . . , s are matrix-valued functions of V ≡ h Ω . ERKN integrators were firstly proposed in [34]. Further discussions about ERKN integratorshave been given recently, including symmetric integrators (see [30]), symplectic integrators (see[33]), energy-preserving integrators (see [32]) and other kinds integrators (see [26, 28]). However, toour knowledge, the volume-preserving property of ERKN integrators has not been researched yetin the literature. With the analysis given in this paper, we get the following VP result of ERKNintegrators.
Corollary 5.5
Consider a kind of s -stage ERKN integrators Q i = φ ( c i V ) q n + hc i φ ( c i V ) p n + h s P j =1 a ij ( c i − c j ) φ (( c i − c j ) V )˜ g ( Q j ) ,i = 1 , . . . , s,q n +1 = φ ( V ) q n + hφ ( V ) p n + h s P i =1 b i (1 − c i ) φ ((1 − c i ) V )˜ g ( Q i ) ,p n +1 = − h Ω φ ( V ) q n + φ ( V ) p n + h s P i =1 b i φ ((1 − c i ) V )˜ g ( Q i ) , (31)14 here c = ( c , . . . , c s ) ⊺ , b = ( b , . . . , b s ) ⊺ and A = ( a ij ) s × s are given in Definition 2.2. Underthe condition that b j = 0 for all j = 1 , . . . , s, all one-stage and two-stage (with c = c ) ERKNintegrators (31) , and compositions thereof, are of volume preservation for solving the separablepartitioned system (29) . Proof In the light of Definition 2.2 and the result (30), we adapt the SSEI methods to the system(29) and then get the scheme (31). Based on Theorem 4.4, the volume preserving result of (31) isimmediately obtained.
Remark 5.6
We note that this is a novel result which studies the volume-preserving ERKN inte-grators for (29) . Moreover, from the scheme (31) , it can be observed that all one-stage and two-stagewith c = c ERKN integrators are explicit, which means that we obtain explicit volume preservingERKN integrators for the separable partitioned system (29) . If Ω is a symmetric and positive semi-definite matrix, and ˜ g ( q ) = −∇ U ( q ), the system (29) isan oscillatory Hamiltonian system (cid:18) qp (cid:19) ′ = (cid:18) I − I (cid:19) ∇ H ( q, p ) with the Hamiltonian H ( q, p ) = 12 p ⊺ p + 12 q ⊺ Ω q + U ( q ) . (32)It has been noted in Subsection 5.1 that this vector field falls into the set H . Thus Theorem 4.3provides another way to prove the well-known fact that all symplectic ERKN integrators (31) areof volume preservation for the oscillatory Hamiltonian system (32).In what follows, we study the volume-preserving property of RKN methods for second-orderODEs. Consider Ω = 0 for the above analysis and under this situation, ERKN integrators reduceto RKN methods. Therefore, we are now in a position to present the following volume-preservingproperty for RKN methods. Corollary 5.7
Consider the following s -stage RKN methods Q i = q n + hc i q ′ n + h s P j =1 a ij ( c i − c j )˜ g ( Q j ) , i = 1 , . . . , s,q n +1 = q n + hq ′ n + h s P i =1 b i (1 − c i )˜ g ( Q i ) ,q ′ n +1 = q ′ n + h s P i =1 b i ˜ g ( Q i ) (33) with the coefficients c = ( c , . . . , c s ) ⊺ , b = ( b , . . . , b s ) ⊺ and A = ( a ij ) s × s of an s -stage RK method.If this RK method is symplectic and b j = 0 for all j = 1 , . . . , s, all one-stage and two-stage RKNmethods (33) (and compositions thereof ) are of volume preservation for solving the second-ordersystem q ′′ = ˜ g ( q ) . Remark 5.8
It is noted that the fact of this corollary can be derived in a different way. Hairer,Lubich and Wanner have proved in [13] that any symplectic RK method with at most two stages (andany composition of such methods) is volume preserving for separable divergence free vector fields.Rewriting the second-order equation q ′′ = ˜ g ( q ) as a first-order system and applying symplectic RKmethods to it implies the result of Corollary 5.7. It is seen that the analysis of volume-preservingERKN integrators provides an alternative derivation of the volume preservation of RKN methods. .3 Other applications It has been shown in [1] that F ( ∞ ) contains the affine vector fields f ( y ) = Ly + d such that (cid:12)(cid:12) I + h L (cid:12)(cid:12) = (cid:12)(cid:12) I − h L (cid:12)(cid:12) for all h >
0. For solving the system in the affine vector fields, the exponentialintegrator (3) becomes k i = e c i hL y n + h s P j =1 ¯ a ij ( hL ) d, i = 1 , , . . . , s,y n +1 = e hL y n + h s P i =1 ¯ b i ( hL ) d. In the light of Theorem 4.6, this SSEI method is of volume preservation for the affine vector fields.It was also noted in [1] that F ( ∞ ) contains the vector fields f ( y ) such that f ′ ( y ) = JS ( y ) witha skew-symmetric matrix J and the symmetric matrix S ( y ). Assume that K = JM, g ′ ( y ) = JS ( y ) , (34)where M is a symmetric matrix. The system (1) with the vector field (34) falls into F ( ∞ ) . Thusall SSEI methods are of volume preservation for the vector field (34). In order to show the performance of SSEI methods, the solvers chosen for comparison are: • SSRK1: the Gauss-Legendre method of order two; • SSEI1: the one-stage SSEI method with the coefficients (4), where the RK method (5) ischosen as SSRK1; • SSRK2: the Gauss-Legendre method of order four; • SSEI2: the two-stage SSEI method with the coefficients (4), where the RK method (5) ischosen as SSRK2.It is noted that all these methods are general implicit and we use fixed-point iteration here. We set10 − as the error tolerance and 100 as the maximum number of each fixed-point iteration. Problem 1.
As the first numerical example, we consider the Duffing equation defined by (cid:18) qp (cid:19) ′ = (cid:18) − ω − k (cid:19) (cid:18) qp (cid:19) + (cid:18) k q (cid:19) , (cid:18) q (0) p (0) (cid:19) = (cid:18) ω (cid:19) . The exact solution of this system is q ( t ) = sn ( ωt ; k/ω ) with the Jacobi elliptic function sn . Sinceit is a Hamiltonian system, all the methods chosen for comparison are volume preserving. For thisproblem, we choose k = 0 .
07 and ω = 20 and then solve it on the interval [0 , h = 1 / , / , / , / { i } i =1 ,..., of thefour methods are given in Figure 2. From the results, it can be observed clearly that the integratorsSSEI1 and SSEI2 perform better than Runge-Kutta methods in the preservation of volume. Finally,we integrate this problem in [0 , t end ] with h = 0 . / i for i = 1 , . . . , . The relative global errors16 p The flow of SSEI1 with h=1/2 −1 0 1−20−1001020 q p The flow of SSEI1 with h=1/10 −1 0 1−20−1001020 q p The flow of SSEI1 with h=1/50 −1 −0.5 0 0.5 1−20−1001020 q p The flow of SSEI1 with h=1/200 −1 0 1−20−1001020 q p The flow of SSRK1 with h=1/2 −100 0 100 200−2000−1000010002000 q p The flow of SSRK1 with h=1/10 −1 0 1−20−1001020 q p The flow of SSRK1 with h=1/50 −1 −0.5 0 0.5 1−20−1001020 q p The flow of SSRK1 with h=1/200 −1 0 1−20−1001020 q p The flow of SSEI2 with h=1/2 −1 0 1−20−1001020 q p The flow of SSEI2 with h=1/10 −1 0 1−20−1001020 q p The flow of SSEI2 with h=1/50 −1 −0.5 0 0.5 1−20−1001020 q p The flow of SSEI2 with h=1/200 −1 0 1−20−1001020 q p The flow of SSRK2 with h=1/2 −1 0 1−20−1001020 q p The flow of SSRK2 with h=1/10 −1 0 1−20−1001020 q p The flow of SSRK2 with h=1/50 −1 −0.5 0 0.5 1−20−1001020 q p The flow of SSRK2 with h=1/200
Figure 2: Problem 1: the flows of different methods. (t/h) l og ( G E ) The global errors with t end =10
SSEI1SSRK1SSEI2SSRK2 3 3.5 4 4.5−12−10−8−6−4−202 log (t/h) l og ( G E ) The global errors with t end =100
SSEI1SSRK1SSEI2SSRK2 4 4.5 5 5.5−10−8−6−4−202 log (t/h) l og ( G E ) The global errors with t end =1000
SSEI1SSRK1SSEI2SSRK2
Figure 3: Problem 1: the relative global errors.17or different t end are presented in Figure 3. These results show clearly again that exponentialintegrators have better accuracy than Runge-Kutta methods. Problem 2.
Consider the following divergence free ODEs xyz ′ = − ω ω − ω ω xyz + sin( x − z )0sin( x − z ) . By choosing P = , it can be checked that the vector field of this problem falls into S . We consider ω = 100 and (0 . , . , . ⊺ as the initial value. This problem is firstly integratedon [0 , h = 1 / , / , / , /
400 and the numerical flows x and y at the time points { i } i =1 ,..., are shown in Figure 4. Then the relative global errors for different t end with h = 0 . / i for i = 2 , . . . , Problem 3.
Consider the damped Helmholtz-Duffing oscillator (see [6]) q ′′ + 2 υq ′ + Aq = − Bq − εq , where q denotes the displacement of the system, A is the natural frequency, ε is a non-linearsystem parameter, υ is the damping factor, and B is a system parameter independent of time. It iswell known that the dynamical behavior of eardrum oscillations, elasto-magnetic suspensions, thinlaminated plates, graded beams, and other physical phenomena all fall into this kind of equations.We choose the parameters υ = 0 . , A = 200 , B = − . , ε = 1and the initial value q (0) = 1 and q ′ (0) = 15 . , h = 1 / , / , / , / q and p = q ′ at the time points { i } i =1 ,..., in Figure 6. We then solve the problem with different t end = 10 , , h = 0 . / i for i = 0 , . . . ,
3. The relative global errors are shown in Figure 7. It follows again fromthe results that SSEI methods perform much better than SSRK methods.
Problem 4.
The last numerical experiment is concerned with the charged particle system witha constant magnetic field (see [11]). The system can be given by (26) with the potential U ( x ) = √ x + x and the constant magnetic field B = (0 , , ⊺ . The initial values are chosen as x (0) =(0 . , , . ⊺ and x ′ (0) = (0 . , . , . ⊺ . We firstly integrate this system on [0 , h =1 / , / , / , /
200 and show the numerical flows x and v = x ′ at the time points { i } i =1 ,..., in Figure 8. Then the problem is solved with t end = 10 , , h = 0 . / i for i = 0 , . . . , y The flow of SSEI1 with h=1/50 0 0.5 1−0.500.5 x y The flow of SSEI1 with h=1/100 0 0.5 1−0.500.5 x y The flow of SSEI1 with h=1/200 0 0.5 1−0.500.5 x y The flow of SSEI1 with h=1/4000 0.5 100.51 x y The flow of SSRK1 with h=1/50 0 0.5 1−0.500.5 x y The flow of SSRK1 with h=1/100 0 0.5 1−0.500.5 x y The flow of SSRK1 with h=1/200 0 0.5 1−0.500.5 x y The flow of SSRK1 with h=1/4000 0.5 1−0.500.5 x y The flow of SSEI2 with h=1/50 0 0.5 1−0.500.5 x y The flow of SSEI2 with h=1/100 0 0.5 1−0.500.5 x y The flow of SSEI2 with h=1/200 0 0.5 1−0.500.5 x y The flow of SSEI2 with h=1/4000 0.5 1−0.500.5 x y The flow of SSRK2 with h=1/50 0 0.5 1−0.500.5 x y The flow of SSRK2 with h=1/100 0 0.5 1−0.500.5 x y The flow of SSRK2 with h=1/200 0 0.5 1−0.500.5 x y The flow of SSRK2 with h=1/400
Figure 4: Problem 2: the flows of different methods. (t/h) l og ( G E ) The global errors with t end =10
SSEI1SSRK1SSEI2SSRK2 3.5 4 4.5 5−7−6−5−4−3−2−101 log (t/h) l og ( G E ) The global errors with t end =100
SSEI1SSRK1SSEI2SSRK2 4.5 5 5.5 6−8−7−6−5−4−3−2−101 log (t/h) l og ( G E ) The global errors with t end =1000
SSEI1SSRK1SSEI2SSRK2
Figure 5: Problem 2: the relative global errors.19
20 0 20−101 p q The flow of SSEI1 with h=1/2 −20 0 20−101 p q The flow of SSEI1 with h=1/10 −20 0 20−101 p q The flow of SSEI1 with h=1/50 −20 0 20−101 p q The flow of SSEI1 with h=1/200−20 0 20−101 p q The flow of SSRK1 with h=1/2 −20 0 20−101 p q The flow of SSRK1 with h=1/10 −20 0 20−101 p q The flow of SSRK1 with h=1/50 −20 0 20−101 p q The flow of SSRK1 with h=1/200−20 0 20−101 p q The flow of SSEI2 with h=1/2 −20 0 20−101 p q The flow of SSEI2 with h=1/10 −20 0 20−101 p q The flow of SSEI2 with h=1/50 −20 0 20−101 p q The flow of SSEI2 with h=1/200−20 0 20−101 p q The flow of SSRK2 with h=1/2 −20 0 20−101 p q The flow of SSRK2 with h=1/10 −20 0 20−101 p q The flow of SSRK2 with h=1/50 −20 0 20−101 p q The flow of SSRK2 with h=1/200
Figure 6: Problem 3: the flows of different methods. (t/h) l og ( G E ) The global errors with t end =10
SSEI1SSRK1SSEI2SSRK2 3 3.5 4−7−6−5−4−3−2−101 log (t/h) l og ( G E ) The global errors with t end =100
SSEI1SSRK1SSEI2SSRK2 4 4.5 5−7−6−5−4−3−2−1012 log (t/h) l og ( G E ) The global errors with t end =1000
SSEI1SSRK1SSEI2SSRK2
Figure 7: Problem 3: the relative global errors.20 x The flow of SSEI1 with h=1/2 −1 0 10.80.91 v x The flow of SSEI1 with h=1/10 −1 0 10.80.91 v x The flow of SSEI1 with h=1/50 −1 0 10.80.91 v x The flow of SSEI1 with h=1/200−1 0 10.80.91 v x The flow of SSRK1 with h=1/2 −1 0 10.80.91 v x The flow of SSRK1 with h=1/10 −1 0 10.80.91 v x The flow of SSRK1 with h=1/50 −1 0 10.80.91 v x The flow of SSRK1 with h=1/200−1 0 10.80.91 v x The flow of SSEI2 with h=1/2 −1 0 10.80.91 v x The flow of SSEI2 with h=1/10 −1 0 10.80.91 v x The flow of SSEI2 with h=1/50 −1 0 10.80.91 v x The flow of SSEI2 with h=1/200−1 0 10.80.91 v x The flow of SSRK2 with h=1/2 −1 0 10.80.91 v x The flow of SSRK2 with h=1/10 −1 0 10.80.91 v x The flow of SSRK2 with h=1/50 −1 0 10.80.91 v x The flow of SSRK2 with h=1/200
Figure 8: Problem 4: the flows of different methods. (t/h) l og ( G E ) The global errors with t end =10
SSEI1SSRK1SSEI2SSRK2 3 3.5 4−12−10−8−6−4−20 log (t/h) l og ( G E ) The global errors with t end =100
SSEI1SSRK1SSEI2SSRK2 4 4.5 5−11−10−9−8−7−6−5−4−3−2 log (t/h) l og ( G E ) The global errors with t end =1000
SSEI1SSRK1SSEI2SSRK2
Figure 9: Problem 4: the relative global errors.21
Conclusions
This paper studied volume-preserving exponential integrators. The necessary and sufficient volume-preserving condition for exponential integrators was derived and volume-preserving properties werediscussed for four kinds of vector fields. It was shown that symplectic exponential integrators can beof volume preservation for a much larger class of vector fields than Hamiltonian systems. It should benoted that a new result has been proved that all a kind of adapted exponential integrators methodsis of volume preservation for the highly oscillatory system (22). Moreover, the volume-preservingproperty of ERKN/RKN methods is discussed for separable partitioned systems. Some new resultson Geometric Integration were presented for second-order highly oscillatory problems and separablepartitioned systems. We also carried out four numerical experiments to demonstrate the remarkablerobustness and superiority of volume-preserving exponential integrators in comparison with volume-preserving Runge-Kutta methods.
References [1]
P. Badera, D. I. McLarena, G.R.W. Quispela, and M. Webb , Volume preservation byRunge–Kutta methods , Appl. Numer. Math., 109 (2016), pp. 123-137.[2]
L. Brugnano, G. Frasca Caccia, and F. Iavernaro , Hamiltonian Boundary ValueMethods (HBVMs) and their efficient implementation , mathematics in engineering, science andaerospace mesa, 5 (2014), pp. 343-411.[3]
E. Celledoni, R.I. McLachlan, D.I. McLaren, B. Owren, G.R.W. Quispel, andW.M. Wright , Energy-preserving Runge–Kutta methods , M2AN Math. Model. Numer. Anal.,43 (2009), pp. 645-649.[4]
P. Chartier and A. Murua , Preserving first integrals and volume forms of additively splitsystems , IMA J. Numer. Anal., 27 (2007), pp. 381-405.[5]
D. Cohen and E. Hairer , Linear energy-preserving integrators for Poisson systems , BITNumer. Math., 51 (2011), pp. 91-101.[6]
A. El´ıas-Z´u˜niga , Analytical solution of the damped Helmholtz-Duffing equation , Appl. Math.Lett. 25, (2012), pp. 2349-2353.[7]
K. Feng and M. Qin , Symplectic Geometric algorithms for Hamiltonian systems , Springer-Verlag, Berlin, Heidelberg, 2010.[8]
K. Feng and Z.J. Shang , Volume-preserving algorithms for source-free dynamical systems ,Numer. Math., 71 (1995), pp. 451-463.[9]
E. Hairer , Energy-preserving variant of collocation methods , J. Numer. Anal. Ind. Appl. Math.,5 (2010), pp. 73-84.[10]
E. Hairer and C. Lubich , Long-time energy conservation of numerical methods for oscilla-tory differential equations , SIAM J. Numer. Anal., 38 (2000), pp. 414-441.[11]
E. Hairer and C. Lubich , Symmetric multistep methods for charged-particle dynamics ,SMAI J. Comput. Math., 3 (2017), pp. 205-218.2212]
E. Hairer and C. Lubich , Long-term analysis of the St¨ormer-Verlet method for Hamiltoniansystems with a solution-dependent high frequency , Numer. Math., 134 (2016), pp. 119-138.[13]
E. Hairer, C. Lubich, and G. Wanner , Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations , 2nd edn. Springer-Verlag, Berlin,Heidelberg, 2006.[14]
Y. He, Y. Sun, J. Liu, and H. Qin , Volume-preserving algorithms for charged particledynamics , J. Comput. Phys., 281 (2015), pp. 135-147.[15]
M. Hochbruck and A. Ostermann , Explicit exponential Runge–Kutta methods for semi-lineal parabolic problems , SIAM J. Numer. Anal., 43 (2005), pp. 1069-1090.[16]
M. Hochbruck and A. Ostermann , Exponential integrators , Acta Numer., 19 (2010), pp.209-286.[17]
M. Hochbruck, A. Ostermann, and J. Schweitzer , Exponential rosenbrock-type meth-ods , SIAM J. Numer. Anal., 47 (2009), pp. 786-803.[18]
A. Iserles, G.R.W. Quispel, and P.S.P. Tse , B-series methods cannot be volume-preserving , BIT Numer. Math., 47 (2007), pp. 351-378.[19]
Y.W. Li and X. Wu , Exponential integrators preserving first integrals or Lyapunov functionsfor conservative or dissipative systems , SIAM J. Sci. Comput., 38 (2016), pp. 1876-1895.[20]
R.I. McLachlan, H.Z. Munthe-Kaas, G.R.W. Quispel, and A. Zanna , Explicit volume-preserving splitting methods for linear and quadratic divergence-free vector fields , Found. Comput.Math., 8 (2008), pp. 335-355.[21]
R.I. McLachlan and G. R. W. Quispel , Discrete gradient methods have an energy con-servation law , Discrete Contin. Dyn. Syst., 34 (2014), pp. 1099-1104.[22]
R.I. McLachlan and C. Scovel , A survey of open problems in symplectic integration , FieldsInst. Commun., 10 (1998), pp. 151-180.[23]
L. Mei and X. Wu , Symplectic exponential Runge–Kutta methods for solving nonlinear Hamil-tonian systems , J. Comput. Phys., 338 (2017), pp. 567-584.[24]
G.R.W. Quispel , Volume-preserving integrators , Phys. Lett. A, 206 (1995), pp. 26-30.[25]
J.M. Sanz-Serna , Symplectic integrators for Hamiltonian problems: An overview , in ActaNumerica 1992, A. Iserles, ed., Cambridge University Press, Cambridge, UK, (1992), pp. 243-286.[26]
B. Wang, A. Iserles, and X. Wu , Arbitrary-order trigonometric Fourier collocation meth-ods for multi-frequency oscillatory systems , Found. Comput. Math., 16 (2016), pp. 151-181.[27]
B. Wang and X. Wu , Functionally-fitted energy-preserving integrators for Poisson systems ,J. Comput. Phys., 364 (2018), pp. 137-152.[28]
B. Wang, X. Wu, and F. Meng , Trigonometric collocation methods based on Lagrange basispolynomials for multi-frequency oscillatory second-order differential equations , J. Comput. Appl.Math., 313 (2017), pp. 185-201. 2329]
B. Wang, X. Wu, F. Meng, and Y. Fang , Exponential Fourier collocation methods forsolving first-order differential equations , J. Comput. Math., 35 (2017), pp. 711-736.[30]
B. Wang, H. Yang, and F. Meng , Sixth order symplectic and symmetric explicit ERKNschemes for solving multi-frequency oscillatory nonlinear Hamiltonian equations , Calcolo, 54(2017), pp. 117-140.[31]
X. Wu and B. Wang , Recent Developments in Structure-Preserving Algorithms for Oscilla-tory Differential Equations , Springer Nature Singapore Pte Ltd, 2018.[32]
X. Wu, B. Wang, and W. Shi , Efficient energy preserving integrators for oscillatory Hamil-tonian systems , J. Comput. Phys., 235 (2013), pp. 587-605.[33]
X. Wu, X. You, and B. Wang , Structure-preserving algorithms for oscillatory differentialequations , Springer-Verlag, Berlin, Heidelberg, 2013.[34]
X. Wu, X. You, W. Shi, and B. Wang , ERKN integrators for systems of oscillatorysecond-order differential equations , Comput. Phys. Comm., 181 (2010), pp. 1873-1887.[35]
H. Xue and A. Zanna , Explicit volume-preserving splitting methods for polynomialdivergence-free vector fields , BIT Numer. Math., 53 (2013), pp. 265-281.[36]