tools/gkdata.hpp
00001 /* gkdata.hpp
00002  *
00003  * Class for computing Gauss-Kronrod data in QUADPACK storage scheme.
00004  * Refactored from the GaussKronrod class used by quadpack++.
00005  *
00006  * Copyright (C) 2010, 2011 Jerry Gagelman
00007  * 
00008  * This program is free software: you can redistribute it and/or modify
00009  * it under the terms of the GNU General Public License as published by
00010  * the Free Software Foundation, either version 3 of the License, or
00011  * (at your option) any later version.
00012  *
00013  * This program is distributed in the hope that it will be useful,
00014  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00015  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00016  * GNU General Public License for more details.
00017  *
00018  * You should have received a copy of the GNU General Public License
00019  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
00020  */
00021 
00022 #ifndef GKDATA_HPP
00023 #define GKDATA_HPP
00024 
00072 template <class Real>
00073 class GKData
00074 {
00075 private:
00076     Real eps_;  // machine epsilon
00077     
00079     Real abs_(Real x)
00080     {
00081         return (x < Real(0)) ? -x : x;
00082     }
00083     
00084     size_t m_;  // Gauss-Legendre degree
00085     size_t n_;  // size of Gauss-Kronrod arrays
00086     Real*  xgk_;  // Gauss-Kronrod abscissae
00087     Real*  wg_;   // Gauss-Legendre weights
00088     Real*  wgk_;  // Gauss-Kronrod weights
00089     
00090     Real *coefs;  // Chebyshev coefficients 
00091     Real *zeros;  // zeros of Legendre polynomial
00092     
00093     void legendre_zeros();
00094     void chebyshev_coefs();
00095     void gauss_kronrod_abscissae();
00096     void gauss_kronrod_weights();
00097     
00098     Real legendre_err(int deg, Real x, Real& err);
00099     Real legendre_deriv(int deg, Real x);   
00100     Real chebyshev_series(Real x, Real& err);
00101     Real chebyshev_series_deriv(Real x);
00102     
00103 public:
00105     GKData(Real MACHEPS, size_t m = 10);
00106     ~GKData();
00107     
00109     size_t size() { return n_; };
00110     
00112     Real xgk(int k) 
00113     { 
00114         return (0 <= k && k < n_) ? xgk_[k] : Real(0); 
00115     }
00116     
00118     Real wgk(int k)
00119     { 
00120         return (0 <= k && k < n_) ? wgk_[k] : Real(0); 
00121     }
00122     
00124     Real wg(int k) 
00125     { 
00126         return (0 <= k && k < n_/2) ? wg_[k] : Real(0); 
00127     }
00128 };
00129 
00130 template <class Real>
00131 GKData<Real>::GKData(Real MACHEPS, size_t m) : eps_(MACHEPS), m_(m)
00132 {
00133     n_ = m_ + 1;
00134     xgk_ = new Real[n_];
00135     wg_  = new Real[n_ / 2];
00136     wgk_ = new Real[n_];    
00137     coefs = new Real[n_ + 1];
00138     zeros = new Real[m_ + 2];
00139     
00140     legendre_zeros();
00141     chebyshev_coefs();
00142     gauss_kronrod_abscissae();
00143     gauss_kronrod_weights();
00144 }
00145 
00146 template <class Real>
00147 GKData<Real>::~GKData()
00148 {
00149     delete[] xgk_;
00150     delete[] wgk_;
00151     delete[] wg_;
00152     delete[] coefs;
00153     delete[] zeros;
00154 }
00155 
00161 template <class Real>
00162 void GKData<Real>::legendre_zeros()
00163 {
00164     Real* temp = new Real[m_+1];
00165     zeros[0] = Real(-1);
00166     zeros[1] = Real(1);
00167     Real delta, epsilon;
00168     
00169     for (int k = 1; k <= m_; ++k) {
00170         // loop to locate zeros of P_k interlacing z_0,...,z_k
00171         for (int j = 0; j < k; ++j) {
00172             // Newton's method for P_k :
00173             // initialize solver at midpoint of (z_j, z_{j+1})
00174             delta = 1;
00175             Real x_j = (zeros[j] + zeros[j+1]) / 2;
00176             Real P_k = legendre_err(k, x_j, epsilon);
00177             while (abs_(P_k) > epsilon && abs_(delta) > eps_) 
00178             {
00179                 delta = P_k / legendre_deriv(k, x_j);
00180                 x_j -= delta;
00181                 P_k = legendre_err(k, x_j, epsilon);
00182             }
00183             temp[j] = x_j;
00184         }
00185         
00186         // copy roots tmp_0,...,tmp_{k-1} to z_1,...z_k:
00187         zeros[k+1] = zeros[k];
00188         for (int j = 0; j < k; ++j)
00189             zeros[j+1] = temp[j];
00190         
00191     }
00192     delete[] temp;
00193 }
00194 
00198 template <class Real>
00199 void GKData<Real>::chebyshev_coefs()
00200 {
00201     size_t ell = (m_ + 1)/2;
00202     Real* alpha = new Real[ell+1];
00203     Real* f = new Real[ell+1];
00204     
00205     /* Care must be exercised in initalizing the constants in the definitions.
00206      * Compilers interpret expressions like "(2*k + 1.0)/(k + 2.0)" as double
00207      * precision before casting to Real.
00208      */
00209     f[1] = Real(m_+1)/Real(2*m_ + 3);
00210     alpha[0] = Real(1); // coefficient of T_{m+1}
00211     alpha[1] = -f[1];
00212     
00213     for (int k = 2; k <= ell; ++k) {
00214         f[k] = f[k-1] * (2*k - 1) * (m_ + k) / (k * (2*m_ + 2*k + 1));
00215         alpha[k] = -f[k];
00216         for (int i = 1; i < k; ++i)
00217             alpha[k] -= f[i] * alpha[k-i];
00218     }
00219     
00220     for (int k = 0; k <= ell; ++k) {
00221         coefs[m_ + 1 - 2*k] = alpha[k];
00222         if (m_  >= 2*k)
00223             coefs[m_ - 2*k] = Real(0);
00224     }
00225     
00226     delete[] alpha;
00227     delete[] f;
00228 }
00229 
00233 template <class Real>
00234 void GKData<Real>::gauss_kronrod_weights()
00235 {
00236     Real err;
00237     /* Gauss-Legendre weights:
00238      */
00239     for (int k = 0; k < n_ / 2; ++k) 
00240     {
00241         Real x = xgk_[2*k + 1];
00242         wg_[k] = (Real(-2) / 
00243                      ((m_ + 1) * legendre_deriv(m_, x) * legendre_err(m_+1, x, err)));
00244     }
00245     
00246     /* The ratio of leading coefficients of P_n and T_{n+1} is computed
00247      * from the recursive formulae for the respective polynomials.
00248      */
00249     Real F_m = Real(2) / Real(2*m_ + 1);
00250     for (int k = 1; k <= m_; ++k)
00251         F_m *= (Real(2*k) / Real(2*k - 1));
00252     
00253     /* Gauss-Kronrod weights:  
00254      */
00255     for (size_t k = 0; k < n_; ++k) 
00256     {
00257         Real x = xgk_[k];
00258         if (k % 2 == 0) 
00259         {
00260             wgk_[k] = F_m / (legendre_err(m_, x, err) * chebyshev_series_deriv(x));
00261         }
00262         else 
00263         {
00264             wgk_[k] = (wg_[k/2] + 
00265                           F_m / (legendre_deriv(m_, x) * chebyshev_series(x, err)));
00266         }
00267     }   
00268 }
00269 
00276 template <class Real>
00277 void GKData<Real>::gauss_kronrod_abscissae()
00278 {
00279     Real delta, epsilon;
00280     
00281     for (int k = 0; k < n_ / 2; ++k) 
00282     {
00283         delta = 1;
00284         // Newton's method for E_{n+1} :
00285         Real x_k = (zeros[m_-k] + zeros[m_+1-k])/Real(2);
00286         Real E = chebyshev_series(x_k, epsilon);
00287         while (abs_(E) > epsilon &&
00288                  abs_(delta) > eps_) 
00289         {
00290             delta = E / chebyshev_series_deriv(x_k);
00291             x_k -= delta;
00292             E = chebyshev_series(x_k, epsilon);
00293         }
00294         xgk_[2*k] = x_k;
00295         // copy adjacent Legendre-zero into the array:
00296         if (2*k+1 < n_)
00297             xgk_[2*k+1] = zeros[m_-k];
00298     }
00299 }
00300 
00310 template <class Real>
00311 Real GKData<Real>::legendre_err(int n, Real x, Real& err)
00312 {
00313     if (n == 0) {
00314         err = Real(0);
00315         return Real(1);
00316     }
00317     else if (n == 1) {
00318         err = Real(0);
00319         return x;
00320     }
00321     
00322     Real P0 = Real(1), P1 = x, P2;
00323     Real E0 = eps_; 
00324     Real E1 = abs_(x) * eps_; 
00325     for (int k = 1; k < n; ++k) 
00326     {
00327         P2 = ((2*k + 1) * x * P1 - k * P0) / (k + 1);
00328         err = ((2*k + 1) * abs_(x) * E1 + k * E0) / (2*(k + 1));
00329         P0 = P1; P1 = P2;
00330         E0 = E1; E1 = err;
00331     }
00332     return P2;  
00333 }
00334 
00340 template <class Real>
00341 Real GKData<Real>::legendre_deriv(int n, Real x)
00342 {
00343     if (n == 0)   
00344         return Real(0);
00345     else if (n == 1)
00346         return Real(1);
00347     
00348     Real P0 = Real(1), P1 = x, P2;
00349     Real dP0 = Real(0), dP1 = Real(1), dP2;
00350     for (int k = 1; k < n; ++k) 
00351     {
00352         P2 = ((2*k + 1) * x * P1 - k * P0) / (k + Real(1));
00353         dP2 = (2*k + 1) * P1 + dP0;
00354         P0 = P1; P1 = P2;
00355         dP0 = dP1; dP1 = dP2;
00356     }
00357     return dP2;
00358 }
00359 
00365 template <class Real>
00366 Real GKData<Real>::chebyshev_series(Real x, Real& err)
00367 {
00368     Real d1(0), d2(0);
00369     Real absc = abs_(coefs[0]); // final term for truncation error
00370     Real y2 = 2 * x; // linear term for Clenshaw recursion
00371     
00372     for (int k = n_; k >= 1; --k) {
00373         Real temp = d1;
00374         d1 = y2 * d1 - d2 + coefs[k];
00375       d2 = temp;
00376         absc += abs_(coefs[k]);
00377     }
00378     
00379     err = absc * eps_;
00380     return x * d1 - d2 + coefs[0]/2;
00381 }
00382 
00389 template <class Real>
00390 Real
00391 GKData<Real>::chebyshev_series_deriv(Real x)
00392 {
00393     Real d1(0), d2(0);
00394     Real y2 = 2 * x; // linear term for Clenshaw recursion
00395     
00396     for (int k = n_; k >= 2; --k) {
00397         Real temp = d1;
00398         d1 = y2 * d1 - d2 + k * coefs[k];
00399       d2 = temp;
00400     }
00401     
00402     return y2 * d1 - d2 + coefs[1];
00403 }
00404 
00405 #endif // GKDATA_HPP
 All Classes Functions Variables