NUMERICAL ANALYSIS OF LEGENDRE-GAUSS-RADAU AND LEGENDRE-GAUSS COLLOCATION METHODS

Abstract

In this paper, we provide numerical analysis of so-called Legendre Gauss-Radau and Legendre-Gauss collocation methods for ordinary differential equations. After recasting these collocation methods as Runge-Kutta methods, we prove that the Legendre-Gauss collocation method is equivalent to the well-known Gauss method, while the Legendre-Gauss-Radau collocation method does not belong to the classes of Radau IA or Radau IIA methods in the Runge-Kutta literature. Making use of the well-established theory of Runge-Kutta methods, we study stability and accuracy of the Legendre-Gauss-Radau collocation method. Numerical experiments are conducted to confirm our theoretical results on the accuracy and numerical stability of the Legendre-Gauss-Radau collocation method, and compare Legendre-Gauss collocation method with the Gauss method.

1. Introduction

Consider the initial value problem for ordinary differential equations (ODEs)

We assume that f ∈ ℝ is sufficiently differentiable on [0, T] and satisfies a global Lipschitz condition, i.e., there is a constant L ≥ 0 such that

Under this assumption, Problem (1) has a unique solution y(t) defined on t ∈ [t0, T].

The method of collocation for the solution of differential equations is given in a general form by Collatz [7]. The principle is very simple: the solution of the differential equation is represented as a linear combination of known functions and the unknown coefficients in this representation are found by satisfying the associated conditions, and the differential equation at an appropriate set of points in the range of interest. Usually, only polynomial forms of solution are considered and the points are taken inside the range of interest. The collocation method also gives a natural way of generating dense polynomial output.

The collocation methods for ODEs are equivalent to a sub-set of Runge-Kutta methods, in the sense that the final values at the end of the step are algebraically the same. The sub-set of implicit Runge-Kutta methods which are equivalent to collocation methods appear to include most of those suggested for practical use. They include methods based on the classical quadrature points such as Gauss-Legendre, Radau where one end point is used, and Lobatto where both end points are used [6, 12, 13, 15]. Wright [21] discussed the relationships between various one-step methods for ODEs and gave a unified treatment of the stability properties of the methods. Butcher [5] has presented an investigation on the special features of the Runge-Kutta collocation methods based on the zeros of various orthogonal polynomials. Barrio [1] investigated the A-stability of the symmetric Runge-Kutta methods based on collocation at the zeros or extrema of the ultraspherical (Gegenbauer) polynomials. Wang and Guo [20] and Guo andWang [11] proposed Legendre-Gauss-Radau collocation method (LGRCM) and Legendre-Gauss collocation method (LGCM) for ODEs and proved that both methods possess the spectral accuracy. The LGRCM and LGRCM also can be regarded as a differential quadrature method in the time domain [2, 3, 4, 8, 9, 10, 17]. Recently, Kiliçman, Hashim, Tavassoli Kajani and Maleki [14] proposed a rational second kind Chebyshev pseudo-spectral method for the solution of the Thomas-Fermi equation over an infinite interval, which is a nonlinear singular ordinary differential equation. Tavassoli Kajani, Maleki and Kiliçman [19] constructed a multiple-step legendre-gauss collocation method for solving volterra’s population growth model which is a nonlinear integro-differential equation. Tohidi and Kiliçman [18] proposed a collocation method based on the Bernoulli operational matrix for solving nonlinear boundary value problems which arise from the problems in calculus of variation.

The manuscript is organised as follows. In Section 2, the LGRCM is briefly reviewed and recast as a Runge-Kutta method. By making use of the theory of the Runge-Kutta methods, the accuracy, A-stability and algebraic stability properties of the LGRCM are discussed as well. In Section 3, the LGCM is briefly reviewed and is proved equivalent to the well-known Gauss method. Numerical examples are then given in Section 4.

 

2. Properties of LGRCM

Let Pn(0, T) be the set of polynomials of degree at most n on [0, T]. It is well-known that a collocation method for solving (1) is to seek a polynomial uN(t) ∈ PN(0, T) such that

where tk, k = 1, 2, ⋯ , N, are collocation points.

2.1. Brief review. Recently, Wang and Guo [20] proposed to seek a collocation solution ∈ PN(0, T) such that

where the nodes = 0, 1, ⋯ , N, are the zeros of LN(t) + LN+1(t) on [0, T]. Here Ll(t) is the shifted Legendre polynomial of degree l on [0, T] defined by

where Pl(t) is the standard Legendre polynomial of degree l on [−1, 1]. The set is a complete L2(0, T)-orthogonal system and thus the collocation solution can expanded as

For any ϕ ∈ P2N(0, T), it follows from the property of the standard Legendre-Gauss-Radau quadrature that with the corresponding

Christoffel numbers , j = 1, 2, ⋯ , N.

Multiplying (5) by Ll(t), l = 0, 1, ⋯ , N, and integrating it over the interval (0, T), one has

and hence

where

and [x] is the integer part of x.

Denote

Then (7) can be written as the following form

and can be evaluated by

The numerical scheme (9)-(10) was proposed by Wang and Guo [20] and called Legendre-Gauss-Radau collocation method.

2.2. LGRCM as a Runge-Kutta Method. We shall rewrite LGRCM (9)-(10) as a Runge-Kutta method. Taking f(t, u(t)) = tk−1 for any fixed k, 1 ≤ k ≤ N, we have

and

Define . Note that , u(t) ∈ PN(0, T), ∈ PN(0, T), , i = 1, 2, ⋯ , N. Thus, ≡ 0, which implies , t ∈ [0, T]. Integrating (11) with respect to t from t = 0 to ti and from t = 0 to T yields

Set

Substituting (13) and f(t, u(t)) = tk−1 into (9) and (10) gives

and

Similarly, we may obtain

It follows from (15) that

Since the nodes are distinct, is nonsingular. Multiplying (9) by , we have

Denote

Using (15) and (17), we can rewrite LGRCM (9)-(10) as

which is an N-stage Runge-Kutta method with step size T.

Remark 2.1. We make some remarks on LGRCM(9)-(10)

(a) To our knowledge, Runge-Kutta method (21) corresponding to LGRCM (9)-(10) does not belong to the classes of Radau IA or Radau IIA methods in the Runge-Kutta literature.

(b) From (19), (17) and , the coefficients of LGRCM(9)-(10) are independent of T.

(c) Comparing to Runge-Kutta method (21), the advantage of LGRCM(9)-(10) is that its coefficients are explicitly given in (8).

2.3. Convergence. From (19) and it follows that satisfy B(N), and , i, j = 1, 2,...,N, satisfy C(N), where

This implies that Runge-Kutta method (21) corresponding to LGRCM (9)-(10) is a collocation method and thus it is of order at least N. In order to obtain its exact order, we need the following result.

Lemma 2.1 ([13]). Let be real and distinct and be determined by Condition B(N). Then Condition B(2N −k) holds if and only if M(t) = is orthogonal to all polynomials of degree less than or equal to N − k − 1.

Theorem 2.2. The Runge-Kutta method (21) corresponding to LGRCM (9)-(10) is of order N.

Proof. Note that the collocation methods with different are of order at least N. Since the collocation points are zeros of . According to the property that the set of Ll(t) is a complete L2(0, T) -orthogonal system, can not be expanded by the polynomials LN,LN−1, ⋯ ,LN−k with k ≥ N. It follows from Lemma 2.1 that Condition B(k) only holds for k ≤ N, which implies that Runge-Kutta method (21) corresponding to LGRCM (9)-(10) has order N. □

2.4. Stability properties. We discuss the A-stability and algebraic-stability of Runge-Kutta method (21) corresponding to LGRCM (9)-(10).

The N-stage Runge-Kutta method (21) applied to u′ = λu yields = R(λT)u0 with R(z) = , where IN is the N-by-N identity matrix. R(z) is called the stability function of Runge-Kutta method (21). The stability function also satisfies

Runge-Kutta method (21) is called A-stable if |R(z)| ≤ 1 for all ℜ(z) ≤ 0. We express the stability function as the ratio of two polynomials

and define the E-polynomial as E(y) = D(iy)D(−iy) − N(iy)N(−iy).

Lemma 2.3 ([6]). A Runge-Kutta method with stability function R(z) = N(z)=D(z) is A-stable if and only if

(a) all poles of R(z) (that is, all zeros of D) are in the positive half plane, and

(b) E(y) ≥ 0 for all real y.

Theorem 2.4. Runge-Kutta method (21) corresponding to LGRCM (9)-(10) is A-stable for N = 1 and N = 2.

Proof. For N = 1, N(z) = and D(z) = . The zero of D(z) is . E(y) = for any real y. It follows from Lemma 2.3 that the 1-stage Runge-Kutta method (21) is A-stable. □

For N = 2, . The zeros of D(z) are z1 = 2 + and z2 = 2 - which lie in the positive half plane. Meanwhile, E(y) = for any real y. It follows from Lemma 2.3 that the 2-stage Runge-Kutta method (21) is A-stable.

It is difficult to get the expression of R(z) and check whether R(z) has no poles in the positive half plane by direct calculation of det and det for N ≥ 3. To overcome these difficulties, we use the explicit expression of R(z) for the collocation methods given by Nørsett [16] and follow the formulation of the Routh-Hurwitz algorithm in Wright [21].

Lemma 2.5 ([16, 21]). The stability function of the collocation method based on the points is given by

where M(t) = and CN is a constant.

According to the construction,

Lemma 2.6 ([21]). The zeros of the polynomial lie in the negative half plane if and only if all the coefficients have the same sign, where

In order to apply the RH algorithm [21] and to reorganize the coefficients, we apply the map . As a result, an alternative form of the stability function that we shall use is

and the condition (a) in Lemma 2.3 is transformed into that R⋆(z) has no poles in the positive half plane.

Theorem 2.7. Runge-Kutta method (21) corresponding to LGRCM (9)-(10) is not A-stable for N = 3, 4, 5.

Proof. For N = 3, we obtain by Lemma 2.5

Following the RH algorithm and using the expression of the coefficients of the polynomial denominator of R⋆(z), we obtain . All these terms have positive sign. Hence, R⋆(z) has no poles in positive half plane. After some calculation, we obtain E(y) = −420y4+60y6 which is negative for some real y. It follows from Lemma 2.3 that the 3-stage Runge-Kutta method (21) corresponding to LGRCM (9)-(10) is not A-stable.

For N = 4, we obtain by Lemma 2.5 the stability function

and

Following the RH algorithm and using the expression of the coefficients of the polynomial denominator of R⋆(z), we may verify that R⋆(z) has no poles in positive half plane. However, after some calculation, we obtain E(y) = −2688y6 + 96y8 which is negative for some real y. It follows from Lemma 2.3 that the 4-stage Runge-Kutta method (21) corresponding to LGRCM (9)-(10) is not A-stable.

For N = 5, we obtain by Lemma 2.5 the stability function

and

Following the RH algorithm and using the expression of the coefficients of the polynomial denominator of R⋆(z), we can verify that R⋆(z) has no poles in positive half plane. However, after some calculation, we obtain E(y) = 140y10−10080y8 +73920y6 which is negative for some real y. It follows from Lemma 2.3 that the 5-stage Runge-Kutta method (21) corresponding to LGRCM (9)-(10) is not A-stable. □

Conjecture. From the above result, we may conjecture that Runge-Kutta method (21) corresponding to LGRCM (9)-(10) is not A-stable for N ≥ 3.

We introduce the concept of algebraic stability of Runge-Kutta methods and quote a sufficient and necessary condition for the algebraic-stability of collocation methods. A Runge-Kutta method is algebraically-stable if , for i = 1, 2,...,N, and if the matrix M, given by

is positive semi-definite.

Lemma 2.8 ([6]). An N-stage algebraically stable collocation Runge-Kutta method must be of order at least 2N − 1.

Theorem 2.9. The Runge-Kutta method (21) corresponding to LGRCM (9)-(10) is algebraically-stable for N = 1 and not algebraically stable for N > 1.

Proof. The algebraic-stability for N = 1 follows from Definition 2.4, and the algebraic-stability for N > 1 follows from Lemma 2.8 and Theorem 2.2. □

 

3. Equivalence of LGCM and Gauss Method

In this section, we verify the equivalence of LGCM and Gauss Method.

3.1. Brief review. Different from the LGRCM, Guo and Wang [11] proposed to expand the collocation solution and the collocation points denoted by are taken as the zeros of the shifted Legendre orthogonal polynomial LN(t) on [0, T].

Denote

where

The numerical scheme LGCM proposed by Guo and Wang [11] is

3.2. Equivalence. We now rewrite the LGCM (22)-(23) as a Runge-Kutta method. Denote . Analogous to the previous analysis in Section 2, we obtain

Denote

The LGCM (22)-(23) can be rewritten as an N-stage Runge-Kutta method with step size T

In order to assert that LGCM (22)-(23) is equivalent to the Gauss method, we need the following result to prove its corresponding Runge-Kutta method (26) is of order 2N.

Lemma 3.1 ([6]). An N-stage Runge-Kutta method (26) has order 2N if and only if its coefficients satisfy

(a) the nodes are the zeros of LN(t),

(b) satisfy B(N),

(c) , i, j = 1, 2,...,N, satisfy C(N).

Theorem 3.2. LGCM (22)-(23) is equivalent to the Gauss method.

Proof. According to the construction of LGCM(22)-(23), the nodes are the zeros of LN(t) and thus (a) holds. From (24) and the invertibility of , C(N) is satisfied and thus (c) holds. Finally, it follows from (25) that B(N) is satisfied and hence (b) holds as well. It follows from Lemma 3.1 that LGCM(22)-(23) is equivalent to the Gauss method. □

We may obtain the following corollary from the equivalence result.

Corollary 3.3. LGCM (22)-(23) is A-stable and algebraically-stable.

 

4. Numerical experiments

We conduct some numerical experiments to illustrate the accuracy and stability of the LGRCM, and compare the LGCM with the Gauss method.

Example 4.1. Consider the following initial value problem [20]

The exact solution is u(t) = + sin 10πt.

In Figure 1, we have plotted the values of Δ of LGRCM with various N and step size T, where Δ = denotes the number of correct digits of the numerical solution at the end point t = 1. We observe a close agreement between theoretical and numerical convergence order.

Figure 1.Global errors of LGRCM with N = 1, 2, 3.

Example 4.2. The A-stability of LGRCM is tested by the following initial value problem

The numerical solutions generated by the LGRCM with step size for different N are plotted in Figure 2. The numerical solution produced by LGRCM with N = 2 tends to zero. However, the numerical solutions produced by LGRCM with N = 3, 4 do not tend to zero, which indicates they are not A-stable.

Figure 2.Numerical solutions by LGRCM with N = 2, 3, 4.

Example 4.3. Consider the following initial value problem [11]

The exact solution u(t) = + 5 sin 2t oscillates and grows to infinity as t → ∞.

The one-step errors of the LGCM and the Gauss method with various N, T are plotted in Figure 4. We observe that both methods perform extremely well and similar despite the influence of roundoff errors, which confirms our theoretical result.

Figure 3.Global errors of LGCM and Gauss method with various N.