Advanced
A New Unified Scheme Computing the Quadrature Weights, Integration and Differentiation Matrix for the Spectral Method
A New Unified Scheme Computing the Quadrature Weights, Integration and Differentiation Matrix for the Spectral Method
Journal of Electrical Engineering and Technology. 2015. May, 10(3): 1188-1200
Copyright © 2015, The Korean Institute of Electrical Engineers
This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/)which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.
  • Received : June 03, 2014
  • Accepted : February 09, 2015
  • Published : May 01, 2015
Download
PDF
e-PUB
PubReader
PPT
Export by style
Article
Author
Metrics
Cited by
TagCloud
About the Authors
Chang-Joo Kim
Dept. of Aerospace Information Engineering, Konkuk University, Korea. (cjkim@konkuk.ac.kr)
Joon-Goo Park
School of Electrical Engineering, Kyungpook National University Korea. (jgpark@ee.knu.ac.kr)
Sangkyung Sung
Corresponding Author: Dept. of Aerospace Information Engineering, Konkuk University, Korea. (sksung@konkuk.ac.kr)

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.
Keywords
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
PPT Slide
Lager Image
given over the computational nodes
PPT Slide
Lager Image
can be represented using the ( M +1)× ( N +1) unknown coefficients
PPT Slide
Lager Image
as
PPT Slide
Lager Image
where N is the leading order of the piecewise continuous interpolating polynomials and τm ∈[0,1] is a nondimensional independent variable defined by
PPT Slide
Lager Image
where Δ tm = t m+1 tm .
For deriving CQFs corresponding to the set of nodes
PPT Slide
Lager Image
, 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
PPT Slide
Lager Image
as
PPT Slide
Lager Image
Using the assumption, the function can be approximated over t ∈ [ tm , t m+1 ] with the following G-SPIN formula.
PPT Slide
Lager Image
where
PPT Slide
Lager Image
The partial integration and the p th derivative of the function at any points over t ∈ [ tm , t m+1 ] can be approximated using (6) - (8), which come from their exact derivation using (4) and (5).
PPT Slide
Lager Image
PPT Slide
Lager Image
where
PPT Slide
Lager Image
The partial integration over t ∈ [ t 0 , tj ] and the differentiation at tj can be represented by
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
Therefore, the CQFs can be represented for the traditional spectral method, where the domain [ t 0 , tM ] of interests is transformed into [.1, 1] using the affine transformation [4 - 7] .
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
where
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
Therefore, once the coefficients
PPT Slide
Lager Image
are known for the set of arbitrary nodes
PPT Slide
Lager Image
, 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 , t m+1 ] using the definition given in (2).
PPT Slide
Lager Image
The required coefficients can be obtained using the following recursive procedures.
PPT Slide
Lager Image
Therefore, the Lagrange interpolation for a data set
PPT Slide
Lager Image
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
PPT Slide
Lager Image
with the order of interpolating polynomials M = N −1. Even in this case, the coefficients
PPT Slide
Lager Image
corresponding to the first interval t ∈ [ t 0 , t 1 ] can be obtained using the procedures shown in (i)-(iii) for each of the Lagrange interpolation polynomials defined with N Radau-type points
PPT Slide
Lager Image
. 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,
PPT Slide
Lager Image
PPT Slide
Lager Image
the first derivative can be expressed by
PPT Slide
Lager Image
Recursive application of (21) results in the conventional formula for the approximation of the p th function-derivative vector corresponding to the collocation points.
PPT Slide
Lager Image
The G-SPIN formula allows multiple differentiation to get the p th derivatives of the function at each node, as shown in (10) and (11). If the components of the matrix
PPT Slide
Lager Image
are denoted by
PPT Slide
Lager Image
, then each component
PPT Slide
Lager Image
can be obtained directly using these equations as
PPT Slide
Lager Image
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
PPT Slide
Lager Image
: the roots of g ( t ) = L N+1 ( t ) = 0.
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
(b) Legendre-Gauss-Lobatto (LGL) points
PPT Slide
Lager Image
: the roots of g ( t ) = L N+1 ( t ) − L N−1 ( t ) = 0 .
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
(c) Chebyshev-Gauss-Lobatto (CGL) points
PPT Slide
Lager Image
:
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
where
PPT Slide
Lager Image
In these equations the orthogonal polynomials Lj ( t ) and Tk ( t ) represent the j th Legendre and k th 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
PPT Slide
Lager Image
where
PPT Slide
Lager Image
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.
Normalized differences in the CQFs
PPT Slide
Lager Image
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
PPT Slide
Lager Image
and the approximation of the function and its first derivatives at each node can be represented by (35) and (36).
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
The corresponding modified Lagrange interpolation functions, the first of which corresponds to the initial end point t 0 = −1 , are typically defined as
PPT Slide
Lager Image
The function and its first derivative at each node can be approximated with
PPT Slide
Lager Image
as
PPT Slide
Lager Image
PPT Slide
Lager Image
where
PPT Slide
Lager Image
, ~ can be defined using (36) and (37) for j = 0,⋯, N and k =1,⋯, M −1 as
PPT Slide
Lager Image
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.
PPT Slide
Lager Image
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 ( t 0 , f 0 ) 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
PPT Slide
Lager Image
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].
PPT Slide
Lager Image
Test functions
PPT Slide
Lager Image
The error function erf( y ) in f 2 ( 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 f 2 ( t ) can be approximated by (45).
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
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 f 1 ( t ) and f 2 ( t ) are given, respectively, by I 2 (-1, 1) ≈ 0.175664900305971 and I 3 (-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 f 1 ( t ).
PPT Slide
Lager Image
Computational errors for f1 (t).
Fig. 3 shows the results for f 2 ( t ) and f 3 ( 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 f 1 ( t ) and f 3 ( 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.
PPT Slide
Lager Image
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 k th order can be approximated at each node by
PPT Slide
Lager Image
PPT Slide
Lager Image
The unknown coefficients
PPT Slide
Lager Image
for spline interpolation are typically determined by imposing the following conditions over the computational nodes.
i) function values at each node
PPT Slide
Lager Image
PPT Slide
Lager Image
ii) continuity of the function and its higher derivatives at each inter-connecting node
PPT Slide
Lager Image
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
PPT Slide
Lager Image
PPT Slide
Lager Image
where N 0 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 N 0 and NM , such as
PPT Slide
Lager Image
or aj,N = 0.
From the conditions given in (48)-(53), the unknown coefficients
PPT Slide
Lager Image
can be obtained by solving the following linear algebraic equation.
PPT Slide
Lager Image
where a and f are defined by
PPT Slide
Lager Image
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
PPT Slide
Lager Image
in (48) ~ (53).
Therefore, the coefficients
PPT Slide
Lager Image
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
PPT Slide
Lager Image
related to a function f ( t ) ∈ C 2 [ t 0 , tf ] can be defined using the solution of the fourth order differential equation [16] for each time interval t ∈ [ tm , t m+1 ] with an adjustable tension parameter
PPT Slide
Lager Image
as given by
PPT Slide
Lager Image
with the following boundary conditions
PPT Slide
Lager Image
and using the definition given in (3), the solution of (55) can be represented by
PPT Slide
Lager Image
The unknown second derivatives
PPT Slide
Lager Image
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.
PPT Slide
Lager Image
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.
PPT Slide
Lager Image
Using the following definition,
PPT Slide
Lager Image
(58) can be rewritten as
PPT Slide
Lager Image
Therefore, the first and third derivatives of the function over t ∈ [ tm , t m+1 ] can be approximated by
PPT Slide
Lager Image
The continuity conditions for the first derivatives can be expressed as
PPT Slide
Lager Image
Two end conditions shown in (58) can be represented by
PPT Slide
Lager Image
As the results of (60) and (61), the unknown second derivatives
PPT Slide
Lager Image
can be obtained by solving the following system of the linear algebraic equations in the same form as given by (54).
PPT Slide
Lager Image
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.
PPT Slide
Lager Image
In addition, the hyperbolic functions e±ρmτm can be approximated using the Taylor series expansion as
PPT Slide
Lager Image
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
PPT Slide
Lager Image
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
PPT Slide
Lager Image
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.
Truncation errors with varying tension parameter
PPT Slide
Lager Image
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
PPT Slide
Lager Image
, the coordinate transformation proposed in [11] is devised to preserve symmetry about τ = 0 and monotonicity in the corresponding computational nodes
PPT Slide
Lager Image
, using the following equation.
PPT Slide
Lager Image
The coefficients aj ( j =0,1,2, …) have been determined by the least squares method to match best with target uniform nodes
PPT Slide
Lager Image
. [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
PPT Slide
Lager Image
PPT Slide
Lager Image
If accurate interpolations over an entire domain of interest are possible, the derivative d f ( ξ ( τj )) / can be estimated with high precision in principle. In such a case, the interpolating functions
PPT Slide
Lager Image
and
PPT Slide
Lager Image
can be defined to approximate the function and its derivatives using the nodes
PPT Slide
Lager Image
, respectively. The following formula can be derived from a simple manipulation to estimate d f ( ξ ( τj )) / with reference to Fig. 4 .
PPT Slide
Lager Image
Schematic diagram to estimatedf (ξ (τj)) / using interpolation
PPT Slide
Lager Image
where
PPT Slide
Lager Image
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 f 1 ( t ) and f 2 ( 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.
PPT Slide
Lager Image
Computational errors for f1(t) and f2(t)
PPT Slide
Lager Image
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 f 1 ( t ) and f 2 ( 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.
PPT Slide
Lager Image
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.
Acknowledgements
This work was support by the Faculty Research Fund of Konkuk University in 2013.
BIO
Chang-Joo Kim He is an associate Professor in the Department of Aerospace Information Engineering at Konkuk University, Korea. His research interests include nonlinear optimal control, helicopter flight mechanics and helicopter system design.
Joon-Goo Park He received his Ph.D. degree in Electronic Engineering from Seoul National University in 2001. He is presently a professor at the School of Electrical Engineering at Kyungpook National University. His research interests include GPS, mobile software navigation and algorithm analysis.
Sangkyung Sung He is an associate Professor in the Department of Aerospace Information Engineering at Konkuk University, Korea. His research interests include inertial sensors, aerospace system control, navigation and estimation.
References
Kim H.-H. , Lee C.-J. , Park J.-B. , Shin J.-R. , Jeong S.-Y. 2008 ‘Two-Terminal Numerical Algorithm for Single-Phase Arcing Fault Detection and Fault Location Estimation Based on the Spectral Information,’ Journal of Electrical Engineering & Technology 3 (4) 460 - 467    DOI : 10.5370/JEET.2008.3.4.460
Kien Le Chi 2013 ‘Numerical Calculations and Analyses in Diagonal Type Magneto-hydro-dynamic Generator,’ Journal of Electrical Engineering & Technology 8 (6) 1365 - 1370    DOI : 10.5370/JEET.2013.8.6.1365
Kim M.-H. , Lee H.-B. , Kim H.-S. , Byun J.-K. 2011 ‘A New Unified Design Environment for Optimization of Electric Machines Based on Continuum Sensitivity and B-Spline Parametrization,’ Journal of Electrical Engineering & Technology 6 (4) 513 - 518    DOI : 10.5370/JEET.2011.6.4.513
Ross I. M. , Karpenko M. 2012 “A Review of Pseudospectral Optimal Control: From theory to flight,” Annual Reviews in Control 36 (2) 182 - 197    DOI : 10.1016/j.arcontrol.2012.09.002
Benson D. 2004 “A Gauss Pseudo-spectral Transcription for Optimal Control,” Ph.D. Thesis MIT, Department of Aeronautics and Astronautics
Huntington G. T. 2007 “Advancement and Analysis of a Gauss Pseudo-spectral Transcription for Optimal Control Problems,” Ph. D. Thesis MIT
Grandclèment P. 2006 “Introduction to Spectral Methods,” European Astronomical Society Publications Series 21 153 - 180
Maleki M. , Hashim I. , Abbasbandy S. 2012 “Analysis of IVPs and BVPs on Semi-Infinite Domains via Collocation Methods,” Journal of Applied Mathematics Article ID 696574 2012 21 -
Williams P. 2005 “A Gauss-Lobatto Quadrature Method for Solving Optimal Control Problems,” Australian and New Zealand Industrial and Applied Mathematics Journal 47 101 - 115
El-gendi S. E. 1969 “Chebyshev Solution of Differential, Integral and Integro-Differential Equations,” The Computer Journal 12 (3) 282 - 287    DOI : 10.1093/comjnl/12.3.282
Kim C.-J. , Sung S. K. , Park S. H. , Jung S. N. , Park T. S. 2014 “Numerical Time Scale Separation for Nonlinear Optimal Control Analyses with Applications to Rotorcrafts,” AIAA Journal of Guidance, Control and Dynamics 37 (2) 658 - 673    DOI : 10.2514/1.59557
Gong Q. , Ross I. M. , Fahroo F. 2009 “Pseudospectral Optimal Control on Arbitrary Grids,” Proc. of AAS / AIAA Astrodynamics Specialist Conference Pittsburgh, PA August 9-13
Waldvogel J. 2006 “Fast Construction of the Fejér and Clenshaw-Curtis Quadrature Rules,” BIT Numerical Mathematics 46 (1) 195 - 202    DOI : 10.1007/s10543-006-0045-4
Barden J. M. 2013 “A Modified Clenshaw-Curtis Quadrature Algorithm,” Master Thesis Worcester Polytechnic Institute, Department of Science In Applied Mathematics
Kim C. -J. , Sung S. 2015 “Efficient Spline Transcription Techniques for Nonlinear Optimal Control Analyses Using a Pseudospectral Framework,” IEEE Trans. Control Systems Technology in press
Pruess S. 1976 “Properties of Splines in Tension,” Journal of Approximation Theory 17 86 - 96    DOI : 10.1016/0021-9045(76)90113-1