libs/math/doc/sf/chebyshev.qbk

   1 [/
   2   Copyright 2017, Nick Thompson
   3   Distributed under the Boost Software License, Version 1.0.
   4   (See accompanying file LICENSE_1_0.txt or copy at
   5   http://www.boost.org/LICENSE_1_0.txt).
   6 ]
   7
   8 [section:chebyshev Chebyshev Polynomials]
   9
  10 [h4 Synopsis]
  11
  12 ``
  13 #include <boost/math/special_functions/chebyshev.hpp>
  14 ``
  15
  16    namespace boost{ namespace math{
  17
  18    template<class Real1, class Real2, class Real3>
  19    ``__sf_result`` chebyshev_next(Real1 const & x, Real2 const & Tn, Real3 const & Tn_1);
  20
  21    template<class Real>
  22    ``__sf_result`` chebyshev_t(unsigned n, Real const & x);
  23
  24    template<class Real, class ``__Policy``>
  25    ``__sf_result`` chebyshev_t(unsigned n, Real const & x, const ``__Policy``&);
  26
  27    template<class Real>
  28    ``__sf_result`` chebyshev_u(unsigned n, Real const & x);
  29
  30    template<class Real, class ``__Policy``>
  31    ``__sf_result`` chebyshev_u(unsigned n, Real const & x, const ``__Policy``&);
  32
  33    template<class Real>
  34    ``__sf_result`` chebyshev_t_prime(unsigned n, Real const & x);
  35
  36    template<class Real1, class Real2>
  37    ``__sf_result`` chebyshev_clenshaw_recurrence(const Real* const c, size_t length, Real2 x);
  38
  39    }} // namespaces
  40
  41
  42 ['"Real analysts cannot do without Fourier, complex analysts cannot do without Laurent, and numerical analysts cannot do without Chebyshev"] --Lloyd N. Trefethen
  43
  44 The Chebyshev polynomials of the first kind are defined by the recurrence /T/[sub n+1](/x/) := /2xT/[sub n](/x/) - /T/[sub n-1](/x/), /n > 0/,
  45 where /T/[sub 0](/x/) := 1 and /T/[sub 1](/x/) := /x/.
  46 These can be calculated in Boost using the following simple code
  47
  48     double x = 0.5;
  49     double T12 = boost::math::chebyshev_t(12, x);
  50
  51 Calculation of derivatives is also straightforward:
  52
  53     double T12_prime = boost::math::chebyshev_t_prime(12, x);
  54
  55 The complexity of evaluation of the /n/-th Chebyshev polynomial by these functions is linear.
  56 So they are unsuitable for use in calculation of (say) a Chebyshev series, as a sum of linear scaling functions scales quadratically.
  57 Though there are very sophisticated algorithms for the evaluation of Chebyshev series,
  58 a linear time algorithm is presented below:
  59
  60     double x = 0.5;
  61     std::vector<double> c{14.2, -13.7, 82.3, 96};
  62     double T0 = 1;
  63     double T1 = x;
  64     double f = c[0]*T0/2;
  65     unsigned l = 1;
  66     while(l < c.size())
  67     {
  68        f += c[l]*T1;
  69        std::swap(T0, T1);
  70        T1 = boost::math::chebyshev_next(x, T0, T1);
  71        ++l;
  72     }
  73
  74 This uses the `chebyshev_next` function to evaluate each term of the Chebyshev series in constant time.
  75 However, this naive algorithm has a catastrophic loss of precision as /x/ approaches 1.
  76 A method to mitigate this way given by [@http://www.ams.org/journals/mcom/1955-09-051/S0025-5718-1955-0071856-0/S0025-5718-1955-0071856-0.pdf Clenshaw],
  77 and is implemented in boost as
  78
  79     double x = 0.5;
  80     std::vector<double> c{14.2, -13.7, 82.3, 96};
  81     double f = chebyshev_clenshaw_recurrence(c.data(), c.size(), Real x);
  82
  83
  84 N.B.: There is factor of /2/ difference in our definition of the first coefficient in the Chebyshev series from Clenshaw's original work.
  85 This is because two traditions exist in notation for the Chebyshev series expansion,
  86
  87 [:/f/(/x/) \u2248 \u2211[sub n=0][super N-1] /a/[sub n]/T/[sub n](/x/)]
  88
  89 and
  90
  91 [:/f/(/x/) \u2248 /c/[sub 0]/2 + \u2211[sub n=1][super N-1] /c/[sub n]/T/[sub n](/x/)]
  92
  93 ['*boost math always uses the second convention, with the factor of 1/2 on the first coefficient.*]
  94
  95 Chebyshev polynomials of the second kind can be evaluated via `chebyshev_u`:
  96
  97     double x = -0.23;
  98     double U1 = boost::math::chebyshev_u(1, x);
  99
 100 The evaluation of Chebyshev polynomials by a three-term recurrence is known to be
 101 [@https://link.springer.com/article/10.1007/s11075-014-9925-x mixed forward-backward stable] for /x/ \u220A \[-1, 1\].
 102 However, the author does not know of a similar result for /x/ outside \[-1, 1\].
 103 For this reason, evaluation of Chebyshev polynomials outside of \[-1, 1\] is strongly discouraged.
 104 That said, small rounding errors in the course of a computation will often lead to this situation,
 105 and termination of the computation due to these small problems is very discouraging.
 106 For this reason, `chebyshev_t` and `chebyshev_u` have code paths for /x > 1/ and /x < -1/ which do not use three-term recurrences.
 107 These code paths are /much slower/, and should be avoided if at all possible.
 108
 109 Evaluation of a Chebyshev series is relatively simple.
 110 The real challenge is /generation/ of the Chebyshev series.
 111 For this purpose, boost provides a /Chebyshev transform/, a projection operator which projects a function onto a finite-dimensional span of Chebyshev polynomials.
 112 But before we discuss the API, let's analyze why we might want to project a function onto a span of Chebyshev polynomials.
 113
 114 * We want a numerically stable way to evaluate the function's derivative.
 115 * Our function is expensive to evaluate, and we wish to find a less expensive way to estimate its value.
 116 An example are the standard library transcendental functions:
 117 These functions are guaranteed to evaluate to within 1 ulp of the exact value, but often this accuracy is not needed.
 118 A projection onto the Chebyshev polynomials with a low accuracy requirement can vastly accelerate the computation of these functions.
 119 * We wish to numerically integrate the function.
 120
 121 The API is given below.
 122
 123 ``
 124 #include <boost/math/special_functions/chebyshev_transform.hpp>
 125 ``
 126
 127    namespace boost{ namespace math{
 128
 129    template<class Real>
 130    class chebyshev_transform
 131    {
 132    public:
 133        template<class F>
 134        chebyshev_transform(const F& f, Real a, Real b, Real tol=500*std::numeric_limits<Real>::epsilon());
 135
 136        Real operator()(Real x) const
 137
 138        Real integrate() const
 139
 140        const std::vector<Real>& coefficients() const
 141
 142        Real prime(Real x) const
 143    };
 144
 145    }}// end namespaces
 146
 147
 148 The Chebyshev transform takes a function /f/ and returns a /near-minimax/ approximation to /f/ in terms of Chebyshev polynomials.
 149 By /near-minimax/, we mean that the resulting Chebyshev polynomial is "very close" the polynomial /p/[sub n]  which minimizes the uniform norm of /f/ - /p/[sub n].
 150 The notion of "very close" can be made rigorous; see Trefethen's "Approximation Theory and Approximation Practice" for details.
 151
 152 The Chebyshev transform works by creating a vector of values by evaluating the input function at the Chebyshev points, and then performing a discrete cosine transform on the resulting vector.
 153 In order to do this efficiently, we have used [@http://www.fftw.org/ FFTW3].
 154 So to compile, you must have `FFTW3` installed, and link with `-lfftw3` for double precision, `-lfftw3f` for float precision, `-lfftw3l` for long double precision, and -lfftwq for quad (`__float128`) precision.
 155 After the coefficients of the Chebyshev series are known, the routine goes back through them and filters out all the coefficients whose absolute ratio to the largest coefficient are less than the tolerance requested in the constructor.
 156
 157 [endsect] [/section:chebyshev Chebyshev Polynomials]
 158