{'fraction_non_alphanumeric': 0.046751863684771036, 'fraction_numerical': 0.017252396166134186, 'mean_word_length': 4.045853483790077, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A characterization of energy-preserving methods and the construction of parallel integrators for Hamiltonian systems

High order energy-preserving methods for Hamiltonian systems are presented. For this aim, an energy-preserving condition of continuous stage Runge--Kutta methods is proved. Order conditions are simplified and parallelizable conditions are also given. The computational cost of our high order methods is comparable to that of the average vector field method of order two.

1. Introduction. This paper is concerned with the numerical integration of Hamiltonian systems of the form where f (y) = S∇H(y) with a non-singular skew-symmetric constant matrix S and sufficiently differentiable Hamiltonian H : R N → R. Although the dimension N of the system is even, in the case of Hamiltonian systems, we do not assume this because the results are also applicable to Poisson systems with constant Poisson structure, where S is not necessarily non-singular and N need not be even. Several geometric properties are known in classical mechanics. For example, the exact flow is energy-preserving in the sense that d dt H(y) = ∇H(y) ⊤ dy dt = ∇H(y) ⊤ S∇H(y) = 0, as well as being symplectic. Numerical integrators exactly inheriting such properties are called geometric numerical integrators [11], and the most remarkable advantage of adopting geometric numerical integrators is that they often give qualitatively correct numerical solutions over an extremely long period of time. Although integrators inheriting such geometric properties as much as possible would be preferable, it is known to be impossible to construct methods inheriting both symplecticity and energy-preservation [7,23]. This is because if such methods exist the numerical flow coincides with the exact flow. In the last decades, methods inheriting one of these properties, i.e. symplectic methods and energy-preserving methods, have been developed. Although both symplectic methods and energy-preserving methods have their own advantages, much more attention has been devoted to symplectic methods. One of the most remarkable advantages of symplectic methods is that they also nearly preserve the energy without any drift. It should also be noted that there is a symplectic condition for Runge-Kutta (RK) methods [15,19,20]: where s is the stage number, and a ij and b i are the coefficients of RK methods. This condition is not only sufficient but also necessary under a certain assumption. In fact, symplectic implicit RK methods exist whose typical examples are Gauss methods. For more details, we refer the reader to [11]. While symplectic methods are useful for most Hamiltonian systems, it is sometimes mandatory to adopt energy-preserving methods. For example, errors in the Hamiltonian obtained by symplectic methods can sometimes be fatal for chaotic dynamics. Thus, it is strongly hoped that the study on energy-preserving methods reaches to the same maturity as symplectic methods. However, the study on energypreserving methods is more challenging and, in fact, newer than that on symplectic methods. One reason is that no RK method is energy-preserving for general Hamiltonian systems [5].
The projection method is a relatively simple method of constructing energypreserving integrators. The idea is very simple, but the implementation and the construction of the projection is specific to individual problems. An even more serious disadvantage of projection-based integrators is that the long term behaviour usually deteriorates. Compared with the projection method, the discrete gradient method, which was first proposed by Gonzalez [9] (see also McLachlan et al. [16]), is rather a more sophisticated and systematic approach.
After some pioneering studies, substantial progress was made in the last decade. From a theoretical viewpoint, Faou-Hairer-Pham [8] and Chartier-Faou-Murua [7] revealed that there exist energy-preserving B-series integrators other than the exact flow by showing an energy-preserving condition. Furthermore, concrete energy-preserving methods have also been developed, Quispel-McLaren proposed the average vector field (AVF) method [18]. This method is energy-preserving, a subclass of the discrete gradient method, a B-series method and of order two. A more recent interest is to find a high order energy-preserving method as an extension of the AVF method. Hairer proposed the AVF collocation method [10], and Brugnano et al. proposed the Hamiltonian boundary value method [1]. These methods are essentially the same methods, although the latter method has been developed mainly for polynomial Hamiltonian systems. These methods are based on so called continuous stage RK (CSRK) methods. It should be noted that roughly speaking the computational cost of the AVF collocation method is almost the same as the same order Gauss method.
The main aim of this paper is to construct a more efficient energy-preserving method than exists at present. To achieve this aim we focus our attention on a matrix M which characterizes CSRK methods. In terms of M , we consider (i) an energy preserving condition, (ii) order conditions and (iii) criteria for parallel implementation using real arithmetic. We then put together these three perspectives to construct high order energy-preserving methods which can be computed in parallel. We then derive concrete integrators of up to order six.
The paper is organized as follows. In Section 2, some existing energy-preserving methods are briefly reviewed with the introduction of CSRK methods. In Section 3, an energy-preserving condition of CSRK methods is discussed. In addition to the sufficient condition already given in [17], we show that the condition is also necessary under a small assumption. The condition is expressed in terms of the matrix M . In Section 4, order conditions of energy-preserving CSRK methods are characterized in terms of the matrix, where simplifying assumptions play an important role. In Section 5, parallel energy-preserving methods are constructed. We first consider a parallelizable condition of CSRK methods, and then construct parallel energy-preserving methods. Furthermore, we investigate fourth-order methods in detail. The derived fourth order methods can be implemented faster than the same order AVF collocation method, but an actual error of our methods per each time step becomes bigger than that of the existing method. We therefore have to discuss to what extent our method is efficient. Finally, concluding remarks are given in Section 6.
2. Energy-preserving continuous stage Runge-Kutta methods. In this section, we review some existing energy-preserving methods and continuous stage Runge-Kutta methods.
The discrete gradient method proposed by Gonzalez [9] (see also McLachlan et al. [16]) is a first systematic method for deriving energy-preserving integrators. The key idea of the method is to introduce a sophisticated approximation of the gradient. We first define a discrete gradient, a map ∇H : for any x, y ∈ R N , and then define an integrator by where h denotes the stepsize. Here, the condition (2.1) corresponding to the chain rule d dt H(x) = ∇H(x) ⊤ẋ is called the discrete chain rule, and the condition (2.2) requires the consistency. Due to the discrete chain rule and skew-symmetry of the matrix S, the energy-preservation for the discrete gradient method (2.3) can be easily verified: The discrete gradient is not generally unique and several approaches have been proposed. The average vector field (AVF) gradient proposed by Quispel-McLaren [18] ∇H(y 1 , y 0 ) = 1 0 ∇H(τ y 0 + (1 − τ )y 1 ) dτ is one of several possible approaches. In general, the AVF method reads The AVF method is a symmetric, second order and furthermore B-series method as will be shown later.
After the proposition of the AVF method, the AVF method was extended to higher orders. Hairer extended the AVF method by slightly modifying the idea of the standard collocation methods [10] (see also the Hamiltonian boundary value method [1]). In this paper, we call the extended method the AVF collocation method. As is the case with the standard collocation method which can be interpreted as a RK method, the AVF collocation method can be interpreted as a continuous stage Runge-Kutta (CSRK) method.
Definition 2.1 (Continuous stage Runge-Kutta methods). Let A τ,ζ be a polynomial in τ and ζ. Assume that A 0,ζ = 0. We denote by s the polynomial degree of A τ,ζ in τ . Let B ζ be defined by B ζ = A 1,ζ . We search for an s-degree polynomial Y τ (τ ∈ [0, 1]) and y 1 such that they satisfy A one-step method y 0 → y 1 is called an s-degree continuous stage Runge-Kutta (CSRK) method. A CSRK method is said to be consistent if 1 0 B ζ dζ = 1. Remark 1. CSRK methods were originally due to Butcher [2] (see also [3]), where A τ,ζ was not assumed to be a polynomial but to be bounded. In this paper, we focus on a polynomial A τ,ζ because this restriction simplifies discussion below and more general functions do not seem to be of greater benefit to our aim than polynomial functions.
In the above definition, A 0,ζ = 0 is assumed so that Y 0 = y 0 . Y τ can be seen as an approximation of y(t 0 + C τ h), where C τ = 1 0 A τ,ζ dζ. The final stage y 1 is also expressed by y 1 = Y 1 due to the relation B ζ = A 1,ζ . Since in this paper we shall focus on A τ,ζ which is polynomial of degree s in τ and s − 1 in ζ, we often express A τ,ζ as with a constant matrix M ∈ R s×s . In the case of the AVF collocation method, the polynomial A τ,ζ is given by and c 1 , . . . , c s are distinct, real numbers. It was proved in [10] that the AVF collocation method is energy-preserving for Hamiltonian systems independently of the c i values.
If the c i values are the zeros of the s-th shifted Legendre polynomial, then the method has order p = 2s. Here lists concrete expressions of A τ,ζ for s = 1, 2, 3: Note that in the case of s = 1 and A τ,ζ = τ , the corresponding CSRK method coincides with the AVF method (2.4).
The following theorem ensures the unique solvability of CSRK methods. Theorem 2.2 (Unique solvability [21]). Suppose that f satisfies a Lipschitz condition there exists a unique solution of (2.5). This theorem can be proved in a similar way to the proof of the unique solvability of RK methods (see, e.g. [3,12]) as mentioned in [21]. For the convenience of the reader, this proof is included below.
Proof. First, we note that the function space of s-degree polynomials with the property Y 0 = y 0 can be identified with R sN . We define a metric on R sN by and consider the map If the stepsize h satisfies hL max τ ∈[0,1] 1 0 |A τ,ζ | dζ < 1, the contraction mapping theorem ensures the existence and uniqueness of the fixed point of the map φ.
As pointed out in [10], the numerical solution of a CSRK method can be interpreted as a B-series. Below we introduce some notation to be needed in the subsequent sections. See [3,11] for more details on B-series.
Let T be the set of rooted trees We denote by |t| the number of vertices of a tree t. ]. The symmetry coefficient σ : T → R is defined recursively by For a map a : T ∪ {∅} → R, a formal series of the form is called a B-series. The exact time-h flow of a differential equationẏ = f (y) can be expressed as y(t 0 + h) = B(e, y(t 0 )), where the coefficients e are given by As is the case with Runge-Kutta methods, every CSRK method is a B-series integrator y 1 = B(φ, y 0 ), where the coefficients φ, called the elementary weights, are given by Here are the first four elementary weights: In addition to the set of trees T , we define a set of forests F and a set of free trees F T . The forest is an unordered finite collection of trees from T : An element of t ∈ F is written as t = t r1 1 · · · t rp p (t 1 , . . . , t p ∈ T ), and the sum of the vertices are denoted by |t|, e.g.
Next, we introduce the set of free trees F T . Note that there is a group of trees whose elements have topologically the same shapes, e.g. the shapes of , , and are the same. We define the set F T by using the Butcher product, for two trees u = [u 1 , . . . , u q ] ∈ T and v ∈ T , Multiple shift of root induces an equivalence. All trees in the same equivalent class has the same number of vertices, and have topologically the same shape. We call each equivalent tree a free tree and denote the set of free trees by F T . The canonical projection is denoted by π : T → F T . This map is not bijective, e.g. π −1 ( ) = { , , , }. For two equivalent trees u and v, we define a distance κ(u, v) by the minimum number of root shifts necessary to obtain u from v.
3. Necessary and Sufficient condition for energy-preservation. In order to construct efficient energy-preserving CSRK methods, it is mandatory to characterize an energy-preserving condition for CSRK methods. The sufficient condition was already given in [17]. Although one can proceed to the subsequent sections by using the condition, we also show that the sufficient condition is also necessary under a small assumption. The necessity indicates that it is almost impossible to find energy-preserving CSRK methods outside the framework presented in this paper.

Sufficient condition. Theorem 3.1 ([17]). A CSRK method is energy
For readers' convenience, we here show the proof of this theorem. Proof. If the matrix M is symmetric, it can be diagonalized by an orthogonal matrix P : M = P ⊤ ΛP , where Λ is a diagonal matrix whose elements are given by where ϕ(τ ) = [1, τ, . . . , τ s−1 ] ⊤ and (·) i denotes the i-th component. By using these notation, we have The last equality follows from the skew-symmetry of S. Remark 2. Tang and Sun have also provided essentially the same sufficient condition in a different manner [22].

Necessary condition.
The idea of proving the necessity is to use a necessary and sufficient condition of energy-preserving B-series methods shown by Chartier et al. [7]. We will also use some lemmas which are motivated by Celledoni et al. [6].
The necessary and sufficient condition of energy-preserving B-series methods is summarized in the following theorem.
where t is a member of π −1 (t).
As corollaries, we get the following properties. Lemma 3.3. If a consistent B-series method is energy-preserving, then it follows that Proof. The idea of the proof is similar to that of Theorem 2.2 in [6]. This theorem is proved by considering a free tree such that k vertices are connected to a common vertex, e.g.
Lemma 3.4. If a consistent CSRK method is energy-preserving, then it follows that Proof. The idea of the proof is similar to that of Lemma 2.3 in [6]. To prove this condition, we introduce double bush trees. As illustrated in the picture below, a double bush tree t p,q is a free tree having p leaves on the left side and r leaves on the right side (this figure is an example with p = 2 and q = 3).
It is clear that t p,p is superfluous, and t p,q = t q,p .
There are four distinct but equivalent rooted trees for t p,q : Here are examples for t 2,3 : Thus, it follows that The second and third are due to Theorem 3.3. The symmetry coefficients are calculated to be Substituting κ, a and σ into (3.1), we obtain the condition (3.3).
We consider an s × ∞ matrix Φ defined by We now show the following theorem. Proof. We see by contradiction that the consistency and the full rank assumption indicate B τ = C ′ τ . Under the assumption B τ = C ′ τ , we can rewrite the condition (3.3) as follows. Using the integration-by-parts formula, we have Similarly, Substituting these relations into the condition (3.3), we have Note that Since we assume that the matrix Φ is of full rank s, it follows that M = M ⊤ . Remark 3. Let us discuss the assumption about the rank of Φ. The technique of the proof of Theorem 3.5 is motivated by the proof of the necessary condition of symplectic RK methods, especially Theorem 7.10 in [11, Chapter VI], where a similar s×∞ matrix characterizes the reducibility of RK methods. In contrast to RK methods, however, we believe that for any consistent CSRK methods the rank of the matrix Φ is always of full rank s. This conjecture can be easily verified for monotonic C τ or the case the function space spanned by , v(0) = 0}, but the conjecture for more general cases has not been proved. We leave the proof of this conjecture to our future work.

Families of energy-preserving continuous stage
Runge-Kutta methods and order conditions. In this section, we characterize order conditions of energy-preserving CSRK methods in terms of the matrix M (2.6), based on the simplifying assumptions and construct several families of energy-preserving methods. The order of a CSRK method satisfying the simplifying assumption B(ρ), C(η) and D(ξ) is at least min(ρ, 2η + 2, η + ξ + 1). The following two theorems can further simplify the order conditions. Proof. We prove Recall that we always assume that C 0 = 1 0 A 0,ζ dζ = 0, and the consistency indicates C 1 = 1. Let so that ∂ ∂ζÃ τ,ζ = A τ,ζ andÃ τ,ζ =Ã ζ,τ . We start with (4.1). By using the integration-by-parts formula, we see that (4.1) is equivalent to Here, we note that Ã τ,ζ C k ζ 1 0 =Ã τ,1 =Ã 1,τ . By changing τ and ζ each other, we obtainÃ Differentiating the both sides with respect to ζ leads to which is equivalent to (4.2). Note that (4.3) is equivalent to (4.4) due to C 0 = 0.
In the following subsections, we characterize the order conditions of energypreserving CSRK methods in terms of the matrix M .  where all components of i k ∈ R s are zero except for the k-th component which is 1.
Proof. We obtain B ζ = 1 by substituting k = 1 into (4.5), and this indicates that C τ = τ . We consider an integrator satisfying the C(η) condition. Since this coincides with C k τ /k = τ k /k for any τ if and only if Remark 4. For the condition (4.5), it is possible to choose η = s. In this case, and thus M can be written as M = H −1 , where H is the Hilbert matrix defined by H ij = 1/(i+j −1). In this case, the CSRK method coincides with the AVF collocation method of order 2s.

Energy-preserving integrators with
If the C(2) assumption is satisfied, we have As illustrated above, energy-preserving integrators can be obtained even for the cases B ζ = 1. However, such cases require larger degrees, and thus are less practical than the cases B ζ = 1.

Parallel energy-preserving methods.
In this section, we construct new energy-preserving methods which can be implemented more efficiently than the AVF collocation method. In the context of RK methods, it is known that the computational cost of solving implicit RK methods can be reduced if a RK matrix A has only real, distinct eigenvalues (see, e.g. [4,13]). The key of the construction of new integrators is to apply similar idea to CSRK methods and put it together with the energy-preserving condition and order conditions.

Implementation of Runge-Kutta methods with real eigenvalues.
In general, in order to proceed with a RK process, we have to solve a system of nonlinear equations of size sN per each time step. Here we explain that if the matrix A has only real, distinct eigenvalues, the computational costs can be reduced, especially in a parallel architecture.
We consider solving an s-stage implicit RK method . . , s, In order to compute the internal stages Y 1 , . . . , Y s ∈ R N , we have to solve a system of nonlinear equations where ⊗ denotes the Kronecker product, I N ∈ R N ×N is an identity matrix and The size of the system (5.1) is sN . We usually apply the simplified Newton method to solve the nonlinear system (5.1), and then obtain the iteration formula where J 0 denotes the Jacobian matrix, i.e. J 0 = S∇ 2 H(y 0 ) for Hamiltonian systems.
Below, we show that the linear system (5.2) of size sN can be computed efficiently if the matrix A has only real, distinct eigenvalues.
Let M = (I sN − hA ⊗ J 0 ) and define M by M = (T −1 ⊗ I N )M (T ⊗ I N ). Then it follows that where the right hand side expresses a block diagonal matrix. Therefore, the formula where ρ l is calculated based on the following relation The key point of this idea is that the linear system (5.4) of size sN consists of s linear systems of size N .
Below we consider to what extent the computational cost is reduced. For both approaches, we equally pay N 2 operations to compute J 0 .
• Standard iteration based on (5.2). We compute the LU factors of (I sN − hA ⊗ J 0 ), where a small constant times (sN ) 3 operations are needed. Then we solve (5.2) iteratively until convergence, where O(s 2 N 2 ) operations are needed per each iteration. Remark 5. Note that the cost of the LU factorisation can be further reduced by making use of complex conjugate eigenvalues pairs [13]. For example, the cost of the LU factorisation for the 2-stage Gauss method can be reduced from

Efficient continuous stage Runge-Kutta methods and their implementation.
We apply the concept of the previous subsection for RK methods with real eigenvalues to CSRK methods.
First, let us express Y τ in (2.5) as where l i (τ ) is defined by We regard (2.5) as a system of nonlinear equations in terms of Y c1 , . . . , Y cs : Here we merely evaluate (2.5) at the nodes τ = c 1 , . . . , c s . As this was done for RK methods, the expression (5.6) is further re-expressed in a form Since the derivation of Φ with respect to Y becomes the application of the simplified Newton method for solving (5.7) gives the iteration formula where the matrix E ∈ R s×s is defined by Therefore, if the matrix E of a CSRK method has only real, distinct eigenvalues, the CSRK method can be implemented efficiently. The following theorem indicates that one does not have to be concerned about the c i values when evaluating the eigenvalues of E.
Theorem 5.1. The eigenvalues of the matrix E defined in (5.8) are independent of the c i values.
Proof. We define a Vandermonde matrix Note that the matrix V is non-singular, if the c i values are distinct. Below, we show that V −1 EV is independent of the c i values. We express A τ,ζ as A τ,ζ = s i=1 τ i P i (ζ) where P 1 , . . . , P s are polynomials in ζ. Since A ci,ζ = s k=1 c k i P k (ζ), the matrix E is expressed as Therefore, it follows that In the second equality, the relation s k=1 c j k l k (ζ) = ζ j is used. The expression of V −1 EV (5.9) can be further simplified. Since where the matrix M is defined in (2.6), the component of V −1 EV can be written as Hence, we have the following formula In the following subsections, we construct efficient fourth and sixth order energypreserving integrators by making use of the idea presented in this subsection.

Three-degree fourth order integrators.
5.3.1. Derivation. We derive fourth order energy-preserving integrators. Since the only fourth order CSRK method with degree two is the AVF method, we set the degree of CSRK methods to s = 3. We also assume that B ζ = 1. According to The characteristic polynomial of the matrix (5.10): It follows that ϕ(λ) = 0 has three real, distinct roots if − α 1 300 > 1 6 2 2/3 + 5 24 Calculating the range of α is easy, and exactly the same evaluation was done by Butcher-Imran [4] in the context of symplectic methods.

Long time behaviour.
Since the new method surely preserves the energy, one may expect a good long-time behaviour. As a criterion of the long-time behaviour, we here consider the conjugate-symplecticity [11].
A one-step method Φ h of order p is said to be conjugate-symplectic up to order p + r (r ≥ 0), if there exists such a change of coordinates z = χ(y) that is O(h p ) close to the identity and . The modified equation of such a method is Hamiltonian up to the terms of size O(h p+r ), and thus the behaviour of the method is like that of a symplectic integrator on intervals of length O(h −r ). See [11,Chapter VI.8] for more details.
In [10,14], it is proved that the AVF collocation method of order 4 (2s in general) is conjugate-symplectic up to order 6 (2s + 2), but not conjugate-symplectic up to a higher-order. It follows from Theorem 5.11 in Hairer-Zbinden [14] that the new method (of order 4) is conjugate-symplectic up to order 6 but not conjugatesymplectic up to a higher-order, independently of the choice of α. Therefore, one could expect that the new fourth-order method behaves like the fourth-order AVF collocation method in terms of conjugate-symplecticity.

Efficiency.
Although the computational cost in each time step of the derived integrators is reduced to about one-eighth (or one-fourth if the AVF collocation method is implemented by making use of complex conjugate eigenvalues pairs), the actual local truncation error of the new integrators is worse than the AVF collocation method. Therefore, we here carefully consider the efficiency of the proposed integrators by taking both truncation errors and computational costs into consideration. Below, we denote −α 1 /300 by θ. Table 1 shows coefficients of the elementary differentials with trees of order 5 for the B-series expansions of the exact solution, the two-degree fourth-order AVF collocation method and the proposed method. It is estimated that when we use the same stepsize the error of the proposed method is approximately (60θ + 1) times bigger than that of the fourth-order AVF collocation method. On the other hand, for large problems, the computation of the proposed method is 8 (or 4) times faster than the AVF collocation method per each time step. Thus it is fair to compare the local errors of the AVF collocation method with the stepsize 8h and the proposed method with stepsize h at t = t 0 + 8h so that the overall computational costs are almost the same. The former is denoted by err AVF and the latter by err New . Then the ratio of the errors at t = t 0 + 8h is roughly estimated to be 8 · err New err AVF ≈ 8 · (60θ + 1)h 5 (8h) 5 = 60θ + 1 4096 .
Therefore, the proposed method is more efficient than the fourth-order AVF collocation method (in the sense that the above ratio is less than 1) if 0.7770503941 < θ < 273 4 = 68.25.
Note that if we consider the new method is only 4 times faster than the AVF collocation method, this range of θ changes to 0.7770503941 < θ < 17 4 = 4.25, and thus the proposed method is still efficient. Table 1 Coefficients of the elementary differentials with trees of order 5 for the B-series expansions of the exact solution, the fourth-order AVF collocation method and the proposed method. We write α = 1 7 + α * , β = 1 8 + β * , γ = 1 9 + γ * .
Let ϕ(λ) = det(λI − N ). It is possible to find a parametric solution to the problem ϕ(0) = ϕ ′ (0) = 0 in terms of α * , β * and γ * . One way of doing this is explained in Step 1 below. Here are three steps that need to be carried out to find a real root solution to ϕ.
We choose the parameters a and c so that δ ∈ [0, 0.2144863035]. A possible choice is (a, c, δ) = ( 1 12 , 1, 1 8 ). 3. We add small perturbation to α or β or γ in the correct direction so that the two zero roots remain real. For example, reducing γ by 10 −10 gives the real roots −0.8831274189, −0.06591290801, 0.06591290335, 0.9185056830, 1.839542361 Remark 6. Although the five-degree sixth order integrator has been derived above, it is still difficult to give a necessary and sufficient condition for the parameters α, β and γ such that the matrix (5.10) has only real eigenvalues. Therefore, we here leave the discussion about the efficiency untouched.
6. Concluding remarks. In this paper, we have shown the energy-preserving condition of CSRK methods for Hamiltonian systems, characterized the order conditions and parallelizable condition in terms of the matrix M . As applications, we have derived new fourth and sixth order integrators which can be implemented more efficiently than the standard method, i.e. the AVF collocation method. In the derivation, the eigenvalues of the matrix E played an important role.
We note several directions for future work. From a theoretical viewpoint, it is hoped that the conjecture discussed in Remark 3 will be solved. Considering other applications of the energy-preserving condition of CSRK methods would be interesting. We are currently trying to give a more systematic approach to deriving higher order efficient energy-preserving integrators.