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
11 namespace boost{ namespace math{ namespace interpolators{ namespace detail{
14 Real b2_spline(Real x) {
19 Real y = absx - 1/Real(2);
20 Real z = absx + 1/Real(2);
23 if (absx < Real(3)/Real(2))
25 Real y = absx - Real(3)/Real(2);
32 Real b2_spline_prime(Real x) {
34 return -b2_spline_prime(-x);
41 if (x < Real(3)/Real(2))
43 return x - Real(3)/Real(2);
50 class cardinal_quadratic_b_spline_detail
53 // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
54 // y[0] = y(a), y[n -1] = y(b), step_size = (b - a)/(n -1).
56 cardinal_quadratic_b_spline_detail(const Real* const y,
58 Real t0 /* initial time, left endpoint */,
59 Real h /*spacing, stepsize*/,
60 Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
61 Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN())
64 throw std::logic_error("Spacing must be > 0.");
70 throw std::logic_error("The interpolator requires at least 3 points.");
75 if (isnan(left_endpoint_derivative)) {
76 // http://web.media.mit.edu/~crtaylor/calculator.html
77 a = -3*y[0] + 4*y[1] - y[2];
80 a = 2*h*left_endpoint_derivative;
84 if (isnan(right_endpoint_derivative)) {
85 b = 3*y[n-1] - 4*y[n-2] + y[n-3];
88 b = 2*h*right_endpoint_derivative;
91 m_alpha.resize(n + 2);
93 // Begin row reduction:
94 std::vector<Real> rhs(n + 2, std::numeric_limits<Real>::quiet_NaN());
95 std::vector<Real> super_diagonal(n + 2, std::numeric_limits<Real>::quiet_NaN());
98 rhs[rhs.size() - 1] = b;
100 super_diagonal[0] = 0;
102 for(size_t i = 1; i < rhs.size() - 1; ++i) {
104 super_diagonal[i] = 1;
107 // Patch up 5-diagonal problem:
108 rhs[1] = (rhs[1] - rhs[0])/6;
109 super_diagonal[1] = Real(1)/Real(3);
110 // First two rows are now:
112 // 0 1 1/3| (8y0+2hy0')/6
115 // Start traditional tridiagonal row reduction:
116 for (size_t i = 2; i < rhs.size() - 1; ++i) {
117 Real diagonal = 6 - super_diagonal[i - 1];
118 rhs[i] = (rhs[i] - rhs[i - 1])/diagonal;
119 super_diagonal[i] /= diagonal;
122 // 1 sd[n-1] 0 | rhs[n-1]
123 // 0 1 sd[n] | rhs[n]
126 rhs[n+1] = rhs[n+1] + rhs[n-1];
127 Real bottom_subdiagonal = super_diagonal[n-1];
130 // 1 sd[n-1] 0 | rhs[n-1]
131 // 0 1 sd[n] | rhs[n]
134 rhs[n+1] = (rhs[n+1]-bottom_subdiagonal*rhs[n])/(1-bottom_subdiagonal*super_diagonal[n]);
136 m_alpha[n+1] = rhs[n+1];
137 for (size_t i = n; i > 0; --i) {
138 m_alpha[i] = rhs[i] - m_alpha[i+1]*super_diagonal[i];
140 m_alpha[0] = m_alpha[2] + rhs[0];
143 Real operator()(Real t) const {
144 if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) {
145 const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work.";
146 throw std::domain_error(err_msg);
148 // Let k, gamma be defined via t = t0 + kh + gamma * h.
149 // Now find all j: |k-j+1+gamma|< 3/2, or, in other words
150 // j_min = ceil((t-t0)/h - 1/2)
151 // j_max = floor(t-t0)/h + 5/2)
154 Real x = (t-m_t0)*m_inv_h;
155 size_t j_min = ceil(x - Real(1)/Real(2));
156 size_t j_max = ceil(x + Real(5)/Real(2));
157 if (j_max >= m_alpha.size()) {
158 j_max = m_alpha.size() - 1;
163 for (size_t j = j_min; j <= j_max; ++j) {
164 y += m_alpha[j]*detail::b2_spline(x - j);
169 Real prime(Real t) const {
170 if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) {
171 const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work.";
172 throw std::domain_error(err_msg);
174 // Let k, gamma be defined via t = t0 + kh + gamma * h.
175 // Now find all j: |k-j+1+gamma|< 3/2, or, in other words
176 // j_min = ceil((t-t0)/h - 1/2)
177 // j_max = floor(t-t0)/h + 5/2)
180 Real x = (t-m_t0)*m_inv_h;
181 size_t j_min = ceil(x - Real(1)/Real(2));
182 size_t j_max = ceil(x + Real(5)/Real(2));
183 if (j_max >= m_alpha.size()) {
184 j_max = m_alpha.size() - 1;
189 for (size_t j = j_min; j <= j_max; ++j) {
190 y += m_alpha[j]*detail::b2_spline_prime(x - j);
196 return m_t0 + (m_alpha.size()-3)/m_inv_h;
200 std::vector<Real> m_alpha;