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

Second derivative for septic Hermite spline.

Browse Source
This commit is contained in:
NAThompson
2020-02-24 07:06:21 -05:00
parent 8343534042
commit 38c2dd376b
3 changed files with 118 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files

108
include/boost/math/interpolators/detail/septic_hermite_detail.hpp
View File

@@ -336,9 +336,57 @@ public:
return dydx;
}
inline Real double_prime(Real x) const
inline Real double_prime(Real x) constq
{
return std::numeric_limits<Real>::quiet_NaN();
Real xf = x0_ + (y_.size()-1)/inv_dx_;
if (x < x0_ || x > xf)
{
std::ostringstream oss;
oss.precision(std::numeric_limits<Real>::digits10+3);
oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
<< x0_ << ", " << xf << "]";
throw std::domain_error(oss.str());
}
if (x == xf)
{
return d2y_.back()*2*inv_dx_*inv_dx_;
}
return this->unchecked_double_prime(x);
}
inline Real unchecked_double_prime(Real x) const
{
using std::floor;
Real s3 = (x-x0_)*inv_dx_;
Real ii = floor(s3);
auto i = static_cast<decltype(y_.size())>(ii);
Real t = s3 - ii;
Real y0 = y_[i];
Real y1 = y_[i+1];
Real dy0 = dy_[i];
Real dy1 = dy_[i+1];
Real a0 = d2y_[i];
Real a1 = d2y_[i+1];
Real j0 = d3y_[i];
Real j1 = d3y_[i+1];
Real t2 = t*t;
Real z0 = 420*t2*(1 + t*(-4 + t*(5 - 2*t)));
Real z1 = 60*t2*(-4 + t*(15 + t*(-18 + 7*t)));
Real z2 = 60*t2*(-3 + t*(13 + t*(-17 + 7*t)));
Real z3 = (1 + t2*(-60 + t*(200 + t*(-225 + 84*t))));
Real z4 = t2*(30 + t*(-140 + t*(195 - 84*t)));
Real z5 = t*(1 + t*(-8 + t*(20 + t*(-20 + 7*t))));
Real z6 = t2*(-2 + t*(10 + t*(-15 + 7*t)));
Real d2ydx2 = z0*(y1-y0)*inv_dx_*inv_dx_;
d2ydx2 += (z1*dy0 + z2*dy1)*inv_dx_*inv_dx_;
d2ydx2 += (z3*a0 + z4*a1)*2*inv_dx_*inv_dx_;
d2ydx2 += 6*(z5*j0 + z6*j1)/(inv_dx_*inv_dx_);
return d2ydx2;
}
private:
@@ -437,7 +485,7 @@ public:
}
if (x == xf)
{
return data_.back()[1];
return data_.back()[1]*inv_dx_;
}
return this->unchecked_prime(x);
@@ -449,7 +497,7 @@ public:
Real ii = floor(s3);
auto i = static_cast<decltype(data_.size())>(ii);
Real t = s3 - ii;
Real y0 = data_[i][0];
Real y1 = data_[i+1][0];
Real dy0 = data_[i][1];
@@ -475,9 +523,57 @@ public:
return dydx;
}
inline Real double_prime(Real x) const
inline Real double_prime(Real x) const
{
return std::numeric_limits<Real>::quiet_NaN();
Real xf = x0_ + (data_.size()-1)/inv_dx_;
if (x < x0_ || x > xf)
{
std::ostringstream oss;
oss.precision(std::numeric_limits<Real>::digits10+3);
oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
<< x0_ << ", " << xf << "]";
throw std::domain_error(oss.str());
}
if (x == xf)
{
return data_.back()[2]*2*inv_dx_*inv_dx_;
}q
return this->unchecked_double_prime(x);
}
inline Real unchecked_double_prime(Real x) const
{
using std::floor;
Real s3 = (x-x0_)*inv_dx_;
Real ii = floor(s3);
auto i = static_cast<decltype(data_.size())>(ii);
Real t = s3 - ii;
Real y0 = data_[i][0];
Real y1 = data_[i+1][0];
Real dy0 = data_[i][1];
Real dy1 = data_[i+1][1];
Real a0 = data_[i][2];
Real a1 = data_[i+1][2];
Real j0 = data_[i][3];
Real j1 = data_[i+1][3];
Real t2 = t*t;
Real z0 = 420*t2*(1 + t*(-4 + t*(5 - 2*t)));
Real z1 = 60*t2*(-4 + t*(15 + t*(-18 + 7*t)));
Real z2 = 60*t2*(-3 + t*(13 + t*(-17 + 7*t)));
Real z3 = (1 + t2*(-60 + t*(200 + t*(-225 + 84*t))));
Real z4 = t2*(30 + t*(-140 + t*(195 - 84*t)));
Real z5 = t*(1 + t*(-8 + t*(20 + t*(-20 + 7*t))));
Real z6 = t2*(-2 + t*(10 + t*(-15 + 7*t)));
Real d2ydx2 = z0*(y1-y0)*inv_dx_*inv_dx_;
d2ydx2 += (z1*dy0 + z2*dy1)*inv_dx_*inv_dx_;
d2ydx2 += (z3*a0 + z4*a1)*2*inv_dx_*inv_dx_;
d2ydx2 += 6*(z5*j0 + z6*j1)/(inv_dx_*inv_dx_);
return d2ydx2;
}
private:

4
include/boost/math/special_functions/daubechies_scaling.hpp
View File

@@ -491,6 +491,10 @@ public:
{
return m_qh->unchecked_double_prime(x);
}
if constexpr (p >= 10)
{
return m_sh->unchecked_double_prime(x);
}
}
std::pair<Real, Real> support() const

13
test/septic_hermite_test.cpp
View File

@@ -51,9 +51,11 @@ void test_constant()
d2ydx2.resize(128, 0);
d3ydx3.resize(128, 0);
auto csh = cardinal_septic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3), x0, dx);
for (Real t = x0; t <= 127; t += 0.25) {
for (Real t = x0; t <= 127; t += 0.25)
{
CHECK_ULP_CLOSE(Real(7), csh(t), 24);
CHECK_ULP_CLOSE(Real(0), csh.prime(t), 24);
CHECK_ULP_CLOSE(Real(0), csh.double_prime(t), 24);
}
std::vector<std::array<Real, 4>> data(128);
@@ -69,6 +71,7 @@ void test_constant()
{
CHECK_ULP_CLOSE(Real(7), csh_aos(t), 24);
CHECK_ULP_CLOSE(Real(0), csh_aos.prime(t), 24);
CHECK_ULP_CLOSE(Real(0), csh_aos.double_prime(t), 24);
}
}
@@ -130,6 +133,7 @@ void test_linear()
{
CHECK_ULP_CLOSE(t, csh(t), 15);
CHECK_ULP_CLOSE(Real(1), csh.prime(t), 15);
CHECK_ULP_CLOSE(Real(0), csh.double_prime(t), 15);
}
std::vector<std::array<Real, 4>> data(10);
@@ -145,6 +149,7 @@ void test_linear()
{
CHECK_ULP_CLOSE(t, csh_aos(t), 15);
CHECK_ULP_CLOSE(Real(1), csh_aos.prime(t), 15);
CHECK_ULP_CLOSE(Real(0), csh_aos.double_prime(t), 15);
}
}
@@ -231,6 +236,7 @@ void test_quadratic()
{
CHECK_ULP_CLOSE(t*t/2, csh(t), 24);
CHECK_ULP_CLOSE(t, csh.prime(t), 24);
CHECK_ULP_CLOSE(Real(1), csh.double_prime(t), 24);
}
std::vector<std::array<Real, 4>> data(10);
@@ -246,6 +252,7 @@ void test_quadratic()
{
CHECK_ULP_CLOSE(t*t/2, csh_aos(t), 24);
CHECK_ULP_CLOSE(t, csh_aos.prime(t), 24);
CHECK_ULP_CLOSE(Real(1), csh_aos.double_prime(t), 24);
}
}
@@ -303,6 +310,7 @@ void test_cubic()
{
CHECK_ULP_CLOSE(t*t*t, csh(t), 151);
CHECK_ULP_CLOSE(3*t*t, csh.prime(t), 151);
CHECK_ULP_CLOSE(6*t, csh.double_prime(t), 151);
}
std::vector<std::array<Real, 4>> data(8);
@@ -319,6 +327,7 @@ void test_cubic()
{
CHECK_ULP_CLOSE(t*t*t, csh_aos(t), 151);
CHECK_ULP_CLOSE(3*t*t, csh_aos.prime(t), 151);
CHECK_ULP_CLOSE(6*t, csh_aos.double_prime(t), 151);
}
}
@@ -377,6 +386,7 @@ void test_quartic()
{
CHECK_ULP_CLOSE(t*t*t*t, csh(t), 117);
CHECK_ULP_CLOSE(4*t*t*t, csh.prime(t), 117);
CHECK_ULP_CLOSE(12*t*t, csh.double_prime(t), 117);
}
std::vector<std::array<Real, 4>> data(10);
@@ -393,6 +403,7 @@ void test_quartic()
{
CHECK_ULP_CLOSE(t*t*t*t, csh_aos(t), 117);
CHECK_ULP_CLOSE(4*t*t*t, csh_aos.prime(t), 117);
CHECK_ULP_CLOSE(12*t*t, csh_aos.double_prime(t), 117);
}
}