libs/math/test/chebyshev_transform_test.cpp

   1 /*
   2  * Copyright Nick Thompson, 2017
   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 #define BOOST_TEST_MODULE chebyshev_transform_test
   8
   9 #include <boost/cstdfloat.hpp>
  10 #include <boost/type_index.hpp>
  11 #include <boost/test/included/unit_test.hpp>
  12 #include <boost/test/tools/floating_point_comparison.hpp>
  13 #include <boost/math/special_functions/chebyshev.hpp>
  14 #include <boost/math/special_functions/chebyshev_transform.hpp>
  15 #include <boost/math/special_functions/sinc.hpp>
  16 #include <boost/multiprecision/cpp_bin_float.hpp>
  17 #include <boost/multiprecision/cpp_dec_float.hpp>
  18
  19 #if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4)
  20 #  define TEST1
  21 #  define TEST2
  22 #  define TEST3
  23 #  define TEST4
  24 #endif
  25
  26 using boost::multiprecision::cpp_bin_float_quad;
  27 using boost::multiprecision::cpp_bin_float_50;
  28 using boost::multiprecision::cpp_bin_float_100;
  29 using boost::math::chebyshev_t;
  30 using boost::math::chebyshev_t_prime;
  31 using boost::math::chebyshev_u;
  32 using boost::math::chebyshev_transform;
  33
  34
  35 template<class Real>
  36 void test_sin_chebyshev_transform()
  37 {
  38     using boost::math::chebyshev_transform;
  39     using boost::math::constants::half_pi;
  40     using std::sin;
  41     using std::cos;
  42     using std::abs;
  43
  44     Real tol = 10*std::numeric_limits<Real>::epsilon();
  45     auto f = [](Real x) { return sin(x); };
  46     Real a = 0;
  47     Real b = 1;
  48     chebyshev_transform<Real> cheb(f, a, b, tol);
  49
  50     Real x = a;
  51     while (x < b)
  52     {
  53         Real s = sin(x);
  54         Real c = cos(x);
  55         if (abs(s) < tol)
  56         {
  57             BOOST_CHECK_SMALL(cheb(x), 100*tol);
  58             BOOST_CHECK_CLOSE_FRACTION(c, cheb.prime(x), 100*tol);
  59         }
  60         else
  61         {
  62             BOOST_CHECK_CLOSE_FRACTION(s, cheb(x), 100*tol);
  63             if (abs(c) < tol)
  64             {
  65                 BOOST_CHECK_SMALL(cheb.prime(x), 100*tol);
  66             }
  67             else
  68             {
  69                 BOOST_CHECK_CLOSE_FRACTION(c, cheb.prime(x), 100*tol);
  70             }
  71         }
  72         x += static_cast<Real>(1)/static_cast<Real>(1 << 7);
  73     }
  74
  75     Real Q = cheb.integrate();
  76
  77     BOOST_CHECK_CLOSE_FRACTION(1 - cos(static_cast<Real>(1)), Q, 100*tol);
  78 }
  79
  80
  81 template<class Real>
  82 void test_sinc_chebyshev_transform()
  83 {
  84     using std::cos;
  85     using std::sin;
  86     using std::abs;
  87     using boost::math::sinc_pi;
  88     using boost::math::chebyshev_transform;
  89     using boost::math::constants::half_pi;
  90
  91     Real tol = 500*std::numeric_limits<Real>::epsilon();
  92     auto f = [](Real x) { return boost::math::sinc_pi(x); };
  93     Real a = 0;
  94     Real b = 1;
  95     chebyshev_transform<Real> cheb(f, a, b, tol/50);
  96
  97     Real x = a;
  98     while (x < b)
  99     {
 100         Real s = sinc_pi(x);
 101         Real ds = (cos(x)-sinc_pi(x))/x;
 102         if (x == 0) { ds = 0; }
 103         if (s < tol)
 104         {
 105             BOOST_CHECK_SMALL(cheb(x), tol);
 106         }
 107         else
 108         {
 109             BOOST_CHECK_CLOSE_FRACTION(s, cheb(x), tol);
 110         }
 111
 112         if (abs(ds) < tol)
 113         {
 114             BOOST_CHECK_SMALL(cheb.prime(x), 5 * tol);
 115         }
 116         else
 117         {
 118             BOOST_CHECK_CLOSE_FRACTION(ds, cheb.prime(x), 300*tol);
 119         }
 120         x += static_cast<Real>(1)/static_cast<Real>(1 << 7);
 121     }
 122
 123     Real Q = cheb.integrate();
 124     //NIntegrate[Sinc[x], {x, 0, 1}, WorkingPrecision -> 200, AccuracyGoal -> 150, PrecisionGoal -> 150, MaxRecursion -> 150]
 125     Real Q_exp = boost::lexical_cast<Real>("0.94608307036718301494135331382317965781233795473811179047145477356668");
 126     BOOST_CHECK_CLOSE_FRACTION(Q_exp, Q, tol);
 127 }
 128
 129
 130
 131 //Examples taken from "Approximation Theory and Approximation Practice", by Trefethen
 132 template<class Real>
 133 void test_atap_examples()
 134 {
 135     using std::sin;
 136     using boost::math::constants::half;
 137     using boost::math::sinc_pi;
 138     using boost::math::chebyshev_transform;
 139     using boost::math::constants::half_pi;
 140
 141     Real tol = 10*std::numeric_limits<Real>::epsilon();
 142     auto f1 = [](Real x) { return ((0 < x) - (x < 0)) - x/2; };
 143     auto f2 = [](Real x) { Real t = sin(6*x); Real s = sin(x + exp(2*x));
 144                            Real u = (0 < s) - (s < 0);
 145                            return t + u; };
 146
 147     auto f3 = [](Real x) { return sin(6*x) + sin(60*exp(x)); };
 148
 149     auto f4 = [](Real x) { return 1/(1+1000*(x+half<Real>())*(x+half<Real>())) + 1/sqrt(1+1000*(x-.5)*(x-0.5));};
 150     Real a = -1;
 151     Real b = 1;
 152     chebyshev_transform<Real> cheb1(f1, a, b);
 153     chebyshev_transform<Real> cheb2(f2, a, b, tol);
 154     //chebyshev_transform<Real> cheb3(f3, a, b, tol);
 155
 156     Real x = a;
 157     while (x < b)
 158     {
 159         //Real s = f1(x);
 160         if (sizeof(Real) == sizeof(float))
 161         {
 162            BOOST_CHECK_CLOSE_FRACTION(f1(x), cheb1(x), 4e-3);
 163         }
 164         else
 165         {
 166            BOOST_CHECK_CLOSE_FRACTION(f1(x), cheb1(x), 1.3e-5);
 167         }
 168         BOOST_CHECK_CLOSE_FRACTION(f2(x), cheb2(x), 6e-3);
 169         //BOOST_CHECK_CLOSE_FRACTION(f3(x), cheb3(x), 100*tol);
 170         x += static_cast<Real>(1)/static_cast<Real>(1 << 7);
 171     }
 172 }
 173
 174 //Validate that the Chebyshev polynomials are well approximated by the Chebyshev transform.
 175 template<class Real>
 176 void test_chebyshev_chebyshev_transform()
 177 {
 178     Real tol = 500*std::numeric_limits<Real>::epsilon();
 179     // T_0 = 1:
 180     auto t0 = [](Real) { return 1; };
 181     chebyshev_transform<Real> cheb0(t0, -1, 1);
 182     BOOST_CHECK_CLOSE_FRACTION(cheb0.coefficients()[0], 2, tol);
 183
 184     Real x = -1;
 185     while (x < 1)
 186     {
 187         BOOST_CHECK_CLOSE_FRACTION(cheb0(x), 1, tol);
 188         BOOST_CHECK_SMALL(cheb0.prime(x), tol);
 189         x += static_cast<Real>(1)/static_cast<Real>(1 << 7);
 190     }
 191
 192     // T_1 = x:
 193     auto t1 = [](Real x) { return x; };
 194     chebyshev_transform<Real> cheb1(t1, -1, 1);
 195     BOOST_CHECK_CLOSE_FRACTION(cheb1.coefficients()[1], 1, tol);
 196
 197     x = -1;
 198     while (x < 1)
 199     {
 200         if (x == 0)
 201         {
 202             BOOST_CHECK_SMALL(cheb1(x), tol);
 203         }
 204         else
 205         {
 206             BOOST_CHECK_CLOSE_FRACTION(cheb1(x), x, tol);
 207         }
 208         BOOST_CHECK_CLOSE_FRACTION(cheb1.prime(x), 1, tol);
 209         x += static_cast<Real>(1)/static_cast<Real>(1 << 7);
 210     }
 211
 212
 213     auto t2 = [](Real x) { return 2*x*x-1; };
 214     chebyshev_transform<Real> cheb2(t2, -1, 1);
 215     BOOST_CHECK_CLOSE_FRACTION(cheb2.coefficients()[2], 1, tol);
 216
 217     x = -1;
 218     while (x < 1)
 219     {
 220         BOOST_CHECK_CLOSE_FRACTION(cheb2(x), t2(x), tol);
 221         if (x != 0)
 222         {
 223             BOOST_CHECK_CLOSE_FRACTION(cheb2.prime(x), 4*x, tol);
 224         }
 225         else
 226         {
 227             BOOST_CHECK_SMALL(cheb2.prime(x), tol);
 228         }
 229         x += static_cast<Real>(1)/static_cast<Real>(1 << 7);
 230     }
 231 }
 232
 233 BOOST_AUTO_TEST_CASE(chebyshev_transform_test)
 234 {
 235 #ifdef TEST1
 236     test_chebyshev_chebyshev_transform<float>();
 237     test_sin_chebyshev_transform<float>();
 238     test_atap_examples<float>();
 239     test_sinc_chebyshev_transform<float>();
 240 #endif
 241 #ifdef TEST2
 242     test_chebyshev_chebyshev_transform<double>();
 243     test_sin_chebyshev_transform<double>();
 244     test_atap_examples<double>();
 245     test_sinc_chebyshev_transform<double>();
 246 #endif
 247 #ifdef TEST3
 248     test_chebyshev_chebyshev_transform<long double>();
 249     test_sin_chebyshev_transform<long double>();
 250     test_atap_examples<long double>();
 251     test_sinc_chebyshev_transform<long double>();
 252 #endif
 253 #ifdef TEST4
 254 #ifdef BOOST_HAS_FLOAT128
 255     test_chebyshev_chebyshev_transform<__float128>();
 256     test_sin_chebyshev_transform<__float128>();
 257     test_atap_examples<__float128>();
 258     test_sinc_chebyshev_transform<__float128>();
 259 #endif
 260 #endif
 261 }
 262
 263