Continuous trigonometric collocation polynomial approximations with geometric and superconvergence analysis for efficiently solving semi-linear highly oscillatory hyperbolic systems

In this paper, based on the continuous collocation polynomial approximations, we derive and analyse a class of trigonometric collocation integrators for solving the highly oscillatory hyperbolic system. The symmetry, convergence and energy conservation of the continuous collocation polynomial approximations are rigorously analysed in details. Moreover, we also proved that the continuous collocation polynomial approximations could achieve at superconvergence by choosing suitable collocation points. Numerical experiments verify our theoretical analysis results, and demonstrate the remarkable superiority in comparison with the traditional temporal integration methods in the literature.


Introduction
The recent growth in extended Runge-Kutta-Nyström (ERKN) methods (see, e.g. [12, 16, 33-35, 38-40, 42]) for second-order multi-frequency highly oscillatory systems and the collocation methods for stiff differential equations (see, e.g. [16,18,32,36]) has let to the development of numerical integrators which systematically incorporate qualitative information transmitted from the underlying problem into their structure. In this paper, we focus on the design and analysis of temporal approximations through continuous trigonometric collocation polynomials for the semi-linear highly oscillatory hyperbolic system: where A ∈ ℝ d×d is supposed to be a symmetric positive definite or skew-Hermitian time-invariant matrix with the decomposition A = Ω 2 , f ∶ [t 0 , T] × ℝ d → ℝ d is a nonlinear function, and are the given initial values of the position u(t) and velocity ̇u(t) , respectively. This hyperbolic system plays an important role in a wide variety of practical application areas in science and engineering, including nonlinear optics, solid state physics and quantum field theory (see, e.g. [1,5,11,29]). As is known, with suitable spatial discretisation strategies, such as the finite difference method, the finite element method and the spectral method (see, e.g. [3, 14, 17, 20, 21, 26-28, 30, 31]), the semi-linear wave equations can be converted into the highly oscillatory system (1.1). Therefore, it will be significant to further develop new approaches for efficiently solving the semi-linear highly oscillatory system.
In recent works, based on the operator spectrum theory and the Duhamel's principle, the evolutionary partial differential equations (PDEs) have been formulated as abstract ordinary differential equations (ODEs) in suitable Banach spaces (see, e.g. [4,19,[22][23][24]). Therefore, in this paper, we will consider both wave equations in highdimensions and second-order highly oscillatory ordinary differential equations as an semi-linear hyperbolic system (1.1) in a suitable Banach space (X, ‖ ⋅ ‖) , and present a new temporal integration strategy and related theoretical analyses. If the nonlinear function f ∶ [t 0 , T] × X → X is continuous, then the solution of the semi-linear hyperbolic system (1.1) and its derivative satisfy the following variation-of-constants formula (see, e.g. [4, 22-24, 38, 40, 42]) Since A = Ω 2 , the operators 0 (t − t 0 ) 2 A and 1 (t − t 0 ) 2 A are defined as (1.3) 0 (t − t 0 ) 2 A = cos (t − t 0 )Ω and 1 (t − t 0 ) 2 A = sinc (t − t 0 )Ω , respectively. Furthermore, for any z ∈ [0, 1] , the solution of the hyperbolic system (1.1) and its derivative at time t = t n + zh ∈ [t n , t n + h] ⊆ [t 0 , T], n = 0, 1, 2, … can be expressed in where V = h 2 A and g t n + h = f t n + h, u(t n + h) . A considerable amount of attention has been paid to the variation-of-constants formula (1.2) or (1.4) to design effective and efficient numerical integrators for solving the highly oscillatory system (1.1). Approximating the nonlinear integrals appearing in formula (1.2) or (1.4) by suitable numerical quadrature formulas, a variety of numerical integrators have been proposed and analysed, including the Gautschi-type method of order two (see, e.g. [8-10, 12, 15]), the adapted Runge-Kutta-Nyström (ARKN) method (see, e.g. [7,41]), the extended Runge-Kutta-Nyström (ERKN) method (see, e.g. [6,25,34,35,38,40,42]), the arbitrarily high-order Birkhoff-Hermite (BH) method [22], and the trigonometric Fourier collocation (TFC) method [32,36]. The integrators derived by using the variation-of-constants formula (1.2) or (1.4), such as Gautschitype method, ARKN method, ERKN method, BH method and TFC method etc., are also termed trigonometric integrators. Therefore, the trigonometric integrators can exactly integrate the multi-frequency unperturbed highly oscillatory system ü(t) + Au(t) = 0 associated with (1.1). In other words, the trigonometric integrators can preserve the oscillatory structure [39] of the hyperbolic system (1.1), since the highly oscillatory behaviour of (1.1) is brought by the linear term Au. Moreover, it is well known that the collocation integrators not only provide a discrete set of approximations, but also a continuous approximation to the original solution (see [16]). By choosing suitable collocation nodes, the collocation methods could achieve a highest convergence order. In the recent research work (see [32,36]), Wang et al. took advantage of the shifted Legendre polynomials and the Lagrange polynomials to derive two kinds of trigonometric collocation methods for solving the system (1.1). However, these two kinds of trigonometric collocation methods are discrete approximations to original solutions. In addition, the local error bounds and the long-term convergence analysis are insufficient. Therefore, in this paper, using the continuous polynomial approximation, we will focus on combing the superior performance of the collocation integrators with the trigonometric integrators to design the continuous trigonometric collocation polynomial approximation for solving the highly oscillatory hyperbolic system (1.1), which could continuously approximate the original solution and preserve the oscillatory structure of the underlying system. Furthermore, the structure-preserving behaviour, single step error bounds, the long-term convergence and the superconvergence of the continuous trigonometric collocation polynomial approximations will be rigourously investigated in this work. This paper is organised as follows. In Sect. 2, we formulate the trigonometric collocation integrators based on the continuous polynomial approximation, and investigate the symmetry of the trigonometric collocation integrator. Section 3 is concerned with the rigorous analysis of single step error bounds, the long-term convergence, and the superconvergence of the continuous trigonometric collocation polynomial approximations. By choosing suitable collocation nodes, the proposed trigonometric collocation integrators can achieve a highest order of convergence. Furthermore, the long-term energy conservation is also analysed in this section. In Sect. 4, we display preliminary numerical results which demonstrate the advantages and efficiency of our new algorithms in comparison with the existing numerical methods in the literature. The last section is devoted to brief conclusions.

Formulation of the continuous collocation polynomial approximation
An essential feature of the collocation integrator is that we not only obtain a discrete numerical solution, but also a continuous polynomial approximation to the solution of the original system (Hairer et. al. [16]). In this section, using continuous collocation polynomial approximations, we will derive a class of trigonometric collocation integrators for solving the highly oscillatory hyperbolic system (1.1) and analyse its symmetry. Our analysis for the collocation integrator will be based on the abstract formulation of the hyperbolic equation (1.1) as an evolution equation in a Banach space (X, ‖ ⋅ ‖).

Definition 2.1
Let 0 ≤ c 1 < c 2 < ⋯ < c s ≤ 1 be distinct nodes. The continuous collocation polynomial approximation is to find a polynomial y(t) of degree s + 1 in the Banach space (X, ‖ ⋅ ‖) such that, Then the numerical solution of the hyperbolic system (1.1) is defined by u n+1 = y(t n + h) and ̇u n+1 =̇y(t n + h).
At first glance, it is difficult to seek such a polynomial in practice. Fortunately, following the Lagrange interpolation formula, the continuous collocation polynomial y(t) satisfies Here and in what follows, we always assume that z is in [0, 1]. Applying the variation-of-constants formula to (2.2), we obtain This seems like that the numerical solution of the continuous collocation polynomial approximations determined by Definition 2.1 could coincide with a class of ERKN integrators (see, e.g. [38,40,42]). The details will be confirmed in the following theorem. Proof It straightforwardly follows from the formula (2.3) that and a ij (V) are given by (2.5) and (2.6), respectively. The statement of the theorem is proved. ◻ Hairer et al. [16] have pointed out that all second order differential equation ü =f (u) written as ̇u = v and ̇v =f (u) are time reversible. The system (1.1) could be expressed as Therefore, the system (1.1) is clearly time reversible. Moreover, as is known, symmetric methods have excellent long-time behavior when solving reversible differential systems (see Hairer et al. [16]). The design and analysis of the symmetric integrators will be significant in the spirit of geometric integration. Thus, a most welcome feature of the continuous collocation polynomial approaches (2.1) is that it preserves the temporal symmetry by choosing suitable collocation nodes. We are now in a position to prove the symmetry of the trigonometric collocation integrators determined by Definition 2.1.
The symmetric nodes c 1 , c 2 , … , c s demonstrate an interesting fact that (2.16) Replacing all indicates i and j in (2.13), (2.14) and (2.16) by s − i + 1 and s − j + 1 , respectively, we find that the ERKN method (2.4) with weights (2.5) and (2.6) is identical to its adjoint method consisted of (2.13), (2.14) and (2.16). Therefore, the ERKN methods determined by (2.4) are symmetric. This completes the proof of the theorem. ◻

Main theoretical results
In this section, we state the error bounds of our continuous trigonometric collocation polynomial approximations for solving the hyperbolic system (1.1). The single step error bounds and the long-term convergence of continuous trigonometric collocation polynomial approximations (2.1) will be rigorously analysed. Our analysis will be based on the variation-of-constants formula (1.2) or (1.4) provided for the hyperbolic system (1.1). Before presenting our theoretical analysis, we need the following hypotheses on the regularity of u and the nonlinearity of f.
We also assume that all occurring derivatives are uniformly bounded.

Assumption 3.2 Let f(t, u)
be Lipschitz-continuous with respect to u, i.e., there exists a real number L such that

Remark 3.1 Under Assumptions 3.1, it should be noted that the composition:
is also Fréchet differentiable in a strip along the exact solution u(t). Therefore, when u ∈ C l [t 0 , T], X and g (l) ∈ L ∞ ([t 0 , T];X) for l = s, s + 1, … , 2s , the following assumptions are valid where K m , m = 0, 1, 2, … , s are constants and independent of h and A.

Remark 3.2
Here, we remark that the essence of the continuous collocation polynomial approximation is to find the polynomial y(t) in the Banach space (X, ‖ ⋅ ‖) . The (i) Under the local assumptions of y(t n ) = u(t n ),̇y(t n ) =u(t n ) , the single step error bounds of the continuous trigonometric collocation polynomial y(t) defined by Definition 2.1 satisfy: (ii) Furthermore, if the collocation polynomial y(t) satisfies y (m) (t n ) = u (m) (t n ) for m = 1, … , s , then the derivatives of y(t) satisfy the following estimations: Here, we point out that the constants C 1 and C 2 satisfy and are obviously independent of h and A.
Proof According to the definition of the continuous trigonometric collocation polynomial defined by Definition 2.1, it is clear that the polynomial y(t) satisfies Moreover, it is evident that the exact solution of (1.1) satisfies In the light of Assumption 3.1, E n (z, h) and its derivatives satisfy the estimations (i) Applying the variation-of-constants formula to the Eq. (3.6), we obtain Taking norms on both sides of Eqs. (3.7) and using Assumption 3.2 and the local assumptions y(t n ) = u(t n ),̇y(t n ) =̇u(t n ) , we obtain Thus, let C 1 ≥ 2K 0 in the above results, and then the proof of the first statement is complete. (ii) The proof of the second statement follows from applying the variation-ofconstants formula to i.e., Inserting the results (3.2) into the left sides of (3.11) and (3.12) and choosing the constant C 2 to satisfy we obtain This completes the proof. In what follows, we will investigate the longterm convergence of the continuous collocation polynomial approximations. For this purpose, we rewritten the Eq. (3.7) as the following compact form: where Ψ(z, , V) is defined as: that is, Ψ(z, , V) is the rotation operator It is easy to see that the norm of Ψ(z, , V) satisfies: This will be repeatedly used in our following theoretical analysis. For simplicity of analysis, we introduce the energy norm ||| u(t),̇u(t) ||| for the solutions u(t),̇u(t) ∈ L ∞ ([t 0 , T];X) , i.e., Moreover, we quote the following Gronwall inequality, which will play an important role in the convergence analysis.   Obviously, the constant C 4 is independent on n, h and A. Therefore, the theorem is confirmed. ◻ ||| e(t n + zh),̇e(t n + zh) ‖e(t n + zh)‖ ≤ C 4 h s and ‖̇e(t n + zh)‖ ≤ C 4 h s , ∀z ∈ [0, 1].

Superconvergence of continuous trigonometric collocation polynomial approximations
In the light of the results shown in Theorem 3.1 and Theorem 3.2, the term �z − c i � appears in the error bounds. This observation implies that by choosing suitable collocation nodes c i , such as the zeros of the shifted Legendre polynomial, can minimize the error bounds to achieve a higher convergence order. In our numerical experiments, we will choose the collocation nodes c i as the zeros of the shifted Legendre polynomials to generate the Gauss-type and Lobatto-type time integrators, respectively. The numerical results in Sect. 4 will show that the continuous trigonometric collocation polynomial approximation can be of superconvergence. In this subsection, we will present the rigorous theoretical analysis of this superconvergence. To this end, we consider y(t n + zh) as the solution of the following perturbed differential equation: where the remainder (t n + zh) will vanish at the chosen collocation points, i.e., (t n + c i h) = 0 for i = 1, 2, … , s. Subtracting the perturbed differential equation (3.21) from the Eq. (1.1) and applying the variation-of-constant formula to this result, we obtain after the linearisation that where the remainder r(t n + zh) is of magnitude O(‖u(t n + zh) − y(t n + zh)‖ 2 ) . Since the defect (t n + zh) vanishes at the collocation nodes c 1 , … , c s , the integral related to ‖ (t n + zh)‖ is equal to its quadrature error, i.e., where p is the algebraic order of the quadrature formula corresponding to the nodes c 1 , … , c s . In what follows, we will show that the precision of the trigonometric collocation integrators defined by Definition 2.1 could be improved by choosing suitable collocation nodes. (3.21) y(t n + zh) = −Ay(t n + zh) + f t n + zh, y(t n + zh) + (t n + zh), (3.28)

Moreover, the statements of Theorem 3.2 yield that the remainder r(t n + zh) is of the magnitude O(u(t n + zh) − ‖y(t n + zh)‖ 2 ) = O(h 2s
(3.32)      we can derive the fourth-order and sixth-order trigonometric collocation time integrators, and denoted as LTC2s4 and LTC3s6, respectively. In order to demonstrate the superiority of the proposed integrators, we select the following time integrators for comparison: • BH1: the symmetric Birkhoff-Hermite time integrator of order four derived in [22] • BH2: the symmetric Birkhoff-Hermite time integrator of order six derived in [22] • GAS2s4: the two-stage Gauss time integration method of order four presented in [16]; , , , , • GAS3s6: the three-stage Gauss time integration method of order six presented in [16]; • LIIIA3s4: the Labatto IIIA method of order four presented in [16]; • LIIIA4s6: the Labatto IIIA method of order six presented in [16]. • ERKN3s4: the three-stage explicit ERKN method of order four derived in [38]; • ERKN7s6: the seven-stage explicit ERKN method of order six presented in [33]; During the numerical experiments, it should be noted that the fixed-point iteration are used for all the implicit integrators. The iteration procedure will be stopped once the l ∞ norm of the difference between two successive approximations is smaller than the fixed error tolerance 10 −15 . Here, we also point out that if the error produced by a method is too large for some time stepsize h, then the corresponding point will not be plotted in the figure. Moreover, we should point out that the ERKN3s4 and ERKN7s6 methods are quite different with our proposed methods in this paper. The chosen ERKN3s4 and ERKN7s6 methods are derived by solving the order conditions (see [33,38]). Our derived time integrators are yielded by using the continuous collocation polynomial approximation.
All computations in the numerical experiments are carried out by using MAT-LAB 2016b on the computer Lenovo ThinkCentre M8300t (CPU: Intel (R) Core (TM) i7-2400 CPU @ 3.10 GHz, Memory: 8 GB, Os: Microsoft Windows 10 with 64 bit). Here, the dimensionless parameter is chosen as = 0.5 and 0.1, respectively. Moreover, we suppose that the Klein-Gordon equation (4.1) is equipped with the periodic boundary conditions. Therefore, the Fourier pseudo-spectral method will be used to discretise the spatial derivatives, and the Klein-Gordon equation (4.1) can be converted into the following form t) , and A = D 2 + I∕ 2 is an M × M matrix used to approximate the operator −Δ + 1∕ 2 . Here, D 2 is the secondorder symmetric semi-definite spectral differential matrix (see [14,17,26,27]), and I is the unit matrix. The discrete energy is given by where the discrete l p norms ‖ ⋅ ‖ p for p = 2, 4 are defined as the discrete inner product (u, v) is defined by Here, Δx = 2L∕M is the spatial stepsize. In order to test the long-term numerical behaviour of the proposed time integrators, we set L = 30, T = 100 and fixed the spatial mesh size M = 1024 for this problem. In addition, as is known, the exact solution of the Eq. (4.1) cannot be represented explicitly. Therefore, we use the posterior error, i.e., RE = ‖U(h;T) − U(h∕2;T)‖ 2 , to compute the convergence order.
In Tables 1 and 2, we test the numerical precision of the time integrators "GTC2s4", "GTC3s6", "LTC3s4" and "LTC4s6" with different . The numerical data indicate that the continuous trigonometric collocation polynomial approximations with Gauss points and Lobatto points can achieve 2sth order and (2s − 2) th order, respectively. The computational results in Tables 1 and 2 demonstrate that the temporal accuracy is completely consistent with our theoretical analysis. The proposed continuous trigonometric collocation polynomial approximations can be of superconvergence (see Theorem 3.4). The logarithms of the posterior errors, i.e., log 10 (RE) , are plotted in Figs. 1 and 2, in comparison with the classical time integrators "BH1", "BH2", "ERKN3s4", "ERKN7s6", "GAS2s4", "GAS3s6", "LIIIA3s4" and "LIIIA4s6". It can be observed from these two figures that the proposed continuous trigonometric collocation integrators are much more accurate than these traditional methods. The data in Tables 3 and  4 illustrate the precision of the discrete energy conservation, which is consistent with our theoretical analyse results in Theorem 3.5. The proposed temporal integrators have better long-term behaviour of energy conservation.

Problem 4.2
Consider the following sine-Gordon equations in two dimensions (see, e.g. [22,23]): where 2 is a dimensionless parameter. We also suppose that the two dimensional sine-Gordon equation (4.3) is subjected to periodic boundary conditions. The spatial derivatives are approximated by the two dimensional Fourier pseudo-spectral method. Therefore, the Eq. (4.3) can be converted to the following matrix form:  where D x F M y are the second-order Fourier spectral matrices (see, e.g. [13]) in the x direction and y direction, respectively. The discrete energy is given by where Δx = 2∕M x and Δy = 2∕M y are the spatial stepsizes. More details can be found in [13,26,27]. For this problem, we take the dimensionless parameter = 1∕20 and set the mesh sizes M x = M y = 256 . In Fig. 3 and 4, we present the simulation results and the corresponding contour plots at the time points t = 0, 20, 40, 50, 80 and 100 with h = 0.01 . The phenomena is termed circular ring solitons and obviously is periodicity. Figure 3 demonstrates the phenomena of the periodic oscillations and radiations of the circular ring solitons as time goes on. The CPU time required to reach T = 100 is 439.476393 seconds. The numerical data in Table 5 clearly indicate that the convergence order of the derived time integrators, which are consistent with our theoretical analysis results in Theorem 3.4. Figure 5 demonstrates that the time integrators "GTC2s4", "GTC3s6", "LTC3s4" and "LTC4s6" have better long-term computational behaviour than the chosen classical integrators. Table 6 shows the numerical precision of the long-term preservation of the discrete energy (4.4). The numerical results in Table 6 confirm the theoretical result in Theorem 3.5.

Problem 4.3
We consider the Duffing equation (see, e.g. [32]) with 0 ≤ k < . This is a Hamiltonian system with the Hamiltonian The analytic solution is given by where sn means the Jacobian elliptic function. We choose k = 0.03 and different frequencies = 10 and 20 which are similar to those in [32]. The problem is investigated over the interval [0, 1000] to verify the convergence order and the precision of the energy conservation for the constructed time integrators with different . The data in Tables 7 and 8 show the convergence order of the time integrators "GTC2s4", "GTC3s6", "LTC3s4" and "LTC4s6" with different . The tabular data implies that the continuous trigonometric collocation polynomial approximations with Gauss points and Lobatto points could achieve 2sth order and (2s − 2) th order, respectively. The logarithm of the global errors GE = ‖q N − q(1000)‖ 2 against different steps for Problem 4.3 are plotted in Figs. 6 and 7. In comparison with the "GAS2s4", "GAS3s6", "LIIIA3s4", "LIIIA4s6", "BH1" and "BH2" integrators, the proposed time integrators in this paper have much better accuracy and cost less CPU time. However, the proposed integrators and the chosen ERKN integrators have similar precision under the same CPU time. The numerical results in Tables 9 and 10 indicate the precision of the energy conservation.
In conclusion, the numerical results for Problem 4.1, Problem 4.2 and Problem 4.3 are consistent with our theoretical analysis results. The continuous collocation polynomial approximations can be superconvergence and have better numerical behaviour q(t) = sn( t, k∕ ), than the chosen classical integrators. Comparing with the classical implicit methods, i.e., "GAS2s4" method, "GAS3s6" method, "LIIIA3s4" method, "LIIIA4s6" method, "BH1" method and "BH2" method, the proposed time integrators "GTC2s4", "GTC3s6", "LTC3s4" and "LTC4s6" have better precision while cost less CPU time. However, comparing with the explicit ERKN3s4 method and ERKN7s6 method, our proposed time integrators have similar precision under the same CPU time. Moreover, the numerical results in [39] indicated that the ERKN3s4 and ERKN7s6 methods have     better numerical behaviour than the ARKN3s4, TFCr2 and TFCr3 methods. Therefore, we could confirm that the proposed time integrators in this paper would be better than the ARKN3s4, TFCr2 and TFCr3 methods. The Gautschi-type methods are time integrators with second order precision at most. We will not compare with the Gautschitype, ARKN and TFC methods in this paper.

Conclusion
Taking into account the superiority of the continuous collocation methods and the trigonometric integrators, we proposed and analysed the continuous trigonometric collocation polynomial approximations for the highly oscillatory hyperbolic system (1.1) in this paper. The derived trigonometric collocation integrators inherit the superconvergence of the classical continuous collocation methods and also can preserve the oscillatory structure of the underlying highly oscillatory systems. The resulting trigonometric collocation integrators were analysed in details for the local error bounds, long-term convergence, superconvergence, symmetry and long-term energy conservation. Furthermore, the remarkable numerical behaviour of the continuous trigonometric collocation polynomial approximations was demonstrated by the numerical experiments in comparison with the existing numerical methods in the literature.