Compositions of pseudo-symmetric integrators with complex coefficients for the numerical integration of differential equations

In this paper, we are concerned with the construction and analysis of a new class of methods obtained as double jump compositions with complex coefficients and projection on the real axis. It is shown in particular that the new integrators are symmetric and symplectic up to high orders if one uses a symmetric and symplectic basic method. In terms of efficiency, the aforementioned technique requires fewer stages than standard compositions of the same orders and is thus expected to lead to faster methods.

Authors

• 7 publications
• 6 publications
• 4 publications
• 61 publications
• High order integrators obtained by linear combinations of symmetric-conjugate compositions

A new family of methods involving complex coefficients for the numerical...
06/11/2021 ∙ by Fernando Casas, et al. ∙ 0

• On symmetric-conjugate composition methods in the numerical integration of differential equations

We analyze composition methods with complex coefficients exhibiting the ...
01/11/2021 ∙ by Sergio Blanes, et al. ∙ 0

• Diagonal asymptotics for symmetric rational functions via ACSV

We consider asymptotics of power series coefficients of rational functio...
04/29/2018 ∙ by Yuliy Baryshnikov, et al. ∙ 0

• Evolutionary Design of Numerical Methods: Generating Finite Difference and Integration Schemes by Differential Evolution

Classical and new numerical schemes are generated using evolutionary com...
12/30/2013 ∙ by C. D. Erdbrink, et al. ∙ 0

• Applying splitting methods with complex coefficients to the numerical integration of unitary problems

We explore the applicability of splitting methods involving complex coef...
04/06/2021 ∙ by S. Blanes, et al. ∙ 0

• Discovering New Runge-Kutta Methods Using Unstructured Numerical Search

Runge-Kutta methods are a popular class of numerical methods for solving...
10/30/2019 ∙ by David K. Zhang, et al. ∙ 0

• Weak SINDy: A Data-Driven Galerkin Method for System Identification

We present a weak formulation and discretization of the system discovery...
05/09/2020 ∙ by Daniel A. Messenger, et al. ∙ 0

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Given a differential equation

 ˙x≡dxdt=f(x),x(0)=x0, (1)

composition methods constitute a powerful technique to raise the order of a given integrator applied to (1) with time-step , as high as might be required, by considering expressions of the form

 ϕτ=ψγ1τ∘ψγ2τ∘⋯∘ψγsτ, (2)

where the coefficients are appropriately chosen so as to satisfy some universal algebraic conditions [HLW06, MSS99, CM09]. It is known in particular that if is of order , i.e. satisfies

 φτ(x0)−ψτ(x0)=O(τk+1),

where denotes the exact flow of (1), then will be at least of order (i.e., local error ) if the following two conditions are satisfied

 (i)s∑i=1γi=1 and (ii)s∑i=1γk+1i=0. (3)

Given that these two equations have no real solution for odd

and arbitrary , a series of authors (e.g. [Suz90, Yos90]) suggested to start from a second-order method and to consider symmetric compositions only, i.e., schemes with coefficients satisfying the additional condition

 γs+1−i=γi,i=1,…,s.

This has led to so-called triple-jump compositions (, ) obtained by iterating the process described above to construct a sequence of symmetric methods with even orders (see, e.g., [HLW06] pp. 44).

In spite of its simplicity, the triple-jump rationale leads to inefficiencies for high orders as compared to methods obtained by solving directly the order conditions [HLW06]

. On top of this, it also suffers from the occurrence of negative time-steps, although this fact is not specific to triple-jump methods and concerns all composition or splitting methods of orders higher than two. This, of course, is a severe limiting factor for equations where the vector field (usually an operator) is not reversible, the prototypical example of which being the heat equation. To circumvent this difficulty, several authors have suggested to use complex time-steps (or complex coefficients) in the context of parabolic equations

[CCDV09, HO09]. One indeed easily sees that, already for , solutions of equations exist in .

Generally speaking, suppose that is an integrator of order , denoted in the sequel for clarity, and consider the composition (2) with ,

 S[k+1]τ=S[k]γ1τ∘S[k]γ2τ. (4)

Then, if the coefficients verify conditions , that is to say if

 γ1=¯γ2≡γ=12+i2sin(2ℓ+1k+1π)1+cos(2ℓ+1k+1π)for⎧⎨⎩−k2≤ℓ≤k2−1if k is even−k+12≤ℓ≤k−12if k is odd, (5)

then (4) results in a method of order , which can subsequently be used to generate recursively higher order composition schemes by applying the same procedure. The choice ,

 γ=γ[k]:=12+i2sin(πk+1)1+cos(πk+1)=12+i2tan(π2(k+1))=12cos(π2(k+1))eπ2(k+1) (6)

gives the solutions with the smallest phase. If we start with a symmetric method of order 2, , and apply composition (4) with corresponding coefficients (6), we can construct the following sequence of methods:

 S[2]τ⟶S[3]τ⟶S[4]τ⟶S[5]τ⟶S[6]τ,

all of which have coefficients with positive real part [HO09]. The final method of order involves evaluations of the basic scheme . By contrast, there are composition methods of order (both with real and complex coefficients) involving just evaluations of [BCCM13, Yos90]. It is thus apparent that this direct approach does not lead to cost-efficient high-order schemes.

One should remark that the composition (4) does not provide a time-symmetric method, i.e., is not the identity map, even if happens to be symmetric. As we have mentioned before, symmetry allows to raise the order by two at each iteration by considering the triple-jump composition

 S[2k+2]τ=S[2k]γ1τ∘S[2k]γ2τ∘S[2k]γ1τ (7)

starting from a symmetric method. Apart from the real solution, the complex one with the smallest phase is

 γ1=eiπ/(k+1)21/(k+1)−2eiπ/(k+1),γ2=1−2γ1, (8)

and symmetric methods up to order with coefficients having positive real part are possible if one starts with a symmetric second-order scheme555It is actually possible to reach order 14 if, in the construction, one uses formula (7) alternatively with coefficients and coefficients [CCDV09].. These order barriers has been rigorously proved in [BCCM13].

The simple third-order scheme (4) corresponding to has been in fact rediscovered several times in the literature [BS91, CCDV09, Cha03, HO09, Suz92]. In particular, it was shown in [Cha03] that the method, when applied to the two-body Kepler problem, behaves indeed as a fourth-order integrator, the reason being attributed to the fact that the main error term in the asymptotic expansion is purely imaginary. In this note we elaborate further the analysis and provide a comprehensive study of the general composition (4), paying special attention to the qualitative properties the method shares with the continuous system (1). In addition, we show how it is possible to combine compositions and a trivial linear combination to raise the order, while still preserving the qualitative properties of the basic integrator up to an order higher than of the method itself.

2 Composition and pseudo-symmetry or pseudo-symplecticity

In what follows, we will assume for convenience that all values of in (1) lie in a compact set where the function is smooth. Before starting the analysis, it is worth recalling the notions of adjoint method and symplectic flow.

The adjoint method of a given method is the inverse map of the original integrator with reversed time step , i.e., . A symmetric method satisfies [Cha15, HLW06].

The vector field in (1) is Hamiltonian if there exists a function such that , where and is the basic canonical matrix. Then, the exact flow of (1) is a symplectic transformation, for [BC16, SSC94].

It then makes sense introducing the following definitions, taken from [CL98] and [AC98]:

Definition 1

Let be a smooth and consistent integrator:

1. it is pseudo-symmetric of pseudo-symmetry order if for all sufficiently small , the following relation holds true:

 ψ∗τ=ψτ+O(τq+1), (9)

where the constant in the -term depends on bounds of derivatives of on .

2. it is pseudo-symplectic of pseudo-symplecticity order if for all sufficiently small , the following relation holds true when is applied to a Hamiltonian system:

 (ψ′τ)TJψ′τ=J+O(τr+1), (10)

where the constant in the -term depends on bounds of derivatives of on .

Remark 1

A symmetric method is pseudo-symmetric of any order , whereas a method of order is pseudo-symmetric of order . A similar statement holds for symplectic methods.

As a first illustration of Definition 1, let us consider again a symmetric 2nd-order method and form the composition

 ψ[3]τ=S[2]γτ∘S[2]¯γτ

with . Then, if the vector field under consideration is real-valued, its real part

 R(ψ[3]τ)=12(ψ[3]τ+¯¯¯¯ψ[3]τ)=12(S[2]γτ∘S[2]¯γτ+S[2]¯γτ∘S[2]γτ).

is a method of order 4 and pseudo-symmetric of pseudo-symmetry order 7. This result is a consequence of the fact that

 (ψ[3]τ)∗=¯¯¯¯ψ[3]τ

and the following general statement, which lies at the core of the construction procedure described in this paper.

Proposition 1

Let be any consistent smooth method for equation (1) and consider the new method

 Rτ=12(ψτ+ψ∗τ).

Assume also that is pseudo-symmetric of order . Then is of pseudo-symmetry order . If is furthermore of pseudo-symplecticity order , then is of pseudo-symplecticity order .

Proof:

By assumption, there exists a smooth function , defined for all in a compact set and for all sufficiently small real , such that

 ψ∗τ=ψτ+τq+1δτ or% ψ−1−τ=ψτ+τq+1δτ or ψ−1τ=ψ−τ+(−τ)q+1δ−τ, (11)

so that

 Rτ=ψτ+12τq+1δτ.

Composing the third relation of (11) from the left by , we obtain

 id=ψτ∘ψ−τ+(−τ)q+1ψ′τ∘ψ−τ⋅δ−τ+O(τ2(q+1)), (12)

where the -term depends on bounds of the derivatives of and on . Similarly, composing the second relation of (11) from the right by , we get

 id=ψτ∘ψ−τ+τq+1δτ∘ψ−τ. (13)

As a consequence, we have

 τq+1δτ∘ψ−τ=(−τ)q+1ψ′τ∘ψ−τ⋅δ−τ+O(τ2(q+1)).

We are then in position to write

 Rτ∘R−τ = (ψτ+12τq+1δτ)∘(ψ−τ+12(−τ)q+1δ−τ) = ψτ∘ψ−τ+12(−τ)q+1ψ′τ∘ψ−τ⋅δ−τ+12τq+1δτ∘ψ−τ+O(τ2(q+1)) = id+O(τ2(q+1)),

which proves the first statement. Now, if is in addition of pseudo-symplecticity order , then its adjoint is also of pseudo-symplecticity order , so that relation (11) leads to

 J+O(τr+1) =(∂xψ∗τ)TJ∂xψ∗τ=(ψ′τ+τq+1δ′τ)TJ(ψ′τ+τq+1δ′τ) =J+O(τr+1)+τq+1((δ′τ)TJψ′τ+(ψ′τ)TJδ′τ)+O(τ2(q+1)),

which implies that

 τq+1((δ′τ)TJψ′τ+(ψ′τ)TJδ′τ)=O(τmin(2(q+1),r+1)).

As an immediate consequence, we have that

 (R′τ)TJR′τ=J+O(τmin(2(q+1),r+1))

which proves the second statement.

This result can be rendered more specific as follows:

Proposition 2

Let be a smooth method of order and pseudo-symmetry order . Let us consider the composition method

 ψ[2n+1]τ=S[2n]γ1τ∘S[2n]γ2τ, (14)

where the coefficients and satisfy both relations and . Then the method

 ^Rτ=12(ψ[2n+1]τ+¯¯¯¯ψ[2n+1]τ) (15)

is of order

 {2n+1 if q=2n+1,2n+2 if q≥2n+2 (16)

when the vector field in (1) is real, and of pseudo-symmetry order

 {2n+1 if q=2n+1,min(q,4n+3) if q≥2n+2. (17)

If in addition, is a (real) Hamiltonian vector field and is of pseudo-symplecticity order , then is of pseudo-symplecticity order

 {min(r,2n+1) if q=2n+1,min(q,r,4n+3) if q≥2n+2. (18)
Remark 2

Note that in Proposition 2, one has necessarily . This can be seen straightforwardly by a direct computation of with .

Proof:

Noticing that and are complex conjugate and (1) is real, and taking into account that is of pseudo-symmetry order , we have

 ¯¯¯¯ψ[2n+1]τ =S[2n]γ2τ∘S[2n]γ1τ=((S[2n]γ2τ)∗+O(τq+1))∘((S[2n]γ1τ)∗+O(τq+1)) =(S[2n]γ2τ)∗∘(S[2n]γ1τ)∗+O(τq+1)=(ψ[2n+1]τ)∗+O(τq+1). (19)

Moreover, by construction, is at least of order , so that

 ψ[2n+1]τ+O(τ2n+2)=¯¯¯¯ψ[2n+1]τ=(ψ[2n+1]τ)∗+O(τ2n+2), (20)

and altogether

 ^Rτ =Rτ+O(τmax(2n+2,q+1)). (21)

Now, since the pseudo-symmetry order of is at least , the method

 Rτ=12(ψ[2n+1]τ+(ψ[2n+1]τ)∗)

is, according to Proposition 1, of pseudo-symmetry order and of pseudo-symplecticity order . The first (16), second (17) and third (18) statements on orders then follow from (21).

In the Appendix we provide an alternative proof of Proposition 2 based on the Lie formalism, which allow us, in addition, to generalise the previous result on pseudo-symplecticity to other geometric properties the continuous system may possess (such as in volume preserving flows, isospectral flows, differential equations evolving on Lie groups, etc.).

Notice that, according with Proposition 2, if we start from , that is to say from a basic symmetric () and symplectic () method of order , we get a method of order that is pseudo-symmetric and pseudo-symplectic of order just by considering the simple composition (14) and taking the real part of the output at each time step. If this technique is applied to a symmetric and symplectic method of order , i.e. with , then is of order and pseudo-symmetric and pseudo-symplectic of order .

Let us consider, in particular, the -order symmetric scheme (7) with as basic scheme. Then, the resulting -order integrator only requires the evaluation of second-order methods , whereas the corresponding -order scheme obtained by the triple-jump technique involves evaluations. This number is reduced to by considering general compositions of [BCCM13]. If we take this -order composition of schemes as the basic method , the resulting integrator of order , , involves the evaluation of , whereas evaluations are required by pure composition methods. Notice that is pseudo-symmetric and pseudo-symplectic of order 15, so that for values of sufficiently small, it preserves effectively the symmetry up to round-off error while the drift in energy for Hamiltonian systems is hardly noticeable.

3 Families of pseudo-symplectic methods

There is another possibility to increase the order, though, and it consists in applying the technique of Proposition 2 recursively. Thus, if denote by the method of eq. (15), we propose to apply the following recurrence:

 For i=2,3,… (22) Φ(i)τ=^R(i−1)γ[2i]τ∘^R(i−1)¯γ[2i]τ ^R(i)τ=12(Φ(i)τ+¯¯¯¯Φ(i)τ)

where is given by (6). Then, according with Proposition 2, it is possible to raise the order up to the pseudo-symmetry order of the underlying basic method . Thus, in particular, the maximum order one can achieve by applying this technique to the basic symmetric method is 7, whereas if we start with a basic symmetric method of order 4, , the maximum order is . It is from a symmetric method of order and so on and so forth.

To give an assessment of the computational cost of the methods obtained by applying this type of composition, we notice that the computation of and required to form by (22) at the intermediate stages can be done in parallel, whereas at the final stage it only requires taking the real part. Thus, the method of order 6 constructed recursively from only requires the effective computation of 4 basic methods .

Starting from a symmetric second-order method , say Strang splitting for instance, it is important to monitor the sign of the real part of all coefficients involved in the previous iteration. It is immediate to see that in the recursive construction

 S[2]τ→^R(1)τ→^R(2)τ→^R(3)τ

envisaged in the recurrence (22), the basic method is used with the following coefficients

 i=1: γ[2],¯γ[2] i=2: γ[4]γ[2],¯γ[4]γ[2],γ[4]¯γ[2],¯γ[4]¯γ[2] i=3: γ[6]γ[4]γ[2],γ[6]¯γ[4]γ[2],γ[6]γ[4]¯γ[2],γ[6]¯γ[4]¯γ[2],¯γ[6]γ[4]γ[2],¯γ[6]¯γ[4]γ[2],¯γ[6]γ[4]¯γ[2],¯γ[6]¯γ[4]¯γ[2]

Given the expression of (see (6)), these coefficients have arguments of the form

 π2i∑j=1±12j+1=π2(±13±15±⋯±12i+1),i=1,2,3,

so that their maximum argument is

 π23∑j=112j+1.

For all the coefficients to have positive real parts, a necessary and sufficient condition is thus that

 3∑j=112j+1≤1.

It clearly holds for methods , and , of respective orders , and , since . Similarly, starting form a symmetric method of order having real or complex coefficients with maximum argument , the condition becomes

 2θ4π+4∑j=112j+3=2θ4π+18883465≤1.

For instance, suppose that in (1) can be split as , so that the exact -flows and corresponding to and , respectively, can be computed exactly. Then, the following composition

 S[4]τ=φ[b]b1τ∘φ[a]a1τ∘φ[b]b2τ∘φ[a]a2τ∘φ[b]b3τ∘φ[a]a2τ∘φ[b]b2τ∘φ[a]a1τ∘φ[b]b1τ (23)

with

 b1=110−130i,b2=415+215i,b3=415−15i and a1=a2=a3=a4=14

provides a 4th-order symmetric scheme (see [CCDV09]). Taking (23) as basic method we get so that

 2θ4π+18883465<0.409666+18883465<0.96<1

and thus all methods

 S[4]τ→^R(1)τ→^R(2)τ→^R(3)τ→^R(4)τ

of respective orders , , , and obtained by the procedure (22) have all their coefficients with positive real parts. As far as the part is concerned, the maximum argument is less than .

4 Numerical experiments

In this section we illustrate the previous results on several numerical examples, comprising Hamiltonian systems and partial differential equations of evolution previously discretised in space.

4.1 Harmonic oscillator

We first consider the simple harmonic oscillator, with Hamiltonian

 H=T(p)+V(q)=12p2+12q2.

If we denote by the exact matrix evolution associated with the Hamiltonians , and , i.e., , then

 MH(τ)=(cos(τ)sin(τ)−sin(τ)cos(τ)),MT(τ)=([]ll1τ01),MV(τ)=(10−τ1),

respectively. We take as basic symmetric (and symplectic) scheme the leapfrog/Strang splitting:

 S[2]τ=MT(τ/2)MV(τ)MT(τ/2) (24)

and compute the first three iterations in (22). In Table 1 we collect the main term in the truncation error for the resulting integrators , . We also check their time-symmetry and the preservation of the symplectic character of the approximate solution matrix by computing its determinant (a matrix is symplectic iff

). One can observe that these results are in agreement with the previous estimates.

Next we take initial conditions , integrate until the final time with , , and and compute the relative error in energy along the evolution. The result is depicted in Figure 1. We see that for and the error in energy is almost constant for a certain period of time, and then there is a secular growth proportional to .

4.2 Kepler problem

Next, we consider the two-dimensional Kepler problem with Hamiltonian

 H(q,p)=T(p)+V(q)=12pTp−μ1r. (25)

Here , , is the gravitational constant and is the sum of the masses of the two bodies. Taking and initial conditions

 q1(0)=1−e,q2(0)=0,p1(0)=0,p2(0)=√1+e1−e, (26)

if , then the solution is periodic with period , and the trajectory is an ellipse of eccentricity . Note that the gradient function must here be implemented carefully so as to be analytic for complex values of . Here, we define it using the following determination of the complex logarithm (analytic on the complex plane outside the negative real axis):

 ∀(x,y)∈R2 s.t. x+iy∉R−,L(x+iy)=log|x+iy|+2iarctan(yx+|x+iy|).

As a consequence, the analytic continuation of the function writes

 exp(−32L(x+iy)),

where and .

Here, as with the harmonic oscillator, we take as basic method the 2nd-order Strang splitting

 S[2]τ=φ[a]τ/2∘φ[b]τ∘φ[a]τ/2, (27)

where (respectively, ) corresponds to the exact solution obtained by integrating the kinetic energy (resp., potential energy ) in (25).

We take , integrate until the final time with Strang and the schemes obtained by the recursion (22) with for several time steps and compute the relative error in energy at the final time. Figure 2 show this error as a function of the inverse of the step size to illustrate the order of convergence: order 2 for Strang, order 4 for and order 6 for . For , and contrary to what happens to the harmonic oscillator, the observed numerical order is higher than expected, varying between 7 and 8. We do not have at present a theoretical explanation for this phenomenon. Figure 2 (right) depicts the time evolution of this error when the final time is .

4.3 The semi-linear reaction-diffusion equation of Fisher

Our third test-problem is the scalar equation in one-dimension

 ∂u(x,t)∂t=Δu(x,t)+F(u(x,t)), (28)

with periodic boundary conditions on the interval . Here is a nonlinear reaction term. For the purpose of testing our methods, we take Fisher’s potential [Sar15]

 F(u)=u(1−u)

as considered for example in [BCCM13].

The splitting corresponds here to solving, on the one hand, the linear equation with the Laplacian (as

), and on the other hand, the non-linear ordinary differential equation

 ∂u(x,t)∂t=u(x,t)(1−u(x,t)),

with initial condition , whose analytical solution is given by the well-defined (for small enough complex time ) formula

 u(x,t)=u0(x)+u0(x)(1−u0(x))(et−1)1+u0(x)(et−1).

Here we aim to solve Eq. (28) with periodic boundary conditions on the interval , and initial condition . Numerically, the interval is discretised on a uniform grid, i.e., , and is approximated by Fourier pseudo-spectral methods. In this way we construct a vector with components , . If we denote by the whole numerical solution computed by a certain integrator with step size from until the final time, and by the corresponding numerical solution computed by the same integrator with step size , then the quantity is a good indicator of the convergence order.

Numerical simulations were carried out in quadruple precision (with Intel Fortran) such that roundoff errors are suppressed. Figure 3 shows the successive errors , at final time , of the methods obtained with the sequence (22) with the Strang splitting as the basic method (left) and the fourth order scheme given by (23) (right) with different time steps . One can clearly observe that the convergence order matches the previous analysis with a slightly better performance for the highest order, analogously to the Kepler problem. Figure 4 shows the successive errors versus the number of basic integrators in each case.

4.4 The semi-linear complex Ginzburg–Landau equation

Our final test problem is the complex Ginzburg–Landau equation on the domain ,

 ∂u(x,t)∂t = αΔu(x,t)+εu(x,t)−β|u(x,t)|2u(x,t), (29)

with , and initial condition . Here, , and denote real coefficients. In physics, the Ginzburg-Landau appears in the mathematical theory used to model superconductivity. For a broad introduction to the rich dynamics of this equation, we refer to [vS95]. Here, we will use the values , and , for which plane wave solutions establish themselves quickly after a transient phase (see [WMC05]). In addition, we set

 u0(x)=0.8cosh(x−10)2+0.8cosh(x+10)2,

so that the solution can be represented in Figure 5.

To apply the composition methods presented in previous sections, it seems natural to split equation (29) as

 ∂u(x,t)∂t=(1+ic1)Δu(x,t)+εu(x,t), (30)

whose solution is