A New Unified Scheme Computing the Quadrature Weights, Integration and Differentiation Matrix for the Spectral Method

Abstract

A unified numerical method for computing the quadrature weights, integration matrix, and differentiation matrix is newly developed in this study. For this purpose, a spline-like interpolation using piecewise continuous polynomials is converted into a global spline interpolation formula, with which the quadrature formulas can be derived from integration and differentiation of the transformed function in an exact manner. To prove the usefulness of the suggested approach, both the Lagrange and tension spline interpolations are represented in exactly the same form as global spline interpolation. The applicability of the proposed method on arbitrary nodes is illustrated using two different sets of nodes. A series of validations using three test functions is conducted to show the flexibility in selecting computational nodes with the present method.

1. Introduction

Numerical methods for system analysis and design optimization are widely studied in many engineering problem, where much effort is paid on selecting node or boundary conditions in each applications [1-3]. In consideration of both efficiency and accuracy, this paper proposes a new unified numerical scheme in the computation of the quadrature weights, integration matrix, and differentiation matrix for spectral methods, collectively referred to in this paper as the coefficients for quadrature formulas (CQFs). The Gauss-type formulas are frequently used in solving the differential and integro-differential equations arising in many engineering and scientific areas, since they can provide extreme accuracy and high convergence at an exponential rate to iterative analyses for problems with well-behaved solutions [4-7]. CQFs for spectral methods are generally derived using the orthogonal properties among interpolating polynomials. Many of them, expressed in an analytic form, are easily accessible in the literature [8-10]. Derivation of CQFs generally requires rigorous application of complex function theories and is strictly dependent on the types of interpolating functions and quadrature nodes, which may cause inconvenience when they are required for other interpolating functions and quadrature nodes different from the standard ones.

This work proposes a new numerical approach in the computation of CQFs that requires only information on the computational nodes and can retain the same level of accuracy as traditional methods. For this purpose, the spline interpolation approach utilizing the piecewise continuous polynomials is used, rather than the global Lagrange interpolation. The spline-like interpolation function can be converted into a global spline interpolating (G-SPIN) formula that has exactly the same form as in the global Lagrange interpolation. CQF computations are performed using integration and differentiation of the G-SPIN formula in an exact manner. As a result, computed CQFs can attain high precision equivalent to those obtained using the traditional spectral method. This paper shows that Lagrange interpolation can be represented by a G-SPIN formula. Therefore, Gauss-type formulas for any kind of computational node can be computed using the proposed methods. The present method is validated through a series of applications using Legendre-Gauss, Legendre-Gauss- Lobatto, and Chevyshev-Gauss-Lobatto points and by showing that a high accuracy in predicting CQFs is obtained and equivalent to those reached using the traditional quadrature formulas.

The extended applicability of the proposed method to solutions on arbitrary nodes is shown next. For this purpose, two different sets of nodes are considered. One is generated using the coordinate transformation devised by the authors [11] for more uniformly distributed nodes than those obtained using the traditional spectral method. The other is a set of uniform nodes known as the worst possible choice for polynomial interpolation, with which the derived formulas are extremely inaccurate as mentioned in [4] and [12]. The Lagrange and tension spline interpolations over these sets of nodes are investigated by applying the present method to three carefully selected test functions and by comparing the resultant absolute errors in predicting integration and differentiation for those functions. The results of comparative studies demonstrate that the quadrature weights can be retained positivity with the tension spline interpolation regardless on the number of nodes. Whereas, the Lagrange interpolation show highly oscillatory behaviors in the quadrature weights, some of which become negative, as the number of nodes are increased. Therefore, the tension spline is more robust for uniform nodes than the Lagrange interpolation. In addition, it is indicated that the applications of advanced grid topologies, such as the local adaptive grid and multigrid techniques, is possible to enhance both accuracy and convergence for the iterative solutions based on the quadrature formula with the flexibility allowed by the preset method in selecting computational nodes.

 

2. Unified Quadrature Formulas using Spline Interpolation

The spline interpolation for a set of function data given over the computational nodes can be represented using the (M +1)× (N +1) unknown coefficients as

where N is the leading order of the piecewise continuous interpolating polynomials and τm ∈[0,1] is a nondimensional independent variable defined by

where Δtm = tm+1 − tm .

For deriving CQFs corresponding to the set of nodes , this paper makes the following assumption, the rationale of which will become clear in the later part of this section.

Assumption: The unknown coefficient am,n , in (1) can be represented by a linear combination of function values as

Using the assumption, the function can be approximated over t ∈ [tm,tm+1] with the following G-SPIN formula.

where

The partial integration and the pth derivative of the function at any points over t ∈ [tm,tm+1] can be approximated using (6) - (8), which come from their exact derivation using (4) and (5).

where

The partial integration over t ∈ [t0,tj] and the differentiation at tj can be represented by

Therefore, the CQFs can be represented for the traditional spectral method, where the domain [t0,tM] of interests is transformed into [.1, 1] using the affine transformation [4-7].

where

Therefore, once the coefficients are known for the set of arbitrary nodes , the corresponding CQFs can be computed using (12) ~ (17).

To access the usefulness of the assumption and the resultant formulas in obtaining CQFs, we first formulate the coefficients {αm,n,k} corresponding to the Lagrange interpolation as follows. The Lagrange interpolation uses continuous polynomials with global support, but it can be reformulated using piecewise continuous polynomials of the same form as the G-SPIN formula. For this, the following reformulation is applied for t ∈ [tm,tm+1] using the definition given in (2).

The required coefficients can be obtained using the following recursive procedures.

Therefore, the Lagrange interpolation for a data set defined over arbitrarily distributed nodes can be reduced to G-SPIN form as shown in (4) and (5). The corresponding CQFs can be computed using (12) ~ (17).

This procedure is represented for the Lobatto-type quadrature points containing two endpoints in the domain of interest. However, minor modifications are sufficient for its extended application to other types of nodes. As an example with the Radau-type quadrature points, the given data set can be defined by with the order of interpolating polynomials M = N −1. Even in this case, the coefficients corresponding to the first interval t ∈ [t0,t1] can be obtained using the procedures shown in (i)-(iii) for each of the Lagrange interpolation polynomials defined with N Radau-type points . Therefore, the proposed method can be used without limits for other types of nodes, such as Gauss-type and flipped Radau-type points, and regardless of the adopted orthogonal polynomials, such as Legendre and Chebyshev polynomials.

The multiple derivatives of the function can easily be formulated using the G-SPIN formula shown in (4) and (5). The differentiation matrices for obtaining higher function derivatives are typically obtained by recursively applying the differentiation formula given in (14). Using the following definition for vectors and matrices,

the first derivative can be expressed by

Recursive application of (21) results in the conventional formula for the approximation of the pth function-derivative vector corresponding to the collocation points.

The G-SPIN formula allows multiple differentiation to get the pth derivatives of the function at each node, as shown in (10) and (11). If the components of the matrix are denoted by , then each component can be obtained directly using these equations as

respectively.

 

3. Validation of the Proposed Method

The proposed formulas for computing CQFs are validated using existing formulas available for the pseudospectral methods. For this purpose, we consider three sets of quadrature points defined by (a)-(c) with the corresponding formulas for the computation of CQFs.

(a) Legendre-Gauss (LG) points : the roots of g(t) = LN+1(t) = 0.

(b) Legendre-Gauss-Lobatto (LGL) points : the roots of g(t) = LN+1(t) − LN−1(t) = 0 .

(c) Chebyshev-Gauss-Lobatto (CGL) points :

where

In these equations the orthogonal polynomials Lj (t) and Tk (t) represent the jth Legendre and kth Chebyshev polynomials, respectively. The integration matrices for (a) and (b) come from Axelson’s algorithm derived using the Christoffel-Darboux identity theorem for the kernel function defined with orthogonal polynomials [9]. The integration matrix for (c) represents the Clenshaw-Curtis algorithm for Chebyshev polynomials [10, 13]. The accuracy of (12)~(17) is estimated from the root-meansquares (RMS) value of the differences in the components of the CQFs, defined by

where is computed using (10). The normalized norm defined by ||ΔD||n = ||ΔD||RMS/||D||RMS is used for convenience in comparison.

Table 1 shows that the computed results and the proposed method can predict CQFs regardless of the number of nodes at a level of accuracy equivalent to that of spectral methods.

Table 1.Normalized differences in the CQFs

The computational domain typically includes both end points to impose the initial and final conditions. As an example of the differentiation method with LG points, the modified Lagrange interpolation functions are used openly to include the effect of the initial end point. The Lagrange interpolation functions for the LG points and the approximation of the function and its first derivatives at each node can be represented by (35) and (36).

The corresponding modified Lagrange interpolation functions, the first of which corresponds to the initial end point t0 = −1 , are typically defined as

The function and its first derivative at each node can be approximated with as

where , ~ can be defined using (36) and (37) for j = 0,⋯, N and k =1,⋯,M −1 as

The elements corresponding to k = 0 can be computed using the fact that the derivative of a constant function is uniformly zero at all nodes.

These derivations are typically used in the spectral method. The modified differentiation matrix is obtainable directly from the G-SPIN formula by including the data (t0, f0 ) in the data set defined over the LG points. Likewise, the elements related to the final endpoint are easily computed in the same manner, if required. Consequently, the proposed method allows straightforward computing of all elements of the differentiation matrix related to all nodes including two end points, regardless of the types of nodes and without computing the modified differentiation matrix using formulas such as (40) and (41).

The accuracy of the proposed method is validated using the following three test functions selected from [14]. Fig. 1 shows these functions over t ∈[−1,1].

Fig. 1.Test functions

The error function erf(y) in f2 (t) is an entire function defined by (42) and can be approximated using the series expansion of the exponential function as shown in (43). Therefore, the first derivative of f2 (t) can be approximated by (45).

The accuracy of (12)~(17) is estimated using the RMS norms of the absolute error vectors for partial integration and differentiation and by the absolute value of the integration error over t ∈ [-1, 1], presented with ||ΔIj|| , ||ΔDj|| , and ||ΔI|| , respectively. The integrations over the entire domain for f1(t) and f2(t) are given, respectively, by I2 (-1, 1) ≈ 0.175664900305971 and I3 (-1, 1)=- 0.176358246030565 [14]. These values are used to measure the absolute integration errors. Fig. 2 shows the computational errors using different formulas and with varying numbers of nodes for f1 (t).

Fig. 2.Computational errors for f1 (t).

Fig. 3 shows the results for f2 (t) and f3 (t). The present method shows nearly the same level of accuracy as the conventional methods in both the integration and the differentiation varying the number of nodes. The present method is more advantageous in obtaining the highly accurate first derivative for f1 (t) and f3 (t) in the case when the number of nodes is sufficiently large to provide machine precision. The LGL points are the better selection for accurate integration than the CGL points for all test functions.

Fig. 3.Computational errors for f2 (t) and f3 (t)

 

4. Extension to Spline Interpolation

Since the formulas given in (12)-(17) are derived on the basis of the assumption shown in (3), other interpolations rather than the Lagrange interpolation can be used, once they are represented in the G-SPIN form. As shown in the authors’ paper [12] and from the present study, the spline interpolation is a natural choice for the G-SPIN formula. For completeness, we repeat the main results of the earlier paper in this Section with detailed procedures for the derivation of the G-SPIN formula for spline interpolation. From the definition given by (1), the function values and its derivatives up to the kth order can be approximated at each node by

The unknown coefficients for spline interpolation are typically determined by imposing the following conditions over the computational nodes.

i) function values at each node

ii) continuity of the function and its higher derivatives at each inter-connecting node

iii) (N-1) additional conditions to close the solution

There exist many alternatives for defining the additional conditions (iii). The most widely used among them are the natural spline and the not-a-knot spline. The natural spline can be obtained by setting higher order derivatives at two end points to zeroes as

where N0 and NM are positive integers and can be selected under the following condition.

N0 + NM = N − 1

In contrast, the not-a-knot spline imposes continuity conditions for the N th order derivatives near two end points with the same condition in selecting N0 and NM, such as

or aj,N = 0.

From the conditions given in (48)-(53), the unknown coefficients can be obtained by solving the following linear algebraic equation.

where a and f are defined by

The coefficient matrix C has the dimension M(N +1) - by- (M +1) and most of its components are zeroes, except at the rows related to condition (i) when the natural or nota- knot spline is used. The matrix M is easily defined by the coefficients in (48) ~ (53).

Therefore, the coefficients correspondding to the function value fk can be computed from the solution of (54) by substituting the elements of the function vector with fj = δ jk( j = 0, ⋯ M) . Therefore, the CQFs corresponding to the spline interpolation can be derived using (12) ~ (17) regardless of what kinds of (N-1) additional conditions (iii) are imposed.

Next, we consider the application of the proposed method using the tension spline interpolation. The tension spline is not a kind of the piecewise continuous polynomial interpolation. However, it will be shown that the tension spline interpolation can be represented by the G-SPIN formula. The tension spline interpolation for a data set related to a function f(t) ∈C2[t0, tf] can be defined using the solution of the fourth order differential equation [16] for each time interval t ∈ [tm,tm+1] with an adjustable tension parameter as given by

with the following boundary conditions

and using the definition given in (3), the solution of (55) can be represented by

The unknown second derivatives are typically determined by imposing the continuity conditions for the first derivatives at each interconnecting node and by adding two end conditions. The possible selection of two end conditions is exemplified in [13] with the following three types.

Types I and III need the direct specification of the first and second derivatives at two ends, which is inconvenient for general cases without prior information on these values. In addition, Type III is not the usual one for the tension spline as mentioned in [16], and the following modified end condition using the third derivatives is proposed for the application of the present G-SPIN formula.

Using the following definition,

(58) can be rewritten as

Therefore, the first and third derivatives of the function over t ∈ [tm,tm+1] can be approximated by

The continuity conditions for the first derivatives can be expressed as

Two end conditions shown in (58) can be represented by

As the results of (60) and (61), the unknown second derivatives can be obtained by solving the following system of the linear algebraic equations in the same form as given by (54).

The trigonometric hyperbolic functions in (59) can be expressed by the linear combination of eρmτm and e−ρmτm using the following relations.

In addition, the hyperbolic functions e±ρmτm can be approximated using the Taylor series expansion as

where O(⋅) is the magnitude of the truncated terms. An accurate approximation for each of the exponential functions can be obtained by increasing N until O(⋅) reaches machine precision. Since |τm| ≤1 , the magnitude of determines the required N, depending on the non-dimensional tension parameter ρm . Table 2 shows the estimated magnitude of the truncated terms varying the non-dimensional tension parameter. With the tension parameter less than ρ = 0.01, the tension spline approaches the cubic spline (with N=3). For the extreme case with ρ =10.0, more than 50 terms are required to get machine precision accuracy. However, this does not cause a large increase in computing time, because the arithmetic operations for (12)-(17) are simple enough, once is known, and the CQFs can be computed prior to their applications. Therefore, (59) can be reduced to G-SPIN form with high precision, and (12) ~ (17) can be used in an efficient manner even with tension spline interpolation.

Table 2.Truncation errors with varying tension parameter

The traditional quadrature nodes are highly clustered around two end points to preserve extreme accuracy, and there is no flexibility in adjusting their spacing. This kind of distribution can be inconvenient in some problems with rapid variations in function values around the midpoint of the computational domain, as shown in [11]. In addition, conventional quadrature weights and matrices for the integration and differentiation might be inadequate for good convergence in the iterative solution processes as the number of nodes is increased [11, 15]. Spline interpolation combined with the present approaches to the CQFs can be a good alternative allowing an adaptive selection of nodes considering the solution properties. This possibility will be demonstrated using two different kinds of node. First, more evenly distributed nodes generated using the coordinate transformation method shown in [8] are used. Second, uniform nodes are selected to show that arbitrary nodes can be utilized in real applications without limits on the number of nodes and without numerical failures using the spline interpolation. With the given set of conventional quadrature nodes , the coordinate transformation proposed in [11] is devised to preserve symmetry about τ = 0 and monotonicity in the corresponding computational nodes , using the following equation.

The coefficients aj (j =0,1,2, …) have been determined by the least squares method to match best with target uniform nodes . [11] showed the effect of the leading order, K, with the selection of the third-order transformation being recommended, since it provides the best integration accuracy for the series of polynomial test functions. The formulas for integration and differentiation can be defined easily in the transformed domain using (64) as

If accurate interpolations over an entire domain of interest are possible, the derivative df (ξ (τj)) / dτ can be estimated with high precision in principle. In such a case, the interpolating functions and can be defined to approximate the function and its derivatives using the nodes , respectively. The following formula can be derived from a simple manipulation to estimate df (ξ (τj)) / dτ with reference to Fig. 4.

Fig. 4.Schematic diagram to estimatedf (ξ (τj)) / dτ using interpolation

where

However, it is not a simple task to find such interpolating functions with the required accuracy in real applications. Therefore, we use the local Lagrange interpolating functions using four-point data in this study to demonstrate the performance of (66), because the global one generally results in extremely poor accuracy, whereas the CQFs corresponding to the transformed nodes are directly obtainable using (12)~(17), which require only the information on node distribution. The Lagrange interpolation and the tension spline with the non-dimensional tension parameter of ρ = 0.01 are used for the applications of the present method. Fig. 5 compares the results of different approaches for f1(t) and f2(t). The integration and differentiation using the coordinate transformation are estimated using (65) and (66). The coordinate transformation method shows nearly the same performance as the spectral method in predicting the integration as the number of nodes is increased. However, the differentiation computed with (66) contains a larger error than that with the present method when the tension spline interpolation is used. The present results with Lagrange interpolation show rapid increase in computational errors for both the differentiation and the integration with more than 40 nodes. Fig. 6 presents the effect of the number of nodes on the variation of the quadrature weights.

Fig. 5.Computational errors for f1(t) and f2(t)

Fig. 6.Quadrature weight distribution with varying node number.

The oscillatory behaviors are made more apparent by increasing the number of nodes, and some weights even become negative owing to poor interpolation accuracy, whereas the present method using the tension spline shows improvement in accuracy both for the integration and for the differentiation as the number of nodes is increased. Even though the same spectral accuracy cannot be obtained as that achievable with the spectral method using the LGL nodes, this method can be used for the problems, where reasonable accuracy is enough for its applications.

Finally, the applicability of the proposed method on arbitrary nodes is investigated. Efforts to utilize arbitrary nodes include recent studies by Ross et al. [4] and Gong et al. [12] showing that convergence for the solution of nonlinear optimal control problems with pseudospectral methods can be guaranteed only when all quadrature weights are positive. They pointed out that analyses with more than ten uniform nodes might fail, because at least one of the weights becomes negative when Lagrange polynomials are used on uniform nodes. Regarding the positivity requirement for the weights, the spline interpolation is more robust than Lagrange interpolation, as shown in Figs. 5 and 6, because the corresponding weights retain positivity even up to 150 nodes. For these reasons, we use the tension splines in the present investigation. Fig. 7 shows the prediction accuracy with different tension parameters for f1(t) and f2(t). As expected, uniform nodes present much worse accuracy than the LGL points. However, absolute prediction errors using uniform nodes gradually decrease as the number of nodes is increased. In the cases when these errors are allowable with a reasonable number of nodes, uniform nodes can also be used, meaning that arbitrary nodes can be selected for real applications in the mentioned conditions.

Fig. 7.Computational errors for f1(t) and f2(t)

Two aspects are considered for a further work. First, various grid topologies can be used to enhance both accuracy and convergence in the iterative solutions using the flexible selection of nodes allowed by the proposed method. The local adaptive grid and multi-grid techniques developed for computational fluid dynamics can be promising topologies for those purposes.

Second, the proposed enhanced computation method can provide theoretic background in analyzing and resolving a group of optimal control and state estimation problems with the help of NOCP formulation. Application to optimal control problem with practical inequality constraints is already under investigation through pseudo-spectral method [15] and can be further extended to the optimal estimation problems with non-stochastic system models.

 

5. Conclusion

This paper proposed a unified approach to the computation of the coefficients for quadrature formulas, such as the quadrature weights, integration matrix, and differentiation matrix. These formulas were derived from the exact integration and differentiation of the piecewise continuous polynomials in the global spline interpolation formula. To demonstrate the usefulness of the present approach, it was shown that Lagrange and tension spline interpolation can be transformed into the global spline interpolation formula, regardless of the type of nodes. The present method was adopted to predict the coefficients of the quadrature formulas for the Legendre-Gauss, Legendre- Gauss-Lobatto, and Chevyshev-Gauss-Lobatto points with levels of precision equivalent to those attainable with traditional spectral methods. The proposed method was validated through the application of the present method to carefully selected test functions. The results are compared with those using the traditional quadrature formula. The integrations for all functions can be approximated in nearly the same level of accuracy with both methods. However, the differentiations can be more accurately predicted with the present method than with the traditional method.

As an effort to extend the applicability of the proposed method to solutions on arbitrary nodes, nodes generated using an analytic coordinate transformation method and uniform nodes were selected to show the accuracy achievable with the proposed method. The coordinate transformation method can be used to compute the quadrature integration with high precision equivalent to that obtained using the spectral method at the Legendre- Gauss-Lobatto nodes. However, it revealed poor accuracy in predicting the differentiation owing to high interpolation error when Local Lagrange interpolation was adopted.

For the proposed method, spline interpolation outperformed Lagrange interpolation in that its accuracy in predicting the integration and differentiation improved uniformly as the number of nodes increased for both sets of nodes. Therefore, the proposed method using spline interpolation can be used on arbitrary nodes when the incurred errors are allowable for real applications. The flexibility in selecting computational nodes allows the use of advanced grid topologies, such as the local adaptive grid and multi-grid techniques, to enhance both accuracy and convergence in iterative solution processes. As a result of this study, a unified method of computing quadrature weights, the integration matrix, and the differentiation matrix corresponding to the quadrature nodes for the spectral methods has been established. Furthermore, the proposed methods are applicable to arbitrarily distributed nodes.