libs/math/doc/interpolators/cardinal_trigonometric.qbk

   1 [/
   2 Copyright (c) 2019 Nick Thompson
   3 Use, modification and distribution are subject to the
   4 Boost Software License, Version 1.0. (See accompanying file
   5 LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6 ]
   7
   8 [section:cardinal_trigonometric Cardinal Trigonometric interpolation]
   9
  10 [heading Synopsis]
  11
  12 ```
  13 #include <boost/math/interpolators/cardinal_trigonometric.hpp>
  14
  15 namespace boost{ namespace math{ namespace interpolators {
  16
  17 template <class RandomAccessContainer>
  18 class cardinal_trigonometric
  19 {
  20 public:
  21     cardinal_trigonometric(RandomAccessContainer const & y, Real t0, Real h);
  22
  23     Real operator()(Real t) const;
  24
  25     Real prime(Real t) const;
  26
  27     Real double_prime(Real t) const;
  28
  29     Real period() const;
  30
  31     Real integrate() const;
  32
  33     Real squared_l2() const;
  34 };
  35 }}}
  36 ```
  37
  38 [heading Cardinal Trigonometric Interpolation]
  39
  40 The cardinal trigonometric interpolation problem takes uniformly spaced samples /y/[sub j] of a periodic function /f/ defined via /y/[sub /j/] = /f/(/t/[sub 0] + /j/ /h/) and represents them as a linear combination of sines and cosines.
  41 If the period of /f/ is /T/, and the number of data points is /n = 2m/ or /n = 2m+1/, we hope to have
  42
  43 /f/(/t/) \u2248 /a/[sub 0]/2 + \u2211[sub /k/ = 1][super /m/] /a/[sub /k/] cos(2\u03C0 /k/ (/t/-/t/[sub 0]) \/T) + /b/[sub /k/] sin(2\u03C0 /k/ (/t/-/t/[sub 0])/T)
  44
  45 Convergence rates depend on the number of continuous derivatives of /f/; see either Atkinson or Kress for details.
  46
  47 A simple use of this interpolator is shown below:
  48
  49 ```
  50 #include <vector>
  51 #include <boost/math/interpolators/cardinal_trigonometric.hpp>
  52 using boost::math::interpolators::cardinal_trigonometric;
  53 ...
  54 std::vector<double> v(17, 3.2);
  55 auto ct = cardinal_trigonometric(v, /*t0 = */ 0.0, /* h = */ 0.125);
  56 std::cout << "ct(1.3) = " << ct(1.3) << "\n";
  57
  58 // Derivative:
  59 std::cout << ct.prime(1.2) << "\n";
  60 // Second derivative:
  61 std::cout << ct.double_prime(1.2) << "\n";
  62 ```
  63
  64 The period is always given by `v.size()*h`.
  65 Off-by-one errors are common in programming, and hence if you wonder what this interpolator believes the period to be, you can query it with the `.period()` member function.
  66
  67 In addition, the transform into the trigonometric basis gives a trivial way to compute the integral of the function over a period; this is done via the `.integrate()` member function.
  68 Evaluation of the square of the L[super 2] norm is trivial in this basis; it is computed by the `.squared_l2()` member function.
  69
  70 Below is a graph of a /C/[super \u221E] bump function approximated by trigonometric series.
  71 The graphs are visually indistinguishable at 20 samples.
  72
  73 [$../graphs/fourier_bump.svg]
  74
  75
  76 [heading Caveats]
  77
  78 This routine depends on FFTW3, and hence will only compile in float, double, long double, and quad precision, unlike the large bulk of the library which is compatible with arbitrary precision arithmetic.
  79 The FFTW linker flags must be added to the compile step, i.e., `-lm -lfftw3` for double precision, `-lm -lfftw3f` for float, so on.
  80
  81 Evaluation of derivatives is done by differentiation of Horner's method.
  82 As always, differentiation amplifies noise; and because some rounding error is produced by computation of the Fourier coefficients, this error is amplified by differentiation.
  83
  84
  85 [heading References]
  86
  87 * Atkinson, Kendall, and Weimin Han. ['Theoretical numerical analysis.] Vol. 39. Berlin: Springer, 2005.
  88
  89 * Kress, Rainer. ['Numerical Analysis.] 1998. Academic Edition 1.
  90
  91 [endsect]
  92 [/section:cardinal_trigonometric]