intccq.cpp Source File

00001 /*=============================================================================
00002         File: intccq.cpp
00003      Purpose:       
00004     Revision: $Id: intccq.cpp,v 1.2 2002/05/13 21:07:45 philosophil Exp $
00005   Created by: Philippe Lavoie          (3 Oct, 1996)
00006  Modified by: 
00007 
00008  Copyright notice:
00009           Copyright (C) 1996-1998 Philippe Lavoie
00010  
00011           This library is free software; you can redistribute it and/or
00012           modify it under the terms of the GNU Library General Public
00013           License as published by the Free Software Foundation; either
00014           version 2 of the License, or (at your option) any later version.
00015  
00016           This library is distributed in the hope that it will be useful,
00017           but WITHOUT ANY WARRANTY; without even the implied warranty of
00018           MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019           Library General Public License for more details.
00020  
00021           You should have received a copy of the GNU Library General Public
00022           License along with this library; if not, write to the Free
00023           Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00024 =============================================================================*/
00025 
00026 #ifndef intccqq_SOURCES
00027 #define intccqq_SOURCES
00028 
00029 
00030 /*
00031 Clenshaw-Curtis-Quadrature
00032 Numerical Automatic Integrator
00033     method    : Chebyshev Series Expansion
00034     dimension : one
00035     table     : use
00036 function
00037     intcc  : integrator of f(x) over [a,b].
00038 necessary package
00039     fft2f.c  : FFT package
00040 */
00041 
00042 /*
00043 intcc
00044     [description]
00045         I = integral of f(x) over [a,b]
00046     [declaration]
00047         void intccini(int lenw, double *w);
00048         void intcc(double (*f)(double), double a, double b, double eps, 
00049             int lenw, double *w, double *i, double *err);
00050     [usage]
00051         intccini(lenw, w);  // initialization of w
00052         ...
00053         intcc(f, a, b, eps, lenw, w, &i, &err);
00054     [parameters]
00055         lenw      : (length of w[]) - 1 (int)
00056         w         : work area and weights of the quadrature 
00057                     formula, w[0...lenw] (double *)
00058         f         : integrand f(x) (double (*f)(double))
00059         a         : lower limit of integration (double)
00060         b         : upper limit of integration (double)
00061         eps       : relative error requested (double)
00062         i         : approximation to the integral (double *)
00063         err       : estimate of the absolute error (double *)
00064     [remarks]
00065         initial parameters
00066             lenw >= 14 and 
00067             lenw > (maximum number of f(x) evaluations) * 3 / 2
00068             example :
00069                 lenc = 3200;
00070         function
00071             f(x) needs to be analytic over [a,b].
00072         relative error
00073             eps is relative error requested excluding 
00074             cancellation of significant digits.
00075             i.e. eps means : (absolute error) / 
00076                              (integral_a^b |f(x)| dx).
00077             eps does not mean : (absolute error) / I.
00078         error message
00079             err >= 0 : normal termination.
00080             err < 0  : abnormal termination (n > nmax).
00081                        i.e. convergent error is detected :
00082                            1. f(x) or (d/dx)^n f(x) has 
00083                               discontinuous points or sharp 
00084                               peaks over [a,b].
00085                               you must use other routine.
00086                            2. relative error of f(x) is 
00087                               greater than eps.
00088                            3. f(x) has oscillatory factor 
00089                               and frequency of the oscillation 
00090                               is very high.
00091 */
00092 
00093 
00094 #include "integrate.h"
00095 #include <math.h>
00096 
00099 namespace PLib {
00100 
00111 template <class T>
00112 void intccini(BasicArray<T> &w)
00113 {
00114   int lenw = w.n() - 1 ; 
00115   int j, k, l, m;
00116   T cos2, sin1, sin2, hl;
00117     
00118   cos2 = 0;
00119   sin1 = 1;
00120   sin2 = 1;
00121   hl = T(0.5);
00122   k = lenw;
00123   l = 2;
00124   while (l < k - l - 1) {
00125     w[0] = hl * T(0.5);
00126     for (j = 1; j <= l; j++) {
00127       w[j] = hl / (1 - 4 * j * j);
00128     }
00129     w[l] *= T(0.5);
00130     dfct(l, T(0.5) * cos2, sin1, w);
00131     cos2 = sqrt(2 + cos2);
00132     sin1 /= cos2;
00133     sin2 /= 2 + cos2;
00134     w[k] = sin2;
00135     w[k - 1] = w[0];
00136     w[k - 2] = w[l];
00137     k -= 3;
00138     m = l;
00139     while (m > 1) {
00140       m >>= 1;
00141       for (j = m; j <= l - m; j += (m << 1)) {
00142         w[k] = w[j];
00143         k--;
00144       }
00145     }
00146     hl *= T(0.5);
00147     l *= 2;
00148   }
00149 }
00150 
00151 
00207 template <class T, class OPPtr>
00208 T intcc(OPPtr f, T a, T b, T eps, BasicArray<T> &w, T& err)
00209 {
00210   T i ; 
00211   int j, k, l;
00212   T esf, eref, erefh, hh, ir, iback, irback, ba, ss, x, y, fx, 
00213     errir;
00214   int lenw = w.n() - 1 ;
00215   esf = 10;
00216   ba = T(0.5) * (b - a);
00217   ss = 2 * w[lenw];
00218   x = ba * w[lenw];
00219   w[0] = T(0.5) * (*f)(a);
00220   w[3] = T(0.5) * (*f)(b);
00221   w[2] = (*f)(a + x);
00222   w[4] = (*f)(b - x);
00223   w[1] = (*f)(a + ba);
00224   eref = T(0.5) * (fabs(w[0]) + fabs(w[1]) + fabs(w[2]) + fabs(w[3]) + 
00225                    fabs(w[4]));
00226   w[0] += w[3];
00227   w[2] += w[4];
00228   ir = w[0] + w[1] + w[2];
00229   i = w[0] * w[lenw - 1] + w[1] * w[lenw - 2] + w[2] * w[lenw - 3];
00230   erefh = eref * sqrt(eps);
00231   eref *= eps;
00232   hh = T(0.25);
00233   l = 2;
00234   k = lenw - 5;
00235   do {
00236     iback = i;
00237     irback = ir;
00238     x = ba * w[k + 1];
00239     y = 0;
00240     i = w[0] * w[k];
00241     for (j = 1; j <= l; j++) {
00242       x += y;
00243       y += ss * (ba - x);
00244       fx = (*f)(a + x) + (*f)(b - x);
00245       ir += fx;
00246       i += w[j] * w[k - j] + fx * w[k - j - l];
00247       w[j + l] = fx;
00248     }
00249     ss = 2 * w[k + 1];
00250     err = esf * l * fabs(i - iback);
00251     hh *= T(0.25);
00252     errir = hh * fabs(ir - 2 * irback);
00253     l *= 2;
00254     k -= l + 2;
00255   } while ((err > erefh || errir > eref) && k > 4 * l);
00256   i *= b - a;
00257   if (err > erefh || errir > eref) {
00258     err *= -fabs(b - a);
00259   } else {
00260     err = eref * fabs(b - a);
00261   }
00262   return i ;
00263 }
00264 
00320 template <class T, class OPPtr>
00321 T intcc2(OPPtr f, T a, T b, T eps, BasicArray<T> w, T& err)
00322 {
00323   T i ; 
00324   int j, k, l;
00325   T esf, eref, erefh, hh, ir, iback, irback, ba, ss, x, y, fx, 
00326     errir;
00327   int lenw = w.n() - 1 ;
00328   esf = 10;
00329   ba = T(0.5) * (b - a);
00330   ss = 2 * w[lenw];
00331   x = ba * w[lenw];
00332   w[0] = T(0.5) * (*f)(a);
00333   w[3] = T(0.5) * (*f)(b);
00334   w[2] = (*f)(a + x);
00335   w[4] = (*f)(b - x);
00336   w[1] = (*f)(a + ba);
00337   eref = T(0.5) * (fabs(w[0]) + fabs(w[1]) + fabs(w[2]) + fabs(w[3]) + 
00338                    fabs(w[4]));
00339   w[0] += w[3];
00340   w[2] += w[4];
00341   ir = w[0] + w[1] + w[2];
00342   i = w[0] * w[lenw - 1] + w[1] * w[lenw - 2] + w[2] * w[lenw - 3];
00343   erefh = eref * sqrt(eps);
00344   eref *= eps;
00345   hh = T(0.25);
00346   l = 2;
00347   k = lenw - 5;
00348   do {
00349     iback = i;
00350     irback = ir;
00351     x = ba * w[k + 1];
00352     y = 0;
00353     i = w[0] * w[k];
00354     for (j = 1; j <= l; j++) {
00355       x += y;
00356       y += ss * (ba - x);
00357       fx = (*f)(a + x) + (*f)(b - x);
00358       ir += fx;
00359       i += w[j] * w[k - j] + fx * w[k - j - l];
00360       w[j + l] = fx;
00361     }
00362     ss = 2 * w[k + 1];
00363     err = esf * l * fabs(i - iback);
00364     hh *= T(0.25);
00365     errir = hh * fabs(ir - 2 * irback);
00366     l *= 2;
00367     k -= l + 2;
00368   } while ((err > erefh || errir > eref) && k > 4 * l);
00369   i *= b - a;
00370   if (err > erefh || errir > eref) {
00371     err *= -fabs(b - a);
00372   } else {
00373     err = eref * fabs(b - a);
00374   }
00375   return i ;
00376 }
00377 
00378 template <class T, class OPPtr> 
00379 T integrate(OPPtr f, T a, T b, T eps, int n, T &err) {
00380   BasicArray<T> w(n) ; 
00381   intccini(w) ; 
00382   return intcc(f,a,b,eps,w,err) ; 
00383 }
00384 
00385 template <class T, class OPPtr> 
00386 T integrate2(OPPtr f, T a, T b, T eps, int n, T &err) {
00387   BasicArray<T> w(n) ; 
00388   intccini(w) ; 
00389   return intcc2(f,a,b,eps,w,err) ; 
00390 }
00391 
00447 template <class T, class OPvPtr>
00448 T intcc(OPvPtr f, void* userData, T a, T b, T eps, BasicArray<T> &w, T& err)
00449 {
00450   T i ; 
00451   int j, k, l;
00452   T esf, eref, erefh, hh, ir, iback, irback, ba, ss, x, y, fx, 
00453     errir;
00454   int lenw = w.n() - 1 ;
00455   esf = 10;
00456   ba = T(0.5) * (b - a);
00457   ss = 2 * w[lenw];
00458   x = ba * w[lenw];
00459   w[0] = T(0.5) * (*f)(a,userData);
00460   w[3] = T(0.5) * (*f)(b,userData);
00461   w[2] = (*f)(a + x,userData);
00462   w[4] = (*f)(b - x,userData);
00463   w[1] = (*f)(a + ba,userData);
00464   eref = T(0.5) * (fabs(w[0]) + fabs(w[1]) + fabs(w[2]) + fabs(w[3]) + 
00465                    fabs(w[4]));
00466   w[0] += w[3];
00467   w[2] += w[4];
00468   ir = w[0] + w[1] + w[2];
00469   i = w[0] * w[lenw - 1] + w[1] * w[lenw - 2] + w[2] * w[lenw - 3];
00470   erefh = eref * sqrt(eps);
00471   eref *= eps;
00472   hh = T(0.25);
00473   l = 2;
00474   k = lenw - 5;
00475   do {
00476     iback = i;
00477     irback = ir;
00478     x = ba * w[k + 1];
00479     y = 0;
00480     i = w[0] * w[k];
00481     for (j = 1; j <= l; j++) {
00482       x += y;
00483       y += ss * (ba - x);
00484       fx = (*f)(a + x,userData) + (*f)(b - x,userData);
00485       ir += fx;
00486       i += w[j] * w[k - j] + fx * w[k - j - l];
00487       w[j + l] = fx;
00488     }
00489     ss = 2 * w[k + 1];
00490     err = esf * l * fabs(i - iback);
00491     hh *= T(0.25);
00492     errir = hh * fabs(ir - 2 * irback);
00493     l *= 2;
00494     k -= l + 2;
00495   } while ((err > erefh || errir > eref) && k > 4 * l);
00496   i *= b - a;
00497   if (err > erefh || errir > eref) {
00498     err *= -fabs(b - a);
00499   } else {
00500     err = eref * fabs(b - a);
00501   }
00502   return i ;
00503 }
00504 
00560 template <class T, class OPvPtr>
00561 T intcc2(OPvPtr f, void* userData, T a, T b, T eps, BasicArray<T> w, T& err)
00562 {
00563   T i ; 
00564   int j, k, l;
00565   T esf, eref, erefh, hh, ir, iback, irback, ba, ss, x, y, fx, 
00566     errir;
00567   int lenw = w.n() - 1 ;
00568   esf = 10;
00569   ba = T(0.5) * (b - a);
00570   ss = 2 * w[lenw];
00571   x = ba * w[lenw];
00572   w[0] = T(0.5) * (*f)(a,userData);
00573   w[3] = T(0.5) * (*f)(b,userData);
00574   w[2] = (*f)(a + x,userData);
00575   w[4] = (*f)(b - x,userData);
00576   w[1] = (*f)(a + ba,userData);
00577   eref = T(0.5) * (fabs(w[0]) + fabs(w[1]) + fabs(w[2]) + fabs(w[3]) + 
00578                    fabs(w[4]));
00579   w[0] += w[3];
00580   w[2] += w[4];
00581   ir = w[0] + w[1] + w[2];
00582   i = w[0] * w[lenw - 1] + w[1] * w[lenw - 2] + w[2] * w[lenw - 3];
00583   erefh = eref * sqrt(eps);
00584   eref *= eps;
00585   hh = T(0.25);
00586   l = 2;
00587   k = lenw - 5;
00588   do {
00589     iback = i;
00590     irback = ir;
00591     x = ba * w[k + 1];
00592     y = 0;
00593     i = w[0] * w[k];
00594     for (j = 1; j <= l; j++) {
00595       x += y;
00596       y += ss * (ba - x);
00597       fx = (*f)(a + x,userData) + (*f)(b - x,userData);
00598       ir += fx;
00599       i += w[j] * w[k - j] + fx * w[k - j - l];
00600       w[j + l] = fx;
00601     }
00602     ss = 2 * w[k + 1];
00603     err = esf * l * fabs(i - iback);
00604     hh *= T(0.25);
00605     errir = hh * fabs(ir - 2 * irback);
00606     l *= 2;
00607     k -= l + 2;
00608   } while ((err > erefh || errir > eref) && k > 4 * l);
00609   i *= b - a;
00610   if (err > erefh || errir > eref) {
00611     err *= -fabs(b - a);
00612   } else {
00613     err = eref * fabs(b - a);
00614   }
00615   return i ;
00616 }
00617 
00618 template <class T, class OPvPtr> 
00619 T integrate(OPvPtr f, void* userData, T a, T b, T eps, int n, T &err) {
00620   BasicArray<T> w(n) ; 
00621   intccini(w) ; 
00622   return intcc(f,userData,a,b,eps,w,err) ; 
00623 }
00624 
00625 template <class T, class OPvPtr> 
00626 T integrate2(OPvPtr f, void* userData, T a, T b, T eps, int n, T &err) {
00627   BasicArray<T> w(n) ; 
00628   intccini(w) ; 
00629   return intcc2(f,userData,a,b,eps,w,err) ; 
00630 }
00631 
00632 
00633 
00634 }
00635 
00636 
00637 #endif