1[/ 2Copyright (c) 2017 Nick Thompson 3Use, modification and distribution are subject to the 4Boost Software License, Version 1.0. (See accompanying file 5LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 6] 7 8[section:cardinal_cubic_b Cardinal Cubic B-spline interpolation] 9 10[heading Synopsis] 11`` 12 #include <boost/math/interpolators/cardinal_cubic_b_spline.hpp> 13`` 14 15 namespace boost { namespace math { namespace interpolators { 16 17 template <class Real> 18 class cardinal_cubic_b_spline 19 { 20 public: 21 22 template <class BidiIterator> 23 cardinal_cubic_b_spline(BidiIterator a, BidiIterator b, Real left_endpoint, Real step_size, 24 Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(), 25 Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN()); 26 cardinal_cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size, 27 Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(), 28 Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN()); 29 30 Real operator()(Real x) const; 31 32 Real prime(Real x) const; 33 34 Real double_prime(Real x) const; 35 }; 36 37 }}} // namespaces 38 39 40[heading Cardinal Cubic B-Spline Interpolation] 41 42The cardinal cubic /B/-spline class provided by Boost allows fast and accurate interpolation of a function which is known at equally spaced points. 43The cubic /B/-spline interpolation is numerically stable as it uses compactly supported basis functions constructed via iterative convolution. 44This is to be contrasted to one-sided power function cubic spline interpolation which is ill-conditioned as the global support of cubic polynomials causes small changes far from the evaluation point to exert a large influence on the calculated value. 45 46There are many use cases for interpolating a function at equally spaced points. 47One particularly important example is solving ODEs whose coefficients depend on data determined from experiment or numerical simulation. 48Since most ODE steppers are adaptive, they must be able to sample the coefficients at arbitrary points; 49not just at the points we know the values of our function. 50 51The first two arguments to the constructor are either: 52 53* A pair of bidirectional iterators into the data, or 54* A pointer to the data, and a length of the data array. 55 56These are then followed by: 57 58* The start of the functions domain, 59* The step size. 60 61Optionally, you may provide two additional arguments to the constructor, namely the derivative of the function at the left endpoint, and the derivative at the right endpoint. 62If you do not provide these arguments, they will be estimated using one-sided finite-difference formulas. 63An example of a valid call to the constructor is 64 65 std::vector<double> f{0.01, -0.02, 0.3, 0.8, 1.9, -8.78, -22.6}; 66 double t0 = 0; 67 double h = 0.01; 68 boost::math::interpolators::cardinal_cubic_b_spline<double> spline(f.begin(), f.end(), t0, h); 69 70The endpoints are estimated using a one-sided finite-difference formula. 71If you know the derivative at the endpoint, you may pass it to the constructor via 72 73 boost::math::interpolators::cardinal_cubic_b_spline<double> spline(f.begin(), f.end(), t0, h, a_prime, b_prime); 74 75To evaluate the interpolant at a point, we simply use 76 77 double y = spline(x); 78 79and to evaluate the derivative of the interpolant we use 80 81 double yp = spline.prime(x); 82 83Be aware that the accuracy guarantees on the derivative of the spline are an order lower than the guarantees on the original function, 84see [@http://www.springer.com/us/book/9780387984087 Numerical Analysis, Graduate Texts in Mathematics, 181, Rainer Kress] for details. 85 86The last interesting member is the second derivative, evaluated via 87 88 double ypp = spline.double_prime(x); 89 90The basis functions of the spline are cubic polynomials, so the second derivative is simply linear interpolation. 91But the interpolation is not constrained by data (you don't pass in the second derivatives), 92and hence is less accurate than would be naively expected from a linear interpolator. 93The problem is especially pronounced at the boundaries, where the second derivative is totally unconstrained. 94Use the second derivative of the cubic B-spline interpolator only in desperation. 95The quintic /B/-spline interpolator is recommended for cases where second derivatives are needed. 96 97 98[heading Complexity and Performance] 99 100The call to the constructor requires [bigo](/n/) operations, where /n/ is the number of points to interpolate. 101Each call the the interpolant is [bigo](1) (constant time). 102On the author's Intel Xeon E3-1230, this takes 21ns as long as the vector is small enough to fit in cache. 103 104[heading Accuracy] 105 106Let /h/ be the stepsize. If /f/ is four-times continuously differentiable, then the interpolant is ['[bigo](h[super 4])] accurate and the derivative is ['[bigo](h[super 3])] accurate. 107 108[heading Testing] 109 110Since the interpolant obeys [role serif_italic s(x[sub j]) = f(x[sub j])] at all interpolation points, 111the tests generate random data and evaluate the interpolant at the interpolation points, 112validating that equality with the data holds. 113 114In addition, constant, linear, and quadratic functions are interpolated to ensure that the interpolant behaves as expected. 115 116[heading Example] 117 118[import ../../example/cardinal_cubic_b_spline_example.cpp] 119 120[cubic_b_spline_example] 121 122[cubic_b_spline_example_out] 123 124[endsect] [/section:cardinal_cubic_b] 125