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

Accurate septic Hermite derivatives.

Browse Source
This commit is contained in:
NAThompson
2020-02-22 14:42:19 -05:00
parent 7755265029
commit 6edcfd6b78
2 changed files with 97 additions and 27 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -303,18 +303,37 @@ public:
return this->unchecked_prime(x);
}
inline Real unchecked_prime(Real x) const {
//TODO: Get the high accuracy approximation by differentiating the interpolant!
inline Real unchecked_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 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 t3 = t2*t;
Real z0 = 140*t3*(1 + t*(-3 + t*(3 - t)));
Real z1 = 1 + t3*(-80 + t*(225 + t*(-216 + 70*t)));
Real z2 = t3*(-60 + t*(195 + t*(-204 + 70*t)));
Real z3 = 1 + t2*(-20 + t*(50 + t*(-45 + 14*t)));
Real z4 = t2*(10 + t*(-35 + t*(39 - 14*t)));
Real z5 = 3 + t*(-16 + t*(30 + t*(-24 + 7*t)));
Real z6 = t*(-4 + t*(15 + t*(-18 + 7*t)));
// Velocity:
Real v0 = dy_[i];
Real v1 = dy_[i+1];
return (v0+v1)*inv_dx_/2;
Real dydx = z0*(y1-y0)*inv_dx_;
dydx += (z1*dy0 + z2*dy1)*inv_dx_;
dydx += 2*t*(z3*a0 + z4*a1)*inv_dx_;
dydx += t*t*(z5*j0 + z6*j1);
return dydx;
}
inline Real double_prime(Real x) const
@@ -425,7 +444,35 @@ public:
}
inline Real unchecked_prime(Real x) const {
return std::numeric_limits<Real>::quiet_NaN();
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 t3 = t2*t;
Real z0 = 140*t3*(1 + t*(-3 + t*(3 - t)));
Real z1 = 1 + t3*(-80 + t*(225 + t*(-216 + 70*t)));
Real z2 = t3*(-60 + t*(195 + t*(-204 + 70*t)));
Real z3 = 1 + t2*(-20 + t*(50 + t*(-45 + 14*t)));
Real z4 = t2*(10 + t*(-35 + t*(39 - 14*t)));
Real z5 = 3 + t*(-16 + t*(30 + t*(-24 + 7*t)));
Real z6 = t*(-4 + t*(15 + t*(-18 + 7*t)));
Real dydx = z0*(y1-y0)*inv_dx_;
dydx += (z1*dy0 + z2*dy1)*inv_dx_;
dydx += 2*t*(z3*a0 + z4*a1)*inv_dx_;
dydx += t*t*(z5*j0 + z6*j1);
return dydx;
}
inline Real double_prime(Real x) const

61
test/septic_hermite_test.cpp
View File

@@ -53,7 +53,7 @@ void test_constant()
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) {
CHECK_ULP_CLOSE(Real(7), csh(t), 24);
//CHECK_ULP_CLOSE(Real(0), csh.prime(t), 24);
CHECK_ULP_CLOSE(Real(0), csh.prime(t), 24);
}
std::vector<std::array<Real, 4>> data(128);
@@ -65,9 +65,10 @@ void test_constant()
data[i][3] = 0;
}
auto csh_aos = cardinal_septic_hermite_aos(std::move(data), 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_aos(t), 24);
//CHECK_ULP_CLOSE(Real(0), csh.prime(t), 24);
CHECK_ULP_CLOSE(Real(0), csh_aos.prime(t), 24);
}
}
@@ -84,7 +85,8 @@ void test_linear()
auto sh = septic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
for (Real t = 0; t <= 9; t += 0.25) {
for (Real t = 0; t <= 9; t += 0.25)
{
CHECK_ULP_CLOSE(Real(t), sh(t), 2);
CHECK_ULP_CLOSE(Real(1), sh.prime(t), 2);
}
@@ -94,7 +96,8 @@ void test_linear()
x.resize(512);
x[0] = dis(rng);
Real xmin = x[0];
for (size_t i = 1; i < x.size(); ++i) {
for (size_t i = 1; i < x.size(); ++i)
{
x[i] = x[i-1] + dis(rng);
}
Real xmax = x.back();
@@ -118,13 +121,15 @@ void test_linear()
dydx.resize(10, 1);
d2ydx2.resize(10, 0);
d3ydx3.resize(10, 0);
for (size_t i = 0; i < y.size(); ++i) {
for (size_t i = 0; i < y.size(); ++i)
{
y[i] = i;
}
auto csh = cardinal_septic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3), x0, dx);
for (Real t = 0; t <= 9; t += 0.125)
{
CHECK_ULP_CLOSE(t, csh(t), 15);
CHECK_ULP_CLOSE(Real(1), csh.prime(t), 15);
}
std::vector<std::array<Real, 4>> data(10);
@@ -139,6 +144,7 @@ void test_linear()
for (Real t = 0; t <= 9; t += 0.125)
{
CHECK_ULP_CLOSE(t, csh_aos(t), 15);
CHECK_ULP_CLOSE(Real(1), csh_aos.prime(t), 15);
}
}
@@ -224,6 +230,7 @@ void test_quadratic()
for (Real t = x0; t <= 9; t += 0.125)
{
CHECK_ULP_CLOSE(t*t/2, csh(t), 24);
CHECK_ULP_CLOSE(t, csh.prime(t), 24);
}
std::vector<std::array<Real, 4>> data(10);
@@ -238,6 +245,7 @@ void test_quadratic()
for (Real t = x0; t <= 9; t += 0.125)
{
CHECK_ULP_CLOSE(t*t/2, csh_aos(t), 24);
CHECK_ULP_CLOSE(t, csh_aos.prime(t), 24);
}
}
@@ -256,12 +264,14 @@ void test_cubic()
}
std::vector<Real> dydx(x.size());
for (size_t i = 0; i < y.size(); ++i) {
for (size_t i = 0; i < y.size(); ++i)
{
dydx[i] = 3*x[i]*x[i];
}
std::vector<Real> d2ydx2(x.size());
for (size_t i = 0; i < y.size(); ++i) {
for (size_t i = 0; i < y.size(); ++i)
{
d2ydx2[i] = 6*x[i];
}
std::vector<Real> d3ydx3(x.size(), 6);
@@ -280,7 +290,8 @@ void test_cubic()
dydx.resize(8);
d2ydx2.resize(8);
d3ydx3.resize(8,6);
for (size_t i = 0; i < y.size(); ++i) {
for (size_t i = 0; i < y.size(); ++i)
{
y[i] = i*i*i;
dydx[i] = 3*i*i;
d2ydx2[i] = 6*i;
@@ -288,8 +299,10 @@ void test_cubic()
auto csh = cardinal_septic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3), x0, dx);
for (Real t = 0; t <= xmax; t += 0.0078125) {
for (Real t = 0; t <= xmax; t += 0.0078125)
{
CHECK_ULP_CLOSE(t*t*t, csh(t), 151);
CHECK_ULP_CLOSE(3*t*t, csh.prime(t), 151);
}
std::vector<std::array<Real, 4>> data(8);
@@ -302,8 +315,10 @@ void test_cubic()
auto csh_aos = cardinal_septic_hermite_aos(std::move(data), x0, dx);
for (Real t = 0; t <= xmax; t += 0.0078125) {
for (Real t = 0; t <= xmax; t += 0.0078125)
{
CHECK_ULP_CLOSE(t*t*t, csh_aos(t), 151);
CHECK_ULP_CLOSE(3*t*t, csh_aos.prime(t), 151);
}
}
@@ -320,17 +335,20 @@ void test_quartic()
}
std::vector<Real> dydx(x.size());
for (size_t i = 0; i < y.size(); ++i) {
for (size_t i = 0; i < y.size(); ++i)
{
dydx[i] = 4*x[i]*x[i]*x[i];
}
std::vector<Real> d2ydx2(x.size());
for (size_t i = 0; i < y.size(); ++i) {
for (size_t i = 0; i < y.size(); ++i)
{
d2ydx2[i] = 12*x[i]*x[i];
}
std::vector<Real> d3ydx3(x.size());
for (size_t i = 0; i < y.size(); ++i) {
for (size_t i = 0; i < y.size(); ++i)
{
d3ydx3[i] = 24*x[i];
}
@@ -345,7 +363,8 @@ void test_quartic()
dydx.resize(10);
d2ydx2.resize(10);
d3ydx3.resize(10);
for (size_t i = 0; i < y.size(); ++i) {
for (size_t i = 0; i < y.size(); ++i)
{
y[i] = i*i*i*i;
dydx[i] = 4*i*i*i;
d2ydx2[i] = 12*i*i;
@@ -354,12 +373,15 @@ void test_quartic()
auto csh = cardinal_septic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3), Real(0), Real(1));
for (Real t = 1; t <= xmax; t += 0.0078125) {
for (Real t = 1; t <= xmax; t += 0.0078125)
{
CHECK_ULP_CLOSE(t*t*t*t, csh(t), 117);
CHECK_ULP_CLOSE(4*t*t*t, csh.prime(t), 117);
}
std::vector<std::array<Real, 4>> data(10);
for (size_t i = 0; i < data.size(); ++i) {
for (size_t i = 0; i < data.size(); ++i)
{
data[i][0] = i*i*i*i;
data[i][1] = 4*i*i*i;
data[i][2] = 12*i*i;
@@ -367,10 +389,11 @@ void test_quartic()
}
auto csh_aos = cardinal_septic_hermite_aos(std::move(data), Real(0), Real(1));
for (Real t = 1; t <= xmax; t += 0.0078125) {
for (Real t = 1; t <= xmax; t += 0.0078125)
{
CHECK_ULP_CLOSE(t*t*t*t, csh_aos(t), 117);
CHECK_ULP_CLOSE(4*t*t*t, csh_aos.prime(t), 117);
}
}