boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Cubic B-spline second derivatives

Browse Source
This commit is contained in:
NAThompson
2019-08-07 14:10:08 -04:00
parent 44e3ed82fa
commit edfb80e76d
3 changed files with 59 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
include/boost/math/interpolators/cubic_b_spline.hpp
View File

@@ -45,6 +45,8 @@ public:
Real prime(Real x) const;
Real double_prime(Real x) const;
private:
std::shared_ptr<detail::cubic_b_spline_imp<Real>> m_imp;
};
@@ -74,5 +76,12 @@ Real cubic_b_spline<Real>::prime(Real x) const
return m_imp->prime(x);
}
template<class Real>
Real cubic_b_spline<Real>::double_prime(Real x) const
{
return m_imp->double_prime(x);
}
}}
#endif

43
include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp
View File

@@ -32,6 +32,8 @@ public:
Real prime(Real x) const;
Real double_prime(Real x) const;
private:
std::vector<Real> m_beta;
Real m_h_inv;
@@ -79,6 +81,25 @@ Real b3_spline_prime(Real x)
return (Real) 0;
}
template<class Real>
Real b3_spline_double_prime(Real x)
{
if (x < 0)
{
return b3_spline_double_prime(-x);
}
if (x < 1)
{
return (3*boost::math::constants::half<Real>()*x - 2) + x*(3*boost::math::constants::half<Real>());
}
if (x < 2)
{
return (2 - x);
}
return (Real) 0;
}
template <class Real>
template <class BidiIterator>
@@ -283,5 +304,27 @@ Real cubic_b_spline_imp<Real>::prime(Real x) const
return z*m_h_inv;
}
template<class Real>
Real cubic_b_spline_imp<Real>::double_prime(Real x) const
{
Real z = 0;
Real t = m_h_inv*(x - m_a) + 1;
using std::max;
using std::min;
using std::ceil;
using std::floor;
size_t k_min = (size_t) (max)(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2)));
size_t k_max = (size_t) (min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2)));
for (size_t k = k_min; k <= k_max; ++k)
{
z += m_beta[k]*b3_spline_double_prime(t - k);
}
return z*m_h_inv*m_h_inv;
}
}}}
#endif

7
test/test_cubic_b_spline.cpp
View File

@@ -30,6 +30,9 @@ void test_b3_spline()
BOOST_CHECK_SMALL(boost::math::detail::b3_spline<Real>(-2.5), (Real) 0);
BOOST_CHECK_SMALL(boost::math::detail::b3_spline_prime<Real>(2.5), (Real) 0);
BOOST_CHECK_SMALL(boost::math::detail::b3_spline_prime<Real>(-2.5), (Real) 0);
BOOST_CHECK_SMALL(boost::math::detail::b3_spline_double_prime<Real>(2.5), (Real) 0);
BOOST_CHECK_SMALL(boost::math::detail::b3_spline_double_prime<Real>(-2.5), (Real) 0);
// On the boundary of support:
BOOST_CHECK_SMALL(boost::math::detail::b3_spline<Real>(2), (Real) 0);
@@ -52,6 +55,7 @@ void test_b3_spline()
Real arg = i*0.01;
BOOST_CHECK_CLOSE(boost::math::detail::b3_spline<Real>(arg), boost::math::detail::b3_spline<Real>(arg), eps);
BOOST_CHECK_CLOSE(boost::math::detail::b3_spline_prime<Real>(-arg), -boost::math::detail::b3_spline_prime<Real>(arg), eps);
BOOST_CHECK_CLOSE(boost::math::detail::b3_spline_double_prime<Real>(-arg), boost::math::detail::b3_spline_double_prime<Real>(arg), eps);
}
}
@@ -109,6 +113,9 @@ void test_constant_function()
BOOST_CHECK_CLOSE(y, constant, 10*std::numeric_limits<Real>::epsilon());
Real y_prime = spline.prime(i*step + a + 0.002);
BOOST_CHECK_SMALL(y_prime, 5000*std::numeric_limits<Real>::epsilon());
Real y_double_prime = spline.double_prime(i*step + a + 0.002);
BOOST_CHECK_SMALL(y_double_prime, 5000*std::numeric_limits<Real>::epsilon());
}
// Test that correctly specified left and right-derivatives work properly: