Legendre spectral collocation method for second-order nonlinearordinary/partial differential equations

We propose an efficient Legendre-Gauss collocation algorithm for 
second-order nonlinear ordinary differential equations (ODEs). We 
also design a Legendre-Gauss-type collocation 
 algorithm for time-dependent second-order nonlinear partial 
differential equations (PDEs), which can be implemented in a 
synchronous parallel fashion. Numerical results indicate the high 
accuracy and effectiveness of the suggested algorithms.


1.
Introduction. During the past decades, there has been a great deal of interest in the research of numerical integration for initial value problems associated with second-order ODEs. Such problems often arise in different fields of applied sciences and engineering. Particularly, after spatial semi-discretization, a large class of nonlinear wave equations such as Klein-Gorden and sine-Gorden equations, are usually reduced to certain systems of second-order ODEs. There are many numerical methods for solving second-order ODEs, and among those the Runge-Kutta method is now accepted as the most effective one, the interested reader may refer to [6,18,19,23,30] and the references therein.
As we know, spectral method has been widely used for numerical PDEs with smooth solutions due to its high-order accuracy (cf. [4,5,7,8,11,13,14,27,28]). For time-dependent PDEs, one usually uses the spectral scheme in space, and the difference scheme in time. However, this tactic results in an unbalanced scheme, i.e., the approximate solution has infinite accuracy in space and finite accuracy in time. In recent years, the spectral method in time has been developed by some researchers. For instance, TAL-Ezer [31,32] presented spectral methods in time using polynomial approximation of the evolution operator in Chebyshev least-square sense for linear hyperbolic and parabolic PDEs with periodic boundary conditions. Glenn et al. [12] proposed spectral methods in time for nonlinear PDEs with periodic boundary conditions. Bar-Yoseph et al. [3] and Zrahia et al. [40] suggested the space-time spectral element methods for nonlinear advection-diffusion problems and second-order hyperbolic equations. Shen and Wang [29] proposed a space-time spectral method using Fourier-like bases in space and a dual-Petrov-Legendre-Galerkin formulation in time for parabolic equations. Tang and Ma [33,34] developed Legendre spectral methods in time and space for linear hyperbolic and parabolic PDEs. Moreover, the high-order finite element methods were also developed for time discretization. For instance, Babuška and Janik [1,2] analyzed the p and hp finite element method in time and space for parabolic equations. Schötzau and Schwab [25,26] applied the hp-discontinuous Galerkin finite element method for time discretization of parabolic problems and initial value problems. Wihler [36] developed the continuous hp-Galerkin finite element time-stepping method for nonlinear ODEs. In the meanwhile, spectral collocation methods have become increasingly popular for solving nonlinear ODEs. Recently, Guo et al. [15,16,17,35] developed several Legendre-Gauss-type collocation methods for nonlinear ODEs based on Legendre polynomial expansions. In [21], Kanyamee and Zhang conducted a systematic comparison of a Legendre(Chebyshev)-Gauss-Lobatto collocation methods with some symplectic methods in solving Hamiltonian dynamical systems. In [38], we proposed some new Legendre-Gauss-type collocation algorithms for the first-order ODEs and applied them to time-dependent PDEs.
The aim of the present paper is to develop some new Legendre-Gauss-type collocation algorithms for the second-order ordinary/partial differential equations. We start with the initial value problem of second-order ODE: where f is a given function, V 0 and U 0 are the initial data. Such problem can be transformed into a system of two first-order differential equations, and then solve them numerically. However, for saving work, it seems reasonable to solve them directly. For this purpose, we first introduce a single interval Legendre-Gauss collocation scheme for (1), i.e., the solution is directly approximated by a finite Legendre polynomial series, and the numerical scheme is collocated at the Legendre-Gauss points. We construct a stable and efficient algorithm for numerical implementation of the single interval scheme. For a more effective implementation, we also propose a multi-interval Legendre-Gauss collocation scheme due to the following considerations: • For large T , it is necessary to partition the solution interval (0, T ) and solve the subsystems successively. The resulting subsystems can be solved more efficiently with a modest number of unknowns. Hence, the multi-interval scheme can be implemented efficiently and economically. • The multi-interval scheme provides us sufficient flexibility with respect to variable time steps and local approximation orders, e.g., it is able to cope with the loss of regularity of solutions at t = 0.
It is worth noting that Guo and Yan [17] have proposed a Legendre-Gauss collocation algorithm for second-order ODEs. However, our new algorithms are much simpler, and require much less computational cost compared with that in [17]. Numerical experiments show that the propesed algorithms are very effective for ODEs with oscillating, stiff, singular and long-time behaviors.
We then consider the time-dependent nonlinear PDEs, and propose a Legendre-Gauss-type collocation scheme for space-time discretization. More specifically, we use the above Legendre-Gauss collocation method in time and the Legendre-Gauss-Lobatto collocation method in space. The algorithm can be implemented in a parallel fashion. Numerical examples show that our algorithms possess an excellent stability and the spectral accuracy in both time and space.
The paper is organized as follows. In Section 2, we present the single interval and multi-interval Legendre-Gauss collocation algorithms for second-order nonlinear ODEs, and provide some numerical examples with oscillating, singular, stiff and long-time behaviors to exhibit the effectiveness of the suggested algorithms. In Section 3, we design the Legendre-Gauss-type collocation scheme for time-dependent second-order nonlinear PDEs. We also perform some numerical experiments to demonstrate the spectral accuracy of the proposed algorithm. Finally, we end with some concluding remarks in Section 4.
2. Collocation algorithms with time discretization for ODEs. In this section, we shall present an efficient Legendre-Gauss collocation algorithm with time discretization for problem (1). 1] be the standard Legendre polynomial of degree l. Clearly, L l (x) is the eigenfunction of the singular Strum-Liouville problem We recall that the shifted Legendre polynomial L T,l (t) is defined by (cf. [15]) In particular, Due to properties of the standard Legendre polynomials, there hold (cf. [15]) It can be easily verified that the shifted Legendre polynomials satisfy the orthogonality relation where δ l,m is the Kronecker symbol. Thus for any v ∈ L 2 (0, T ), there holds 302 LIJUN YI AND ZHONGQING WANG By (2) we can deduce that Moreover, according to some classical properties of Jacobi polynomials, we also derive readily that Next, let t N j , 0 ≤ j ≤ N be the nodes of the standard Legendre-Gauss interpolation on the interval (−1, 1), and ω N j , 0 ≤ j ≤ N be the corresponding Christoffel numbers. The nodes of the shifted Legendre-Gauss interpolation on the interval (0, T ) are the zeros of L T,N +1 (t), denoted by t N T,j , 0 ≤ j ≤ N . According to [15], we can verify readily that t N T,j = T 2 (t N j + 1), and the corresponding Christoffel numbers are ω N T,j = T 2 ω N j , 0 ≤ j ≤ N . We denote by P N (0, T ) the set of polynomials of degree at most N . Due to the property of the standard Legendre-Gauss quadrature, there holds for any φ ∈ P 2N +1 (0, T ) (cf. [15]), (10) Let (u, v) T and v T be the inner product and norm of the space L 2 (0, T ), respectively, The discrete inner product and the discrete norm are defined as Due to (10), for any φψ ∈ P 2N +1 (0, T ) and ϕ ∈ P N (0, T ), we have (cf. [15]) The shifted Legendre-Gauss interpolation operator Thanks to (11), we have for any φ ∈ P N +1 (0, T ) that The shifted Legendre-Gauss interpolation I T,N v(t) can be expanded as From (7) and (12) we obtain 2.2. The single interval scheme. We first introduce the single interval Legendre-Gauss collocation scheme for solving (1). It is to seek u N (t) ∈ P N +2 (0, T ) such that It is noted that Guo and Yan [17] derived a Legendre-Gauss collocation algorithm for (15) and analyzed the convergence of (15). Here, we shall present a new efficient algorithm for (15), which is much faster and easier to implement. We now describe this process. In view of (15), we have Next, let Clearly We also denote Then, by using (14) and (11) we get (20) Due to (16)- (19), (3) and (5), we have (cf. [38]) Furthermore, due to (3) and (5), a direct computation yields According to (9), {∂ 2 t L T,k (t)} k≥2 are mutually orthogonal polynomials. Hence, we compare the expansion coefficients in terms of ∂ 2 t L T,k (t) to obtain that where { f k } N k=0 is given by (20). On the other hand, substituting t = 0 into (18) and using the fact ∂ t L T,k (0) = 1 T (−1) k−1 k(k + 1), we obtain from (15) that Moreover, taking t = 0 in (17), using (15), (24) and the fact L T,k (0) = (−1) k , we obtain Let U = ( u 0 , u 1 , · · · , u N +2 ) T . Then, in view of (17), (18) and (20), we can rewrite (23)-(25) into a nonlinear system: where an obvious expression. We can employ a simple iterative process to evaluate the expansion coefficients U in (26), as stated below.
The Simple Iterative Algorithm: Provide an initial guess for the coefficients (17) and (18); k=0 by (23), (24) and (25); Renew the data of u N (t) and ∂ t u N (t) at {t N T,j } N j=0 by (17) and (18); Repeat steps -; k=0 u k by (17). Note that the simple iterative algorithm can be described as with U 0 be an initial guess. This iterative method is also called the fixed-point iteration. This algorithm is much simpler, faster and easier to implement than the Newton-Raphson iterative process adopted in [17]. Numerical experiments presented in Section 2.4 confirm the effectiveness and high accuracy of this algorithm.
Remark 1. The convergence of the fixed-point iteration has been studied extensively. We list the standard result as follows (cf. [22]): Let Ω be a closed subset of R N +2 and suppose that F : Ω → R N +2 is Lipschitz continuous on Ω with a constant γ < 1 such that for all x, y ∈ Ω and F (x) ∈ Ω for all x ∈ Ω. Then there exists a unique fixed point x * ∈ Ω of F , and the iteration defined by (27) converges q-linearly to x * with q-factor γ for all initial iterates x 0 ∈ Ω.
2.3. The multi-interval scheme. We now introduce the multi-interval scheme.
Let M and N m , 1 ≤ m ≤ M be any positive integers. We first divide the interval Firstly, replacing T and N by τ 1 and N 1 in (15) and all other formulas in Subsection 2.2, we can obtain the numerical solution u N1 Secondly, we evaluate the numerical solutions u Nm m (t) ∈ P Nm+2 (0, τ m ), 2 ≤ m ≤ M, step by step. Finally, the global numerical solution of (1) is given by More precisely, we denote by t Nm τm,k the nodes of the shifted Legendre-Gauss interpolation on the interval (0, τ m ). The multi-interval scheme for (1) where ∂ t u Nm−1 m−1 (τ m−1 ) can be evaluated by (18).
In view of (29) and (30), we observe that the local numerical solution u Nm m (t) is actually an approximation to the local exact solution U m (t), with the approximate initial data ∂ t u Nm . The proposed method can also be applied to a system of second-order ODEs.

LEGENDRE SPECTRAL COLLOCATION METHOD 307
The multi-interval Legendre-Gauss collocation scheme for (31) The global numerical solution of (31), denoted by − → u N (t), is given by In the following, we take the initial guess u k = 0, 0 ≤ k ≤ N +2, and use uniform time step-size τ m = τ and uniform mode N m = N in the multi-interval scheme for convenience, if there are no special instructions provided.
The exact solution of (34) is given by U (t) = (t + 1) 3 2 + 2 sin(5t), which oscillates and grows to infinity as t → ∞. We first use the single interval scheme (15) to resolve the problem (34). In Figure  1, we plot the point-wise absolute errors log 10 E N p (T ) at T = 0.5, 1, 2 with various values of N, by using the simple iterative algorithm. Clearly, the numerical errors decay exponentially as N increases. In Figure 2, we also present the numerical errors versus the number of iteration steps for the single interval scheme (15).
We next use the multi-interval scheme (29), combined with the simple iterative algorithm at each time step to resolve (34) numerically. In Figure 3, we list the numerical errors of the multi-interval algorithm for problem (34) with moderate τ and N . Clearly, this algorithm is stable and accurate even for large T .
Note that the simple iterative algorithm is still valid if we collocate the numerical scheme at the Legendre-Gauss-Radau points (with the right-end point). In Figure  4, we make a comparison between the Legendre-Gauss collocation method and the Legendre-Gauss-Radau collocation method. It is observed that the former method performs slightly better than the latter one.  The second example is the Duffing equation (cf. [9]): with ω > 0, 0 ≤ k < ω. The analytic solution of this problem is given by U (t) = sn(ωt; k/ω), which represents a periodic motion in terms of the Jacobian elliptic function. We first use the single interval scheme (15) to resolve the problem (35). In Figure  5, we plot the point-wise absolute errors log 10 E N p (T ) at T = 1 with various values of N for different choices of the parameters ω and k. Clearly, the numerical errors     decay exponentially as N increases. In Figure 6, we present the numerical errors versus the number of iteration steps for the single interval scheme (15). In addition, we also list in Figure 7 the number of iteration steps with different initial guesses. It can be seen that the number of iteration steps is not sensitive to initialization.
We next use the multi-interval scheme (29), combined with the simple iterative algorithm at each time step to resolve (35) with ω = 5, k = 0.03. In Figure 8, we list the numerical errors of the multi-interval algorithm for problem (35) with moderate τ and N . Clearly, this algorithm is stable and accurate even for very large T .

Singular solution.
Consider the following problem: where the exact solution is given by U (t) = t r , and r > 1 is a non-integer. Obviously, the solution has finite regularity for non-integer r.
We first use the multi-interval scheme (29) with uniform steps, combined with the simple iterative algorithm at each step, to resolve (36) numerically. Denote by Err the maximum of the absolute errors: In Figure 9, we plot the absolute errors against various N with τ = 0.1 and r = 1.5, 2.5, 3.5, respectively, which shows that the proposed method provides an accurate approximation.
As is well known, to resolve the singular behavior of the solution more efficiently, it would be more appropriate to use geometrically refined steps and linearly increasing degree vectors (cf. [36]). More specifically, we take • Mesh nodes: • Degree of polynomials: where [µm] denotes the largest integer ≤ µm. We next use the multi-interval scheme (29) with geometrically refined steps and linearly increasing degree vectors, combined with the simple iterative algorithm at each step, to resolve (36) numerically. For our computations, we choose µ = 1.5 and γ = 0.2 (we refer to [25] and the references therein for further details on the optimal choice of these parameters). We denote by "DOF" the total degrees of freedom. The errors are plotted in Figure 10 with r = 1.5, 2.5, 3.5, respectively. It was shown that exponential convergence rates were achieved with respect to DOF 1/2 .
We first use the the single interval scheme (15) to resolve the problem (38). In Figure 11, we plot the point-wise absolute errors log 10 E N p (T ) at T = 1 with various values of N for different choices of the parameter ω. Clearly, the numerical errors decay exponentially as N increases.
We next use the multi-interval scheme (29), combined with the simple iterative algorithm at each time step to resolve (38) with different ω. In Figure 12, we list the numerical errors of the multi-interval algorithm for problem (38) with τ = 0.5 and N = 10. Clearly, this algorithm is stable and accurate even for very large T .

2.4.4.
Long-time calculation. Consider the two-body problem (see [37]): where e ∈ [0, 1) is the (constant) eccentricity of elliptical orbit. The Hamiltonian function of the system is given by For description of numerical errors, we denote by q N 1 (t) and q N 2 (t) the numerical solutions, by H N (t) the numerical energy of the Hamiltonian, and by E N (t) the energy error at t, i.e., where H 0 is the initial energy of the Hamiltonian.  Figure 14. Numerical orbit (q N 1 (t), q N 2 (t)) of problem (39).
We use the multi-interval scheme (32) to solve (39) with e = 0.2. We take the initial guesses q N 1 (t) ≡ 1 and q N 2 (t) ≡ 1 at all collocation points for simplicity. In Figure 13, we list the point-wise energy errors E N (t) for t ≤ 10 7 , with uniform τ = 1 and N = 15. It can be observed that our algorithm is stable and provides accurate numerical results. In Figure 14, we plot the numerical orbit (q N 1 (t), q N 2 (t)) for 0 ≤ t ≤ 10 7 , which is virtually indistinguishable with the exact orbit of movement governed by the Hamiltonian system (39).
To exhibit the high efficiency, we also compare our simple iterative algorithm with the Newton-Raphson (N-R) iterative algorithm used in [17] for the problem (39) with e = 0.2. We first test convergence rates of both algorithms for the single interval scheme (15) with different T and N . It can be seen from Figure 15 that our simple iterative algorithm has almost the same convergence rates as the Newton-Raphson iterative algorithm. We next use the multi-interval scheme (32) with τ m ≡ 1 and N m ≡ 15 to solve (39) for 0 ≤ t ≤ 100000. In Table 1, we list the energy errors E N (t) and the corresponding CPU elapsed time (CPUT) at different time t. It can be observed that our algorithm costs less computational time, but provides slightly better approximation results than that in [17]. Another example is the Hénon-Heiles model (see Haier et al. [20]): and the Hamiltonian function of the system is given by We use the multi-interval scheme (32) to solve (41). In Figure 16, we list the point-wise energy errors E N (t) for t ≤ 10 7 , with uniform τ = 1 and N = 15. It can be observed that our algorithm is stable and provides accurate numerical results. In Figure 17, we also plot the numerical orbit (q N 1 (t), q N 2 (t)) of problem (41).
3. Collocation algorithms with space-time discretization for PDEs. In this section, we present the Legendre-Gauss-type collocation algorithms for space-time discretization of PDEs. 3.1. Discretization in space and time. Consider the following nonlinear hyperbolic PDE: with suitable consistent initial and boundary conditions. We shall use the Legendre-Gauss-Lobatto collocation method for spatial discretization, and the Legendre-Gauss collocation method for temporal discretization.
We first describe the semi-discretization in space. Let {x k } M k=0 be the Legendre-Gauss-Lobatto points with x 0 = −1 and x M = 1, and denote by {φ k } M k=0 the Lagrange basis functions associated with the points {x k } M k=0 . Then the Legendre-Gauss-Lobatto collocation method in space is to find U M (x, t) ∈ P M (−1, 1), i.e., Clearly Inserting (43) into (44), we further obtain We next describe the semi-discretization of (45) in time. For this purpose, we introduce the first-order differentiation matrix and the second-order differentiation matrix It was shown in [28] that D (2) = D (1) D (1) . Hence, we can rewrite the ODEs (45) in the following form: We use the multi-interval scheme (32), combined with the simple iterative algorithm at each time step, to resolve (49) numerically.
We measure the numerical errors by the maximum norm: where u N k (t) is the numerical solution of u k (t), 1 ≤ k ≤ M − 1, and l=0û k,l L T,l (t).   In Figures 18 and 19, we test the exact solutions u(x, t) = cos(t) sin(πx) and u(x, t) = e −t sin(πx), respectively, and list the numerical errors log 10 E max (T ) with τ = 0.1. We observe that the suggested algorithm provides accurate and stable numerical results. In Figures 20 and 21, the rates of convergence in both space and time are investigated with fixed N or M , respectively. It can be seen from Figures  20 and 21 that the errors decrease rapidly as M or N increases.

3.2.2.
The Klein-Gordon-Zakharov system. Consider the Klein-Gordon-Zakharov (K-G-Z) equations (cf. [24]): with the solutions The initial and boundary conditions are extracted from the exact solutions (52) and (53). The K-G-Z equations are used to describe the interaction of the Langmuir wave and the ion acoustic wave in plasma.
To solve the problem (51) more efficiently, we shall make minor changes of the numerical scheme as described in Subsection 3.1 such that it can be applicable to the case of x ∈ [−20, 20] in the space direction. For this purpose, we only need to extend the standard Lagrange basis functions based on the Legendre-Gauss-Lobatto points from the interval [−1, 1] to [−20, 20] by using a simple scaling, and then replace the differential matrix D (1) and D (2) as defined in (46)   We use the multi-interval scheme (32), combined with the simple iterative algorithm at each time step, to resolve (51) numerically. For convenience, we denote by U and V the numerical approximations of u and v, respectively, and
In Figures 22 and 23, we plot the curves of the solitary waves at various times computed by the multi-interval scheme (32) with τ = 0.01, M = 300 and N = 10.  It is shown that the waves at t = 5 and t = 10 agree with the one at t = 0 quite well, which implies the high accuracy of our scheme. In Figure 24, we investigate the rate of convergence in space with τ = 0.001, N = 2 and T = 1. It can be seen that the error decays rapidly as M increases. In Figure 25, we plot the numerical errors at T = 10 with uniform τ = 0.001, N = 2 and different M . Again, we find that the suggested algorithm provides accurate and stable numerical results.
We also make a comparison of our scheme with the numerical scheme proposed in [39], where the Legendre spectral Galerkin method is used for space discretization and the Crank-Nicolson finite difference method is used for time discretization. We use the multi-interval scheme (32) with uniform time step τ = 0.001 and N = 2, while one in [39] uses the time step τ = 0.0001. The numerical errors under the maximum norm with different M are listed in Table 2. Obviously, our scheme is more accurate.
4. Concluding remarks. In this paper, we proposed an efficient and accurate algorithm of the Legendre-Gauss collocation method for initial value problems of second-order ODEs, and applied this algorithm to the multi-interval discretization. We also designed an efficient Legendre-Gauss-type collocation algorithm for timedependent second-order nonlinear PDEs. These approaches have several fascinating features: • The single interval Legendre-Gauss collocation algorithm is easy to implement and possesses the spectral accuracy.
• The multi-interval Legendre-Gauss collocation algorithm provides us sufficient flexibility to adapt to the evolutionary process of solutions, with which we can evaluate the numerical solutions step by step with moderate mode N and stepsize τ . This approach not only simplifies actual calculation, but also keeps the global spectral accuracy. Hence, it is particularly suitable for long time calculations. • The proposed Legendre-Gauss-type collocation algorithm for time-dependent second-order nonlinear PDEs possesses an excellent numerical stability and high accuracy in both time and space. Particularly, the algorithm can be implemented in a synchronous parallel fashion.
The numerical results demonstrated the effectiveness and spectral accuracy of the proposed algorithms. These algorithms presented in this paper can be easily extended to a large amount of time-dependent nonlinear PDEs with different boundary conditions.