[/ Copyright (c) 2020 Nick Thompson 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) ] [section:quintic_hermite Quintic Hermite interpolation] [heading Synopsis] `` #include namespace boost::math::interpolators { template class quintic_hermite { public: using Real = typename RandomAccessContainer::value_type; quintic_hermite(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2) inline Real operator()(Real x) const; inline Real prime(Real x) const; inline Real double_prime(Real x) const; std::pair domain() const; friend std::ostream& operator<<(std::ostream & os, const quintic_hermite & m); void push_back(Real x, Real y, Real dydx, Real d2ydx2); }; template class cardinal_quintic_hermite { public: using Real = typename RandomAccessContainer::value_type; cardinal_quintic_hermite(RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2, Real x0, Real dx); inline Real operator()(Real x) const; inline Real prime(Real x) const; inline Real double_prime(Real x) const; std::pair domain() const; }; template class cardinal_quintic_hermite_aos { public: using Point = typename RandomAccessContainer::value_type; using Real = typename Point::value_type; cardinal_quintic_hermite_aos(RandomAccessContainer && data, Real x0, Real dx) inline Real operator()(Real x) const; inline Real prime(Real x) const; inline Real double_prime(Real x) const; std::pair domain() const; } `` [heading Quintic Hermite Interpolation] The quintic Hermite interpolator takes a list of possibly non-uniformly spaced abscissas, ordinates, and their velocities and accelerations which are used to construct a quintic interpolating polynomial between segments. This is useful for taking solution skeletons from ODE steppers and turning them into a continuous function, provided that the right-hand side /f/(/x/, /y/) is differentiable along the solution path. The interpolant is /C/[super 2] and its evaluation has [bigo](log(/N/)) complexity. An example usage is as follows: std::vector x{1,2,3, 4, 5, 6}; std::vector y{7,8,9,10,11,12}; std::vector dydx{1,1,1,1,1,1}; std::vector d2ydx2{0,0,0,0,0,0}; using boost::math::interpolators::quintic_hermite; auto spline = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)); // evaluate at a point: double z = spline(3.4); // evaluate derivative at a point: double zprime = spline.prime(3.4); Periodically, it is helpful to see what data the interpolator has. This can be achieved via std::cout << spline << "\n"; Note that the interpolator is pimpl'd, so that copying the class is cheap, and hence it can be shared between threads. (The call operator and `.prime()` are threadsafe.) The interpolator can be updated in constant time. Hence we can use `boost::circular_buffer` to do real-time interpolation. #include ... boost::circular_buffer initial_x{1,2,3,4}; boost::circular_buffer initial_y{4,5,6,7}; boost::circular_buffer initial_dydx{1,1,1,1}; boost::circular_buffer initial_d2ydx2{0,0,0,0}; auto circular_akima = quintic_hermite(std::move(initial_x), std::move(initial_y), std::move(initial_dydx), std::move(initial_d2ydx2)); // interpolate via call operation: double y = circular_akima(3.5); // add new data: circular_akima.push_back(5, 8, 1, 0); // interpolate at 4.5: y = circular_akima(4.5); [$../graphs/quintic_sine_approximation.svg] For equispaced data, we can use `cardinal_quintic_hermite` or `cardinal_quintic_hermite_aos` to get constant-time evaluation. This is useful in memory-constrained or performance critical applications where data is equispaced. [heading Complexity and Performance] The following google benchmark demonstrates the cost of the call operator for this interpolator: ``` Run on (16 X 4300 MHz CPU s) CPU Caches: L1 Data 32K (x8) L1 Instruction 32K (x8) L2 Unified 1024K (x8) L3 Unified 11264K (x1) Load Average: 0.92, 0.64, 0.35 -------------------------------------------------- Benchmark Time -------------------------------------------------- QuinticHermite/8 8.14 ns QuinticHermite/16 8.69 ns QuinticHermite/32 9.42 ns QuinticHermite/64 9.90 ns QuinticHermite/128 10.4 ns QuinticHermite/256 10.9 ns QuinticHermite/512 11.6 ns QuinticHermite/1024 12.3 ns QuinticHermite/2048 12.8 ns QuinticHermite/4096 13.6 ns QuinticHermite/8192 14.6 ns QuinticHermite/16384 15.5 ns QuinticHermite/32768 17.4 ns QuinticHermite/65536 18.5 ns QuinticHermite/131072 18.8 ns QuinticHermite/262144 19.8 ns QuinticHermite/524288 20.5 ns QuinticHermite/1048576 21.6 ns QuinticHermite/2097152 22.5 ns QuinticHermite/4194304 24.2 ns QuinticHermite/8388608 26.6 ns QuinticHermite/16777216 26.7 ns QuinticHermite_BigO 1.14 lgN CardinalQuinticHermite/256 5.22 ns CardinalQuinticHermite/512 5.21 ns CardinalQuinticHermite/1024 5.21 ns CardinalQuinticHermite/2048 5.32 ns CardinalQuinticHermite/4096 5.33 ns CardinalQuinticHermite/8192 5.50 ns CardinalQuinticHermite/16384 5.74 ns CardinalQuinticHermite/32768 7.74 ns CardinalQuinticHermite/65536 10.6 ns CardinalQuinticHermite/131072 10.7 ns CardinalQuinticHermite/262144 10.6 ns CardinalQuinticHermite/524288 10.5 ns CardinalQuinticHermite/1048576 10.6 ns CardinalQuinticHermite_BigO 7.57 (1) CardinalQuinticHermiteAOS/256 5.27 ns CardinalQuinticHermiteAOS/512 5.26 ns CardinalQuinticHermiteAOS/1024 5.26 ns CardinalQuinticHermiteAOS/2048 5.28 ns CardinalQuinticHermiteAOS/4096 5.30 ns CardinalQuinticHermiteAOS/8192 5.41 ns CardinalQuinticHermiteAOS/16384 5.89 ns CardinalQuinticHermiteAOS/32768 5.97 ns CardinalQuinticHermiteAOS/65536 5.96 ns CardinalQuinticHermiteAOS/131072 5.92 ns CardinalQuinticHermiteAOS/262144 5.94 ns CardinalQuinticHermiteAOS/524288 5.96 ns CardinalQuinticHermiteAOS/1048576 5.93 ns CardinalQuinticHermiteAOS_BigO 5.64 (1) ``` [endsect] [/section:quintic_hermite]