A unified numerical method for computing the quadrature weights, integration matrix, and differentiation matrix is newly developed in this study. For this purpose, a splinelike 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 Gausstype formulas are frequently used in solving the differential and integrodifferential 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 wellbehaved 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 splinelike interpolation function can be converted into a global spline interpolating (GSPIN) formula that has exactly the same form as in the global Lagrange interpolation. CQF computations are performed using integration and differentiation of the GSPIN 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 GSPIN formula. Therefore, Gausstype 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 LegendreGauss, LegendreGauss Lobatto, and ChevyshevGaussLobatto 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 Δ
t_{m}
=
t
_{m+1}
−
t_{m}
.
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
a_{m,n}
, in (1) can be represented by a linear combination of function values
as
Using the assumption, the function can be approximated over
t
∈ [
t_{m}
,
t
_{m+1}
] with the following GSPIN formula.
where
The partial integration and the p
^{th}
derivative of the function at any points over
t
∈ [
t_{m}
,
t
_{m+1}
] can be approximated using (6)  (8), which come from their exact derivation using (4) and (5).
where
The partial integration over
t
∈ [
t
_{0}
,
t_{j}
] and the differentiation at
t_{j}
can be represented by
Therefore, the CQFs can be represented for the traditional spectral method, where the domain [
t
_{0}
,
t_{M}
] 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 GSPIN formula. For this, the following reformulation is applied for
t
∈ [
t_{m}
,
t
_{m+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 GSPIN form as shown in (4) and (5). The corresponding CQFs can be computed using (12) ~ (17).
This procedure is represented for the Lobattotype 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 Radautype 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
∈ [
t
_{0}
,
t
_{1}
] can be obtained using the procedures shown in (i)(iii) for each of the Lagrange interpolation polynomials defined with N Radautype points
. Therefore, the proposed method can be used without limits for other types of nodes, such as Gausstype and flipped Radautype 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 GSPIN 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 p
^{th}
functionderivative vector corresponding to the collocation points.
The GSPIN 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
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) LegendreGauss (LG) points
: the roots of
g
(
t
) =
L
_{N+1}
(
t
) = 0.
(b) LegendreGaussLobatto (LGL) points
: the roots of
g
(
t
) =
L
_{N+1}
(
t
) −
L
_{N−1}
(
t
) = 0 .
(c) ChebyshevGaussLobatto (CGL) points
:
where
In these equations the orthogonal polynomials
L_{j}
(
t
) and
T_{k}
(
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 ChristoffelDarboux identity theorem for the kernel function defined with orthogonal polynomials
[9]
. The integration matrix for (c) represents the ClenshawCurtis algorithm for Chebyshev polynomials
[10
,
13
]. The accuracy of (12)~(17) is estimated from the rootmeansquares (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.
Normalized differences in the CQFs
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
t
_{0}
= −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 GSPIN 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
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].
Test functions
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).
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 Δ
I_{j}
 , Δ
D_{j}
 , 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
).
Computational errors for f_{1} (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.
Computational errors for f_{2} (t) and f_{3} (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 GSPIN form. As shown in the authors’ paper
[12]
and from the present study, the spline interpolation is a natural choice for the GSPIN formula. For completeness, we repeat the main results of the earlier paper in this Section with detailed procedures for the derivation of the GSPIN 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
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 interconnecting 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 notaknot spline. The natural spline can be obtained by setting higher order derivatives at two end points to zeroes as
where
N
_{0}
and
N_{M}
are positive integers and can be selected under the following condition.
N_{0} + N_{M} = N − 1
In contrast, the notaknot spline imposes continuity conditions for the
N
^{th}
order derivatives near two end points with the same condition in selecting
N
_{0}
and
N_{M}
, such as
or
a_{j,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
f_{j}
=
δ _{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 GSPIN formula. The tension spline interpolation for a data set
related to a function
f
(
t
) ∈
C
^{2}
[
t
_{0}
,
t_{f}
] can be defined using the solution of the fourth order differential equation
[16]
for each time interval
t
∈ [
t_{m}
,
t
_{m+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 GSPIN formula.
Using the following definition,
(58) can be rewritten as
Therefore, the first and third derivatives of the function over
t
∈ [
t_{m}
,
t
_{m+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 nondimensional tension parameter
ρ_{m}
.
Table 2
shows the estimated magnitude of the truncated terms varying the nondimensional 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 GSPIN 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
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
a_{j}
(
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 thirdorder 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
d
f
(
ξ
(
τ_{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
d
f
(
ξ
(
τ_{j}
)) /
dτ
with reference to
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 fourpoint 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 nondimensional 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.
Computational errors for f_{1}(t) and f_{2}(t)
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.
Computational errors for f_{1}(t) and f_{2}(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 multigrid 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 pseudospectral method
[15]
and can be further extended to the optimal estimation problems with nonstochastic 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 LegendreGauss, Legendre GaussLobatto, and ChevyshevGaussLobatto 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 GaussLobatto 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 multigrid 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
ChangJoo 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.
JoonGoo 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.
Kim H.H.
,
Lee C.J.
,
Park J.B.
,
Shin J.R.
,
Jeong S.Y.
2008
‘TwoTerminal Numerical Algorithm for SinglePhase 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 Magnetohydrodynamic 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 BSpline Parametrization,’
Journal of Electrical Engineering & Technology
6
(4)
513 
518
DOI : 10.5370/JEET.2011.6.4.513
Benson D.
2004
“A Gauss Pseudospectral Transcription for Optimal Control,” Ph.D. Thesis
MIT, Department of Aeronautics and Astronautics
Huntington G. T.
2007
“Advancement and Analysis of a Gauss Pseudospectral 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 SemiInfinite Domains via Collocation Methods,”
Journal of Applied Mathematics
Article ID 696574
2012
21 
Williams P.
2005
“A GaussLobatto Quadrature Method for Solving Optimal Control Problems,”
Australian and New Zealand Industrial and Applied Mathematics Journal
47
101 
115
Elgendi S. E.
1969
“Chebyshev Solution of Differential, Integral and IntegroDifferential 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 913
Waldvogel J.
2006
“Fast Construction of the Fejér and ClenshawCurtis Quadrature Rules,”
BIT Numerical Mathematics
46
(1)
195 
202
DOI : 10.1007/s1054300600454
Barden J. M.
2013
“A Modified ClenshawCurtis 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