1 // Copyright Nick Thompson, 2019
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0.
4 // (See accompanying file LICENSE_1_0.txt
5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
7 #ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_QUADRATIC_B_SPLINE_DETAIL_HPP
8 #define BOOST_MATH_INTERPOLATORS_CARDINAL_QUADRATIC_B_SPLINE_DETAIL_HPP
13 namespace boost{ namespace math{ namespace interpolators{ namespace detail{
16 Real b2_spline(Real x) {
21 Real y = absx - 1/Real(2);
22 Real z = absx + 1/Real(2);
25 if (absx < Real(3)/Real(2))
27 Real y = absx - Real(3)/Real(2);
30 return static_cast<Real>(0);
34 Real b2_spline_prime(Real x) {
36 return -b2_spline_prime(-x);
43 if (x < Real(3)/Real(2))
45 return x - Real(3)/Real(2);
47 return static_cast<Real>(0);
52 class cardinal_quadratic_b_spline_detail
55 // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
56 // y[0] = y(a), y[n -1] = y(b), step_size = (b - a)/(n -1).
58 cardinal_quadratic_b_spline_detail(const Real* const y,
60 Real t0 /* initial time, left endpoint */,
61 Real h /*spacing, stepsize*/,
62 Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
63 Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN())
66 throw std::logic_error("Spacing must be > 0.");
72 throw std::logic_error("The interpolator requires at least 3 points.");
77 if (isnan(left_endpoint_derivative)) {
78 // http://web.media.mit.edu/~crtaylor/calculator.html
79 a = -3*y[0] + 4*y[1] - y[2];
82 a = 2*h*left_endpoint_derivative;
86 if (isnan(right_endpoint_derivative)) {
87 b = 3*y[n-1] - 4*y[n-2] + y[n-3];
90 b = 2*h*right_endpoint_derivative;
93 m_alpha.resize(n + 2);
95 // Begin row reduction:
96 std::vector<Real> rhs(n + 2, std::numeric_limits<Real>::quiet_NaN());
97 std::vector<Real> super_diagonal(n + 2, std::numeric_limits<Real>::quiet_NaN());
100 rhs[rhs.size() - 1] = b;
102 super_diagonal[0] = 0;
104 for(size_t i = 1; i < rhs.size() - 1; ++i) {
106 super_diagonal[i] = 1;
109 // Patch up 5-diagonal problem:
110 rhs[1] = (rhs[1] - rhs[0])/6;
111 super_diagonal[1] = Real(1)/Real(3);
112 // First two rows are now:
114 // 0 1 1/3| (8y0+2hy0')/6
117 // Start traditional tridiagonal row reduction:
118 for (size_t i = 2; i < rhs.size() - 1; ++i) {
119 Real diagonal = 6 - super_diagonal[i - 1];
120 rhs[i] = (rhs[i] - rhs[i - 1])/diagonal;
121 super_diagonal[i] /= diagonal;
124 // 1 sd[n-1] 0 | rhs[n-1]
125 // 0 1 sd[n] | rhs[n]
128 rhs[n+1] = rhs[n+1] + rhs[n-1];
129 Real bottom_subdiagonal = super_diagonal[n-1];
132 // 1 sd[n-1] 0 | rhs[n-1]
133 // 0 1 sd[n] | rhs[n]
136 rhs[n+1] = (rhs[n+1]-bottom_subdiagonal*rhs[n])/(1-bottom_subdiagonal*super_diagonal[n]);
138 m_alpha[n+1] = rhs[n+1];
139 for (size_t i = n; i > 0; --i) {
140 m_alpha[i] = rhs[i] - m_alpha[i+1]*super_diagonal[i];
142 m_alpha[0] = m_alpha[2] + rhs[0];
145 Real operator()(Real t) const {
146 if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) {
147 const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work.";
148 throw std::domain_error(err_msg);
150 // Let k, gamma be defined via t = t0 + kh + gamma * h.
151 // Now find all j: |k-j+1+gamma|< 3/2, or, in other words
152 // j_min = ceil((t-t0)/h - 1/2)
153 // j_max = floor(t-t0)/h + 5/2)
156 Real x = (t-m_t0)*m_inv_h;
157 size_t j_min = ceil(x - Real(1)/Real(2));
158 size_t j_max = ceil(x + Real(5)/Real(2));
159 if (j_max >= m_alpha.size()) {
160 j_max = m_alpha.size() - 1;
165 for (size_t j = j_min; j <= j_max; ++j) {
166 y += m_alpha[j]*detail::b2_spline(x - j);
171 Real prime(Real t) const {
172 if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) {
173 const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work.";
174 throw std::domain_error(err_msg);
176 // Let k, gamma be defined via t = t0 + kh + gamma * h.
177 // Now find all j: |k-j+1+gamma|< 3/2, or, in other words
178 // j_min = ceil((t-t0)/h - 1/2)
179 // j_max = floor(t-t0)/h + 5/2)
182 Real x = (t-m_t0)*m_inv_h;
183 size_t j_min = ceil(x - Real(1)/Real(2));
184 size_t j_max = ceil(x + Real(5)/Real(2));
185 if (j_max >= m_alpha.size()) {
186 j_max = m_alpha.size() - 1;
191 for (size_t j = j_min; j <= j_max; ++j) {
192 y += m_alpha[j]*detail::b2_spline_prime(x - j);
198 return m_t0 + (m_alpha.size()-3)/m_inv_h;
202 std::vector<Real> m_alpha;