123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130 |
- /*
- * Copyright Nick Thompson, 2019
- * Use, modification and distribution are subject to the
- * Boost Software License, Version 1.0. (See accompanying file
- * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
- */
- #include "math_unit_test.hpp"
- #include <numeric>
- #include <utility>
- #include <boost/math/interpolators/cardinal_quadratic_b_spline.hpp>
- using boost::math::interpolators::cardinal_quadratic_b_spline;
- template<class Real>
- void test_constant()
- {
- Real c = 7.2;
- Real t0 = 0;
- Real h = Real(1)/Real(16);
- size_t n = 512;
- std::vector<Real> v(n, c);
- auto qbs = cardinal_quadratic_b_spline<Real>(v.data(), v.size(), t0, h);
- size_t i = 0;
- while (i < n) {
- Real t = t0 + i*h;
- CHECK_ULP_CLOSE(c, qbs(t), 2);
- CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 100*std::numeric_limits<Real>::epsilon());
- ++i;
- }
- i = 0;
- while (i < n) {
- Real t = t0 + i*h + h/2;
- CHECK_ULP_CLOSE(c, qbs(t), 2);
- CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 300*std::numeric_limits<Real>::epsilon());
- t = t0 + i*h + h/4;
- CHECK_ULP_CLOSE(c, qbs(t), 2);
- CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 150*std::numeric_limits<Real>::epsilon());
- ++i;
- }
- }
- template<class Real>
- void test_linear()
- {
- Real m = 8.3;
- Real b = 7.2;
- Real t0 = 0;
- Real h = Real(1)/Real(16);
- size_t n = 512;
- std::vector<Real> y(n);
- for (size_t i = 0; i < n; ++i) {
- Real t = i*h;
- y[i] = m*t + b;
- }
- auto qbs = cardinal_quadratic_b_spline<Real>(y.data(), y.size(), t0, h);
- size_t i = 0;
- while (i < n) {
- Real t = t0 + i*h;
- CHECK_ULP_CLOSE(m*t+b, qbs(t), 2);
- CHECK_ULP_CLOSE(m, qbs.prime(t), 820);
- ++i;
- }
- i = 0;
- while (i < n) {
- Real t = t0 + i*h + h/2;
- CHECK_ULP_CLOSE(m*t+b, qbs(t), 2);
- CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits<Real>::epsilon());
- t = t0 + i*h + h/4;
- CHECK_ULP_CLOSE(m*t+b, qbs(t), 3);
- CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits<Real>::epsilon());
- ++i;
- }
- }
- template<class Real>
- void test_quadratic()
- {
- Real a = 8.2;
- Real b = 7.2;
- Real c = -9.2;
- Real t0 = 0;
- Real h = Real(1)/Real(16);
- size_t n = 513;
- std::vector<Real> y(n);
- for (size_t i = 0; i < n; ++i) {
- Real t = i*h;
- y[i] = a*t*t + b*t + c;
- }
- Real t_max = t0 + (n-1)*h;
- auto qbs = cardinal_quadratic_b_spline<Real>(y, t0, h, b, 2*a*t_max + b);
- size_t i = 0;
- while (i < n) {
- Real t = t0 + i*h;
- CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 2);
- ++i;
- }
- i = 0;
- while (i < n) {
- Real t = t0 + i*h + h/2;
- CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 47);
- t = t0 + i*h + h/4;
- if (!CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 104)) {
- std::cerr << " Problem abscissa t = " << t << "\n";
- }
- ++i;
- }
- }
- int main()
- {
- test_constant<float>();
- test_constant<double>();
- test_constant<long double>();
- test_linear<float>();
- test_linear<double>();
- test_linear<long double>();
- test_quadratic<double>();
- test_quadratic<long double>();
- return boost::math::test::report_errors();
- }
|