GGEOMETRIC EXPONENTIAL INTEGRATORS
XUEFENG SHEN AND MELVIN LEOK
Abstract.
In this paper, we consider exponential integrators for semilinearPoisson systems. Two types of exponential integrators are constructed, onepreserves the Poisson structure, and the other preserves energy. Numerical ex-periments for semilinear Possion systems obtained by semi-discretizing Hamil-tonian PDEs are presented. These geometric exponential integrators exhibitbetter long time stability properties as compared to non-geometric integrators,and are computationally more efficient than traditional symplectic integratorsand energy-preserving methods based on the discrete gradient method. Introduction
Exponential integrators [4] are a class of numerical integrators for stiff systemswhose vector field can be decomposed into a linear term and a nonlinear term,(1) ˙ q = Aq + f ( q ) . Usually, the coefficient matrix A has a large spectral radius, and is responsible forthe stiffness of the system of differential equations, while the nonlinear term f ( q ) isrelatively smooth. There are various ways to construct an exponential integrator [7].For example, we can perform a change of variables ˜ q ( t ) = e − At q ( t ), and transform(1) to obtain(2) [ e − At q ( t )] (cid:48) = ˜ q (cid:48) ( t ) = e − At f ( e At ˜ q ( t )) . Notice that the Jacobi matrix of (2) equals e − At ∇ f e At , which has a smaller spectralradius than the Jacobi matrix A + ∇ f of (1). A natural idea is to apply a classicalintegrator for the mollified system (2) to obtain an approximation of ˜ q ( t ), theninvert the change of variables to obtain an approximation of the solution q ( t ) of(1). In Section 2, we shall demonstrate how to construct symplectic exponentialintegrators using this approach.Another way of constructing exponential integrators starts from the variation-of-constants formula,(3) q ( t ) = e A ( t − t ) q ( t ) + (cid:90) tt e A ( t − τ ) f ( q ( τ )) dτ, which is the exact solution for (1) with initial condition q ( t ) = q . Then, acomputable approximation can be obtained by approximating the f ( q ( τ )) terminside the integral. If we approximate f ( q ( τ )) by f ( q k ), we arrive at the exponentialEuler method,(4) q k +1 = e Ah q k + (cid:90) h e Aτ dτ · f ( q k ) . a r X i v : . [ m a t h . NA ] M a r XUEFENG SHEN AND MELVIN LEOK
An exponential Runge–Kutta method of collocation type [3] could also be con-structed by approximating f ( q ( τ )) with polynomials. In Section 3, we shall showhow to construct energy-preserving exponential integrators from (3).In this paper, we consider a specific form of (1) which is a Poisson system. Weassume A = JD , f ( q ) = J ∇ V ( q ), where J T = − J, D T = D , and JD = DJ , thusthe coefficient matrix A is also skew-symmetric. Now, the semilinear system (1)can be written as,(5) ˙ q = J ( Dq + ∇ V ( q )) = J ∇ H ( q ) , with Hamiltonian function H ( q ) = q T Dq + V ( q ). Equation (5) describes a con-stant Poisson system, and there are at least two quantities that are preserved bythe flow: the Poisson structure J ij ∂∂x i ⊗ ∂∂x j and Hamiltonian H ( q ). Geometricintegrators that preserve the geometric structure and first integrals of the systemtypically exhibit superior qualitative properties when compared to non-geometricintegrators, and they are an active area of research [2, 5, 6].Here, we construct geometric exponential integrators that either preserve thePoisson structure or Hamiltonian of (5). They exhibit long time stability, allowfor relatively larger timesteps for the stiff problem, and are computationally moreefficient as they can be implemented using fixed point iterations as opposed to New-ton type iterations. For the rest of the paper, Section 2 is devoted to developingsymplectic exponential integrators that preserve the Poisson structure; Section 3 isdevoted to developing energy preserving exponential integrators; numerical meth-ods and experiments are presented in Section 4 and Section 5, respectively.2. Symplectic Exponential Integrator
For constant Poisson systems (5), it was shown in [12] that the midpoint rule anddiagonally implicit symplectic Runge–Kutta methods preserve the Poisson structure J ij ∂∂x i ⊗ ∂∂x j . We first start by constructing an exponential midpoint rule: applythe classical midpoint rule to the transformed system (2) to obtain(6) ˜ q k +1 − ˜ q k h = e − At k +1 / f (cid:16) e At k +1 / ˜ q k +1 + ˜ q k (cid:17) , where t k +1 / = t k + t k +1 , h = t k +1 − t k , ˜ q k = e − At k q k , ˜ q k +1 = e − At k +1 q k +1 . Transform (6) back to q k and q k +1 , and we obtain the exponential midpoint rule,(7) q k +1 = e Ah q k + h · e A h f (cid:16) e A h q k + e − A h q k +1 (cid:17) . Theorem 1.
The exponential midpoint rule (7) preserves the Poisson structurewhen applied to the semilinear Poisson system (5) .Proof.
Recall that a map φ preserves Poisson structure J ij ∂∂x i ⊗ ∂∂x j iff(8) ( ∇ φ ) J ( ∇ φ ) T = J. EOMETRIC EXPONENTIAL INTEGRATORS 3
Differentiating (7), we obtain, dq k +1 = e Ah dq k + h · e A h ∇ f (cid:16) e A h dq k + 12 e − A h dq k +1 (cid:17) , (cid:16) I − h · e A h ∇ f e − A h (cid:17) dq k +1 = (cid:16) e Ah + h · e A h ∇ f e A h (cid:17) dq k . So the map φ ( q k ) = q k +1 has Jacobi matrix ∇ φ = M − N , where M = I − h · e A h ∇ f e − A h = I − h · e A h J ∇ V e − A h ,N = e Ah + h · e A h ∇ f e A h = e Ah + h · e A h J ∇ V e A h . Then, we just need to verify (8), which is equivalent to
M JM T = N JN T ,(9) M JM T = (cid:16) I − h · e A h J ∇ V e − A h (cid:17) J (cid:16) I + h · e A h ∇ V Je − A h (cid:17) = J − h e A h J ∇ V e − A h Je A h ∇ V Je − A h = J − h e A h J ∇ V J ∇ V Je − A h = N JN T . In (9), we used the property that ∇ V is symmetric, that A is skew-symmetricwhich implies that ( e Ah ) T = e − Ah , and the assumptions that JD = DJ and that e A h and J commute. (cid:3) The exponential midpoint rule is a second-order method, and we will now develophigher-order symplectic exponential integrators. Recall that a general diagonallyimplicit Symplectic Runge–Kutta method (DISRK) has a Butcher tableau of thefollowing form,
Table 1.
DISRK c b c b b c b b b ... ... ... ... . . . ...c s b b b · · · b s b b b · · · b s It was shown in [12] that any DISRK method can be regarded as the composi-tion of midpoint rules with timesteps b h, b h, b h, · · · b s h . If we apply the DISRKmethod for the transformed system (2), and then convert back, we obtain the fol-lowing diagonally implicit symplectic exponential (DISEX) integrator,(10) Q i = e Ahc i q k + h · i (cid:88) j =1 a ij e Ah ( c i − c j ) f ( Q j ) ,q k +1 = e Ah q k + h · s (cid:88) i =1 b i e Ah (1 − c i ) f ( Q i ) . XUEFENG SHEN AND MELVIN LEOK where a ij are the coefficients in Table 1. This integrator can be represented interms of the following Butcher tableau, Table 2.
DISEX e Ahc b e Ahc b e Ah ( c − c ) b e Ahc b e Ah ( c − c ) b e Ah ( c − c ) b ... ... ... ... . . . ...e Ahc s b e Ah ( c s − c ) b e Ah ( c s − c ) b e Ah ( c s − c ) · · · b s e Ah b e Ah (1 − c ) b e Ah (1 − c ) b e Ah (1 − c ) · · · b s e Ah (1 − c s ) The DISEX method preserves the Poisson structure, as demonstrated by thefollowing theorem.
Theorem 2.
The DISEX method is equivalent to the composition of exponentialmidpoint rules with timesteps b h, b h, b h, · · · b s h .Proof. The composition of exponential midpoint rules q k b h −−→ Z b h −−→ Z · · · b s h −−→ Z s = q k +1 is represented as follows, Z = e Ahb q k + b h · e Ah b f (cid:16) e Ah b q k + e − Ah b Z (cid:17) , (M.1) Z = e Ahb Z + b h · e Ah b f (cid:16) e Ah b Z + e − Ah b Z (cid:17) , (M.2) ...Z i = e Ahb i Z i − + b i h · e Ah bi f (cid:16) e Ah bi Z i − + e − Ah bi Z i (cid:17) , (M.i) ...Z s = e Ahb s Z s − + b s h · e Ah bs f (cid:16) e Ah bs Z s − + e − Ah bs Z s (cid:17) . (M.s)Introduce Q = e Ah b q k + e − Ah b Z in (M.1) as an intermediate variable. Then, onboth sides of (M.1), multiply by e − Ah b , add e Ah b q k , and divide by two, whichyields an equivalent form of (M.1),(S.1) Q = e Ah b q k + b h · f ( Q ) . EOMETRIC EXPONENTIAL INTEGRATORS 5
For the Runge–Kutta method represented by the Butcher tableau in Table 1 to beconsistent, the coefficients have to satisfy, c = b ,c = b + b ,...c i = b + b + · · · b i − + b i ,...c s = b + b + b + · · · + b s − + b s , b + b + b + · · · + b s − + b s . So equation (S.2) coincides with the first line of the Butcher tableau of the DI-SEX method. Similarly, introduce Q = e Ah b Z + e − Ah b Z . Then, on both sidesof (M.2), multiply by e − Ah b , add e Ah b Z , then divided by two, which yields anequivalent form of (M.2):(S.2) Q = e Ah b Z + b h · f ( Q )= e Ah ( b + b ) q k + b h · e Ah ( b + b ) f ( Q ) + b h · f ( Q )= e Ahc q k + b h · e Ah ( c − c ) f ( Q ) + b h · f ( Q ) . So equation (S.2) coincides with the second line of the Butcher tableau of the DISEXmethod. Then, as before, we introduce Q i = e Ah bi Z i − + e − Ah bi Z i , and apply thesame technique to (M.i), which yields Q i = e Ah bi Z i − + b i h · f ( Q i ) . By induction, Z i − = e Ah ( b i − + ··· + b + b ) q k + b h · e Ah ( b i − + b i − + ··· + b ) f ( Q )+ b h · e Ah ( b i − + b i − + ··· + b ) f ( Q ) + · · · b i − h · e Ah ( bi − ) f ( Q i − ) , so(S.i) Q i = e Ah ( bi + b i − + ··· + b + b ) q k + b h · e Ah ( bi + b i − + b i − + ··· + b ) f ( Q )+ b h · e Ah ( bi + b i − + b i − + ··· + b ) f ( Q ) + · · · + b i − h · e Ah ( bi + bi − ) f ( Q i − ) + b i h · f ( Q i )= e Ahc i q k + b h · e Ah ( c i − c ) f ( Q ) + · · · + b i − h · e Ah ( c i − c i − ) f ( Q i − ) + b i h · f ( Q i ) , XUEFENG SHEN AND MELVIN LEOK which coincides with the i -th row of the Butcher tableau of the DISEX method.Finally, we have q k +1 = Z s = e Ahb s Z s − + b s h · f ( Q s )= e Ah ( b s + ··· b + b ) q k + b h · e Ah ( b s + ··· b + b ) f ( Q ) + · · · + b s − h · e Ah ( b s + bs − ) f ( Q s − ) + b s h · e Ah bs f ( Q s )= e Ah q k + b h · e Ah (1 − c ) f ( Q ) + · · · b s h · e Ah (1 − c s ) f ( Q s ) , which coincides with the last row of the Butcher tableau of the DISEX method. Sothe composition of exponential midpoint rules with timesteps b h, b h, b h, · · · b s h is equivalent to the DISEX method of Table 2. (cid:3) Energy-preserving Exponential Integrator
Though classical symplectic methods exhibit superior long time stability, it wasobserved that symplectic schemes are less competitive for the numerical integrationof stiff systems with high frequency. In sharp contrast, energy-preserving methodsperform much better [9]. A general way to construct an energy-preserving methodfor a Poisson system ˙ q = J ∇ H ( q ) is the discrete gradient method [8]. We design adiscrete gradient ∇ H ( q k , q k +1 ) that satisfies the following property,(11) ∇ H ( q k , q k +1 ) · ( q k +1 − q k ) = H ( q k +1 ) − H ( q k ) . Then, the resulting discrete gradient method is given by,(12) q k +1 − q k h = J ∇ H ( q k , q k +1 ) . Multiplying ∇ H ( q k , q k +1 ) on both sides of (12), we obtain(13) H ( q k +1 ) − H ( q k ) = ∇ H ( q k , q k +1 ) · ( q k +1 − q k )= h · ∇ H ( q k , q k +1 ) J ∇ H ( q k , q k +1 )= 0 . The last equation of (13) holds simply due to the skew-symmetric property of matrix J , which implies that discrete gradient method (12) preserves energy. We shall com-bine exponential integrators with the discrete gradient method to obtain an energy-preserving exponential integrator. This approach was initially proposed in [11] forseparable Hamiltonian systems using the extended discrete gradient method, andwe generalize this to semilinear Poisson systems. Replace the f ( q k ) term in theexponential Euler method (4) by the discrete gradient J ∇ V ( q k , q k +1 ), which yields(14) q k +1 = e Ah q k + (cid:90) h e Aτ dτ · J ∇ V ( q k , q k +1 ) . Theorem 3.
Method (14) preserves the Hamiltonian H ( q ) .Proof. Let S = e Ah , T = (cid:82) h e Aτ dτ , then q k +1 = Sq k + T J ∇ V ( q k , q k +1 ), and wewill show that the following properties hold:(1) S T = S − ;(2) AT = S − I ;(3) AT T = I − S T ; EOMETRIC EXPONENTIAL INTEGRATORS 7 (4) S T T = T T .Property 1 follows from the fact that S T = ( e Ah ) T = e − Ah = S − . Property 2follows from e Ah − I = e Aτ (cid:12)(cid:12)(cid:12)(cid:12) h = (cid:90) h A · e Aτ dτ = AT.
Taking transposes on both sides of Property 2 and using the fact that A and T commute gives Property 3. Property 4 follow from S T T = e − Ah (cid:90) h e Aτ dτ = (cid:90) h e − A ( h − τ ) dτ = (cid:90) h e − Aτ dτ = T T . Recall that the Hamiltonian of the constant Poisson system (5) is H ( q ) = q T Dq + V ( q ), so(15) H ( q k +1 ) = 12 q T k +1 Dq k +1 + V ( q k +1 )= 12 ( Sq k + T J ∇ V ) T D ( Sq k + T J ∇ V ) + V ( q k +1 )= 12 q T k S T DSq k + q T k S T DT J ∇ V + 12 ∇ V T J T T T DT J ∇ V + V ( q k +1 )= 12 q T k Dq k + q T k S T AT ∇ V + 12 ∇ V T J T T T AT ∇ V + V ( q k +1 ) , Applying the properties above yields the following, q T k S T AT ∇ V = q T k T T A ∇ V = q T k ( I − S T ) ∇ V,J T T T AT = J T T T ( S − I ) = J T T + JT T , ∇ V T J T T T AT ∇ V = 12 ∇ V T ( J T T + JT T ) ∇ V = ( ∇ V ) T T T J ∇ V. Substituting these into the expression for H ( q k +1 ) yields(16) H ( q k +1 ) = 12 q T k Dq k + [( I − S ) q k ] T ∇ V − ( ∇ V ) T ( T J ) ∇ V + V ( q k +1 )= 12 q T k Dq k + ( q k − Sq k − T J ∇ V ) T ∇ V + V ( q k +1 )= 12 q T k Dq k + ( q k − q k +1 ) T ∇ V ( q k , q k +1 ) + V ( q k +1 )= 12 q T k Dq k − V ( q k +1 ) + V ( q k ) + V ( q k +1 )= H ( q k ) , which proves that the Hamiltonian is preserved. (cid:3) One significant advantage of exponential integrators is they allow the numericalmethod to be implemented using fixed point iterations as opposed to the more
XUEFENG SHEN AND MELVIN LEOK computationally expensive Newton iterations. Recall that classical implicit Runge-Kutta methods Y i = y n + h · s (cid:88) j =1 a ij f ( t n + c j h, Y j ) ,y n +1 = y n + h · s (cid:88) i =1 b i f ( t n + c i h, Y i ) , have a form that naturally lends itself to fixed point iterations. However, when ∂f∂y has a large spectral radius, the timestep is forced to be very small in orderto guarantee that the fixed point iteration converges. The alternative is to usea Newton type iteration, which is time consuming since we need to perform LUdecomposition ( O ( n ) complexity) during each iteration. This is the problem weface for the stiff semilinear system (1) when the coefficient matrix A has a largespectral radius. In contrast, in both the exponential midpoint rule q k +1 = e Ah q k + h · e A h f (cid:16) e A h q k + e − A h q k +1 (cid:17) , and the energy-preserving exponential integrator q k +1 = e Ah q k + (cid:90) h e Aτ dτ · J ∇ V ( q k , q k +1 ) , the matrix A only appears in the exponential term e Ah . Since A is skew-symmetric,this term is an orthogonal matrix, which has spectral radius 1. Thus, fixed pointiterations can be used to implement the exponential integrator, regardless of thestiffness of A . 4. Numerical Methods
We consider two Hamilton PDEs here, namely the nonlinear Schr¨odinger equa-tion,(17) iψ t + ψ xx − | ψ | ψ = 0 , in which ψ = u + iv is the wave function with real part u and imaginary part v ,and has the following equivalent form,(18) (cid:40) u t = − v xx + 2( u + v ) v,v t = u xx − u + v ) u ;as well as the KdV equation,(19) u t + uu x + u xxx = 0 . To discretize the two PDEs, we impose 2 π periodic boundary conditions. Given asmooth 2 π periodic function f ( x ), on the interval [0 , π ], choose 2 n + 1 equispacedinterpolation points x j = jh, j = 0 , , , . . . , n , h = π n +1 . Given nodal values { v j } nj =0 , there exists a unique trigonometric polynomial v ( x ) with degree less orequal n , such that, v ( x j ) = v j (see, for example, [1]). v ( x ) = n (cid:88) k = − n ˆ v k e ikx , ˆ v k = 12 n + 1 n (cid:88) j =0 v j e − ikx j . EOMETRIC EXPONENTIAL INTEGRATORS 9
By substituting the expression for the coefficients ˆ v k , we obtain(20) v ( x ) = n (cid:88) k = − n (cid:32) n + 1 n (cid:88) j =0 v j e − ikx j (cid:33) e ikx = n (cid:88) j =0 v j (cid:32) n (cid:88) k = − n n + 1 e ik ( x − x j ) (cid:33) = n (cid:88) j =0 v j φ ( x − x j )= n (cid:88) j =0 v j φ j ( x ) , where φ ( x ) = n (cid:88) k = − n n + 1 e ikx = 12 n + 1 sin(( n + ) x )sin( x ) . From this, we see that { e ikx } nk = − n and { φ j } nj =0 are equivalent orthogonal basesfor the trigonometric polynomial function space, and each such function can beparametrized by either the nodal values { v j } nj =0 or Fourier coefficients { ˆ v k } nk = − n .They represent the same function, but with respect to two different bases. Thetransformation between { v j } nj =0 and { ˆ v k } nk = − n can be performed using the FastFourier transformation (FFT), which has O ( n log n ) complexity.The first and second-order differentiation matrices [10] with respect to the rep-resentation in terms of nodal values { v j } nj =0 are given by( D ) kj = , k = j, ( − ( k − j ) ( k − j ) h ) , k (cid:54) = j, ( D ) kj = − n ( n +1)3 , k = j, ( − ( k − j +1) cos( ( k − j ) h )2 sin ( ( k − j ) h ) , k (cid:54) = j, respectively. However, with respect to the representation in terms of Fourier coef-ficients { ˆ v k } nk = − n , they are diagonal,ˆ D = diag( ik ) nk = − n , ˆ D = diag( − k ) nk = − n . Later, we will see that this observation is critical to a fast implementation of theproduct of matrix functions with vectors. We can also define a third-order differen-tiation matrix D , which has the property D = D D = D D , and it is diagonalwith respect to the Fourier coefficients ˆ D = diag( − ik ) nk = − n .4.1. Nonlinear Schr¨odinger equation.
We perform a semi-discretization of (18)by discretizing the solution u , v in space using their corresponding nodal values { q j } and { p j } . Applying the pseudospectral method, we obtain the following system ofODEs,(21) (cid:40) ˙ q = − D p + 2( q + p ) p, ˙ p = D q − q + p ) q, where the nonlinear term ( q + p ) p is computed elementwise, and represents thevector consisting of { ( q j + p j ) p j } entries. We adopt this notation throughout therest of the paper for brevity. Then, (21) can be expressed as,(22) ddt (cid:18) qp (cid:19) = (cid:18) − D D (cid:19) (cid:18) qp (cid:19) + (cid:18) q + p ) p − q + p ) q (cid:19) , where A = (cid:18) − D D (cid:19) = (cid:18) I − I (cid:19) (cid:18) − D − D (cid:19) = J · D,f ( q, p ) = (cid:18) q + p ) p − q + p ) q (cid:19) = (cid:18) I − I (cid:19) (cid:18) q + p ) q q + p ) p (cid:19) = J · ∇ V ( q, p ) , and V ( q, p ) = ( q + p ) . It is easy to verify that J is skew-symmetric, D issymmetric, and JD = DJ . Thus, (22) is a semilinear Poisson system.To apply the exponential midpoint rule, we need to compute the product of amatrix function and a vector, which has the form e Ah (cid:18) qp (cid:19) , e (cid:16) − D D (cid:17) = ∞ (cid:88) k =0 k ! (cid:0) − D D (cid:1) k , = ∞ (cid:88) k =0 k )! (cid:16) ( − k · ( D ) k
00 ( − k · ( D ) k (cid:17) + ∞ (cid:88) k =0 k + 1)! (cid:16) − k +1 · ( D ) k +1 ( − k · ( D ) k +1 (cid:17) = (cid:16) cos( D ) − sin( D )sin( D ) cos( D ) (cid:17) . Theorem 4.
Suppose that D , D are the second and third-order differentiationmatrices, respectively, q = { q j } nj =0 is a vector with Fourier transform F [ q ] = ˆ q = { ˆ q k } nk = − n , and f is an analytic function, then f ( D ) q = F − [diag( f ( − k ))ˆ q ] , f ( D ) q = F − [diag( f ( − ik ))ˆ q ] , where F − is the inverse Fourier transform.Proof. Recall that matrix D is diagonalizable with eigenvalues λ k = − k , andcorresponding eigenvectors e k = { e ikx j } nj =0 , f ( D ) q = f ( D ) (cid:32) n (cid:88) k = − n ˆ q k · e k (cid:33) = n (cid:88) k = − n ˆ q k · f ( D ) e k = n (cid:88) k = − n ˆ q k · f ( λ k ) e k = F − [diag( f ( λ k ))ˆ q ]= F − [diag( f ( − k ))ˆ q ] . EOMETRIC EXPONENTIAL INTEGRATORS 11
Notice D is also diagonalizable with eigenvalues λ k = − ik , and correspondingeigenvectors e k , so the property that f ( D ) q = F − [diag( f ( − ik ))ˆ q ] can be verifiedin the same way. (cid:3) By Theorem 4, we have that(23) e (cid:16) − D D (cid:17) h (cid:18) qp (cid:19) = (cid:18) cos( D h ) − sin( D h )sin( D h ) cos( D h ) (cid:19) (cid:18) qp (cid:19) = (cid:18) cos( D h ) · q − sin( D h ) · p sin( D h ) · q + cos( D h ) · p (cid:19) = (cid:18) F − [cos( k h )ˆ q k + sin( k h )ˆ p k ] F − [cos( k h )ˆ p k − sin( k h )ˆ q k ] (cid:19) . In summary, the exponential midpoint rule for the nonlinear Schr¨odinger equa-tion is given by z k +1 = e Ah z k + h · e A h f (cid:18) e A h z k + e − A h z k +1 (cid:19) , where z k = (cid:18) q k p k (cid:19) , z k +1 = (cid:18) q k +1 p k +1 (cid:19) , A = (cid:18) − D D (cid:19) , and e Ah z k can be efficientlycalculated using (23).For the energy-preserving exponential integrator for the nonlinear Schr¨odingerequation, z k +1 = e Ah z k + (cid:90) h e Aτ dτ · J ∇ V ( z k , z k +1 ) . Here, (cid:90) h e (cid:16) − D D (cid:17) τ dτ = (cid:90) h (cid:16) cos( D τ ) − sin( D τ )sin( D τ ) cos( D τ ) (cid:17) dτ = h · (cid:32) sin( D h ) D h cos( D h ) − D h − cos( D h ) D h sin( D h ) D h (cid:33) , and we can construct the discrete gradient ∇ V ( z k , z k +1 ) = (cid:18) q k + + p k + ) · q k + q k + + p k + ) · p k + (cid:19) , where q k + = q k + q k +1 , p k + = p k + p k +1 , (24) q k + = q k + q k +1 , p k + = p k + p k +1 . (25)Notice that ∇ V ( z k , z k +1 ) is symmetric with respect to z k and z k +1 , and it canbe verified that it satisfies (11). A classical discrete gradient method can be con-structed as follows, z k +1 = z k + hJ ∇ H ( z k , z k +1 ) , (26)where ∇ H ( z k , z k +1 ) = (cid:18) − D − D (cid:19) z k + z k +1 ∇ V ( z k , z k +1 ) . The method described by (26) is very similar to classical midpoint rule, the onlydifference is in ∇ V ( z k , z k +1 ), q k + is used, while ( q k + ) is used in the midpointrule, so (26) can be viewed as a modified midpoint rule.4.2. KdV equation.
Rewrite (19) as u t = (cid:16) − ∂∂x (cid:17)(cid:16) u + u xx (cid:17) , then apply pseudospectral semi-discretization to obtain the following system,˙ q = ( − D ) (cid:16) q + D q (cid:17) = ( − D ) D q + ( − D ) (cid:16) q (cid:17) , which has the form of a semilinear Poisson system (5),˙ q = J ( Dq + ∇ V ( q )) = J ∇ H ( q ) , where J = − D , D = D , A = JD = − D , ∇ V ( q ) = q , H ( q ) = q T D q + q .The exponential midpoint rule for KdV reads as follows, q k +1 = e − D h q k + h · e − D h f (cid:16) e − D h q k + e D h q k +1 (cid:17) , and the energy preserving exponential integrator is given by, q k +1 = e − D h q k + (cid:90) h e − D τ dτ · ( − D ) ∇ V ( q k , q k +1 ) , with discrete gradient ∇ V ( q k , q k +1 ) = ( q k + q k · q k +1 + q k +1 ). A related classicaldiscrete gradient method can be constructed as follows, ∇ H ( q k , q k +1 ) = ( D ) q k + q k +1 ∇ V ( q k , q k +1 ) . (27)By Theorem 4, in each iteration, the matrix function and vector product can beimplemented as e − D h q = F − [ e ik h ˆ q k ] , and (cid:32) (cid:90) h e − D τ dτ (cid:33) q = (cid:16) e − D h − I − D (cid:17) q = F − (cid:104) e ik h − ik ˆ q k (cid:105) . Numerical Experiments
Nonlinear Schr¨odinger equation.
In Table 3 below, n denotes the numberof nodes we discretize the spatial domain with, and the data in the table indi-cates the maximum timesteps for which the nonlinear solver converges. The firsttwo columns correspond to midpoint rule, using fixed point iteration and Newtontype iteration; the third column corresponds to the exponential midpoint rule, thefourth column is discrete gradient method (26), and the last column is the energy-preserving exponential integrator, all implemented using fixed point iterations.We observe that when the midpoint rule is implemented using fixed point itera-tions, there is a significant decay in the allowable timestep as the spatial resolutionis increased, whereas the Newton type iteration allows a relatively large timestepthat is independent of the spatial resolution. In contrast, the midpoint exponentialmethod exhibits a slower rate of decrease in allowable timestep when using fixed EOMETRIC EXPONENTIAL INTEGRATORS 13 midpoint midpoint exp discrete gradient energy exp n fixed point Newton fixed point fixed point fixed point11 0.02 0.1 0.1 0.02 0.121 0.01 0.1 0.08 0.01 0.141 4 × − × − × − × − − − × − × − × − × − × − × − × − × − × − × − Table 3.
Maximum timestep for nonlinear solver to converge asa function of numerical integrator, nonlinear solver, and spatialresolution.point iterations. When using fixed point iterations, the discrete gradient method,which is an energy preserving method, the allowable timestep behaves similarly tothe midpoint rule, and in contrast, the energy preserving exponential integratorhas an allowable timestep that is independent of the spatial resolution. log(n) -11-10-9-8-7-6-5-4-3-2 l og ( t i m e s t ep ) Exponential midpointmidpoint
Figure 1.
Maximum timestep for which fixed point iterationsconverge as a function of the spatial resolution.
In Figure 1, we observe that the allowable timestep when using fixed pointiteration scale like n − for the classical midpoint rule, and n − for the midpointexponential rule. time -0.500.511.522.533.5 ene r g y e rr o r -5 energy error for midpoint ex (a) Energy error time t r a j e c t o r y e rr o r -3 trajectory error for midpoint ex (b) Trajectory error
Figure 2.
Error plots for the exponential midpoint rule, n = 161, timestep = 0 . time -8-7-6-5-4-3-2-10 ene r g y e rr o r -11 energy error for energy preserving ex (a) Energy error time t r a j e c t o r y e rr o r trajectory error for energy preserving ex (b) Trajectory error
Figure 3.
Error plots for the energy preserving exponential inte-grator, n = 161, timestep = 0 . .
1, we see in Figure 3 that the energy is stillpreserved approximately to within machine error.
EOMETRIC EXPONENTIAL INTEGRATORS 15
We now consider the diagonally implicit symplectic exponential (DISEX) inte-grator with six stages: b = 0 . b = 1 . b = 2 . b = 0 . b = − . b = − . . The energy and trajectory error is shown in Figure 4. Observe that the en-ergy error is small and bounded, as expected of a symplectic integrator, and thetrajectory error is small as well, as expected of a higher-order numerical integrator. time -10-505 ene r g y e rr o r -8 energy error for DISEX (a) Energy error time t r a j e c t o r y e rr o r -6 trajectory error for DISEX (b) Trajectory error
Figure 4.
Error plots for the 6 stage diagonally implicit symplec-tic exponential integrator (DISEX), n = 161, timestep = 0 . KdV.
We simulate the KdV equation, u t + uu x + νu xxx = 0 , where ν = 5 × − . As before, we explore how the maximum timestep for whichthe fixed point iteration converges depends on the choice of numerical integrator,and the spatial resolution of the semi-discretization. In Table 4, n is the numberof nodes we use to discretize the spatial domain, the data in the table indicatesthe maximum timestep for which fixed point iterations converge. For the midpointrule, the timestep decreases like n − for fixed point iterations, and a comparabletimestep is required for Newton iterations which is thereby too costly to implement.The classical discrete gradient method for the KdV equation is given by (27), andit exhibits the same timestep restrictions as the midpoint rule. In contrast, boththe exponential midpoint and energy preserving exponential integrator allow ratherlarge timesteps that are independent of the spatial resolution.In Figure 5, we observe that the exponential midpoint rule has an energy errorthat is small and bounded, as is typical for a symplectic integrator, and the trajec-tory error grows linearly. In Figure 6, the energy preserving exponential integratorhas an energy that is preserved to within machine precision, and the trajectoryerror grows linearly. n midpoint midpoint exp discrete gradient energy exp401 4 × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − Table 4.
Maximum timestep for fixed point iteration to convergeas a function of numerical integrator, and spatial resolution. time ene r g y e rr o r -5 energy error for midpoint ex (a) Energy error time t r a j e c t o r y e rr o r trajectory error for midpoint ex (b) Trajectory error
Figure 5.
Error plots for the exponential midpoint rule, n = 401, timestep = 5 × − . time -101234567 ene r g y e rr o r -12 energy error for energy preserving ex (a) Energy error time t r a j e c t o r y e rr o r trajectory error for energy preserving ex (b) Trajectory error
Figure 6.
Error plots for the energy preserving exponential inte-grator, n = 401, timestep = 0 . EOMETRIC EXPONENTIAL INTEGRATORS 17 Summary
For semilinear Poisson system, we have developed two types of geometric expo-nential integrators, one that preserves the Poisson structure, and the other preservesthe energy. They allow fixed point iteration methods to be used for significantlylarger timesteps as compared to non-exponential integrators, such as the classicalmidpoint rule and the standard discrete gradient method, which require the use ofNewton type iterations to converge with comparably large timesteps. This resultsin substantial computational savings, particularly when applied to the simulationof semi-discretized Hamiltonian PDEs. They also exhibit long time stability andconservation of the first integrals. We notice that in practice, energy preservingexponential integrators are more stable and allow for larger timesteps than sym-plectic exponential integrators. In future work, we will develop higher-order energypreserving exponential integrators.
Acknowledgements
This research has been supported in part by NSF under grants DMS-1010687,CMMI-1029445, DMS-1065972, CMMI-1334759, DMS-1411792, DMS-1345013.
References [1] K. Atkinson and W. Han.
Theoretical numerical analysis , volume 39 of
Textsin Applied Mathematics . Springer, Dordrecht, third edition, 2009. A functionalanalysis framework.[2] E. Hairer, C. Lubich, and G. Wanner.
Geometric numerical integration , vol-ume 31 of
Springer Series in Computational Mathematics . Springer-Verlag,Berlin, second edition, 2006. Structure-preserving algorithms for ordinary dif-ferential equations.[3] M. Hochbruck and A. Ostermann. Exponential Runge-Kutta methods forparabolic problems.
Appl. Numer. Math. , 53(2-4):323–339, 2005.[4] M. Hochbruck and A. Ostermann. Exponential integrators.
Acta Numer. , 19:209–286, 2010.[5] B. Leimkuhler and S. Reich.
Simulating Hamiltonian dynamics , volume 14of
Cambridge Monographs on Applied and Computational Mathematics . Cam-bridge University Press, Cambridge, 2004.[6] J. E. Marsden and M. West. Discrete mechanics and variational integrators.
Acta Numer. , 10:357–514, 2001.[7] B. Minchev.
Exponential integrators for semilinear problems . PhD thesis,University of Bergen, 2004.[8] G. R. W. Quispel and G. S. Turner. Discrete gradient methods for solvingODEs numerically while preserving a first integral.
J. Phys. A , 29(13):L341–L349, 1996.[9] J. C. Simo and O. Gonzalez. Assessment of energy-momentum and symplecticschemes for stiff dynamical systems.
Proc. ASME Winter Annual Meeting,New Orleans, Louisiana. , 1993.[10] L. N. Trefethen.
Spectral methods in MATLAB , volume 10 of
Software, Envi-ronments, and Tools . Society for Industrial and Applied Mathematics (SIAM),Philadelphia, PA, 2000.[11] X. Wu, K. Liu, and W. Shi. An extended discrete gradient formula for multi-frequency oscillatory Hamiltonian systems. In
Structure-Preserving Algorithmsfor Oscillatory Differential Equations II , pages 95–115. Springer, 2015.[12] W. J. Zhu and M. Z. Qin. Poisson schemes for Hamiltonian systems on Poissonmanifolds.