diff --git a/README.md b/README.md index dd874a1..65c4e8d 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ the following features. * available spline types: * **cubic C2 splines**: global, twice continuously differentiable * **cubic Hermite splines**: local, continuously differentiable (C1) -* boundary conditions: **first** and **second order** derivatives can be specified, periodic condition is not implemented +* boundary conditions: **first** and **second order** derivatives can be specified, **not-a-knot** condition, periodic condition is not implemented * extrapolation * linear: if first order derivatives are specified or 2nd order = 0 * quadratic: if 2nd order derivatives not equal to zero specified diff --git a/doc/spline.tex b/doc/spline.tex index 78dbe8b..ae3a729 100644 --- a/doc/spline.tex +++ b/doc/spline.tex @@ -209,6 +209,19 @@ \subsubsection*{Boundary conditions} $b_n$ by $\delta_n$ and replacing $b_{n-1}$ and $d_{n-1}$ by the relations in Equation~\eqref{eq:cubic_c2_spline_coeffs}. +The not-a-knot condition, i.e.\ the continuity of the third derivatives +at the two outer segments, gives $d_1=d_2$ and $d_{n-2}=d_{n-1}$: +\begin{equation*} +\begin{aligned} + f_1'''(x_2)&=f_2'''(x_2): & d_1&=d_2 \follows& + -h_2 c_1 + (h_1+h_2) c_2 - h_1 c_3 & = 0,\\ + f_{n-2}'''(x_{n-1})&=f_{n-1}'''(x_{n-1}): & d_{n-2}&=d_{n-1}\follows& + -h_{n-1} c_{n-2} + (h_{n-2}+h_{n-1}) c_{n-1} - h_{n-2} c_n & = 0. +\end{aligned} +\end{equation*} + + + In all cases, $d_n$ and $b_n$ are then determined by: \begin{equation*} \begin{aligned} @@ -267,9 +280,10 @@ \subsubsection*{Boundary conditions} \begin{equation*} \begin{aligned} f'(x_1)&=\delta_1: & b_1 & = \delta_1,\, c_0 = 0,\\ - f'(x_n)&=\delta_n: & b_n & = \delta_n,\, c_n = 0,\\ + f'(x_n)&=\delta_n: & b_n & = \delta_n,\, c_n = 0. \end{aligned} \end{equation*} + Second order conditions imply a value for $c$ but we can re-express this in terms of $b$: \begin{equation*} @@ -278,13 +292,31 @@ \subsubsection*{Boundary conditions} b_1 = \frac{1}{2}\left(-b_2 -\frac{\gamma}{2} h_1 + 3 \frac{y_2-y_1}{h_1}\right),\\ f_n''(x_n)&=\gamma_n: & c_n &= \frac{1}{2}\gamma_n,\\ - f_{n-1}''(x_n)&=\gamma_n: & 6 d_{n-1}h_{n-1} + 2 c_{n-1} &= \frac{1}{2}\gamma_n \equivalent + f_{n-1}''(x_n)&=\gamma_n: & 6 d_{n-1}h_{n-1} + 2 c_{n-1} &= \frac{1}{2}\gamma_n \equivalent b_n = \frac{1}{2}\left(-b_{n-1} +\frac{\gamma}{2} h_{n-1} + 3 \frac{y_n-y_{n-1}}{h_{n-1}}\right).\\ \end{aligned} \end{equation*} - +The not-a-knot contition gives $d_1=d_2$ and $d_{n-2}=d_{n-1}$ and +re-expressing in terms of $b$ yields: +\begin{equation*} +\begin{aligned} + f_1'''(x_2)&=f_2'''(x_2): & d_1&=d_2 \equivalent & + b_1 & = -b_2+2\frac{y_2-y_1}{h_1}+\frac{h_1^2}{h_2^2} + \left(b_2+b_3-2\frac{y_3-y_2}{h_2}\right),\\ + f_{n-2}'''(x_{n-1})&=f_{n-1}'''(x_{n-1}): & d_{n-2}&=d_{n-1}\equivalent& + b_n & = -b_{n-1}+2\frac{y_n-y_{n-1}}{h_{n-1}}\\ + & & & & & +\frac{h_{n-1}^2}{h_{n-2}^2} + \left(b_{n-2}+b_{n-1}-2\frac{y_{n-1}-y_{n-2}}{h_{n-2}}\right). +\end{aligned} +\end{equation*} +To determine a value for $c_n$ we additionally impose second order +continuity at the boundary: +\begin{equation*} + f_{n-1}''(x_n)=f_n''(x_n): \, c_n =3d_{n-1}h_{n-1}+c_{n-1} \equivalent\, + c_n=\frac{b_{n-1}+2 b_n}{h_{n-1}}-3\frac{y_n-y_{n-1}}{h_{n-1}^2}. +\end{equation*} \subsection{Monotonic splines} The spline interpolating function $f$ is monotonic increasing diff --git a/examples/tests/unit_tests.cpp b/examples/tests/unit_tests.cpp index bd8f0f8..7f241ba 100644 --- a/examples/tests/unit_tests.cpp +++ b/examples/tests/unit_tests.cpp @@ -14,6 +14,39 @@ #define BOOST_TEST_MODULE SplineUnitTests #include +// cubic function value and its first derivative +double f(double a, double b, double c, double d, double x) +{ + return (((d*x) + c)*x + b)*x + a; +} +double f1(double b, double c, double d, double x) +{ + return ((3.0*d*x) + 2.0*c)*x + b; +} + +// uniform [0,1] distribution +double rand_unif() +{ + return (double)std::rand()/RAND_MAX; +} +// exponential distribution +double rand_exp(double lambda=1.0) +{ + double u = rand_unif(); + return -log(u)/lambda; +} +// double exponential distribution +double rand_laplace(double lambda=1.0, double mu=0) +{ + double u = rand_unif(); + if(u<=0.5) { + return mu-log(2.0*u)/lambda; + } else { + return mu+log(2.0*(u-0.5))/lambda; + } +} + + // setup different types of splines for testing // which are at least "cont_deriv" times continuously differentiable @@ -35,6 +68,13 @@ std::vector setup_splines(const std::vector& X, splines.back().set_boundary(tk::spline::second_deriv, 0.2, tk::spline::second_deriv, -0.3); splines.back().set_points(X,Y); + // not-a-knot cubic spline + if(X.size()>=4) { + splines.push_back(tk::spline()); + splines.back().set_boundary(tk::spline::not_a_knot, 0.0, + tk::spline::not_a_knot, 0.0); + splines.back().set_points(X,Y); + } } if(cont_deriv<=1) { // C^1 splines @@ -72,6 +112,13 @@ std::vector setup_splines(const std::vector& X, splines.back().set_boundary(tk::spline::first_deriv, 0.0, tk::spline::first_deriv, -1.2); splines.back().set_points(X,Y,tk::spline::cspline_hermite); + if(X.size()>=4) { + // cubic hermite spline: not-a-knot + splines.push_back(tk::spline()); + splines.back().set_boundary(tk::spline::not_a_knot, 0.0, + tk::spline::not_a_knot, -1.2); + splines.back().set_points(X,Y,tk::spline::cspline_hermite); + } } if(cont_deriv<=0) { // C^0 splines @@ -119,8 +166,8 @@ std::vector evaluation_grid(const std::vector& X, // checks spline is C^0, i.e. continuous and goes through all grid points BOOST_AUTO_TEST_CASE( Continuity ) { - const double dx = 1e-15; // step size to check continuity - const double max_deriv = 10.0; // max 1st deriv of below spline + const double dx = 5e-15; // step size to check continuity + const double max_deriv = 150.0; // max 1st deriv of below spline // using well posed grid data for spline interpolation std::vector X = {-10.2, 0.1, 0.8, 2.4, 3.1, 3.2, 4.7, 19.1}; std::vector Y = { 2.7, -1.2, -0.5, 1.5, -1.0, -0.7, 1.2, -1.3}; @@ -151,7 +198,7 @@ BOOST_AUTO_TEST_CASE( Continuity ) BOOST_AUTO_TEST_CASE( Differentiability ) { const double dx = 1e-15; // step size to check continuity - const double h = 2e-8; // step size for finite differences + const double h = 8e-8; // step size for finite differences const double max_deriv = 100.0; // max 2nd deriv of below spline std::vector X = {-10.2, 0.1, 0.8, 2.4, 3.1, 3.2, 4.7, 19.1}; std::vector Y = { 2.7, -1.2, -0.5, 1.5, -1.0, -0.7, 1.2, -1.3}; @@ -187,7 +234,7 @@ BOOST_AUTO_TEST_CASE( Differentiability ) BOOST_AUTO_TEST_CASE( Smoothness ) { const double dx = 1e-15; // step size to check continuity - const double h = 3e-6; // step size for finite differences + const double h = 6e-6; // step size for finite differences const double max_deriv = 1000.0; // max 3nd deriv of below spline std::vector X = {-10.2, 0.1, 0.8, 2.4, 3.1, 3.2, 4.7, 19.1}; std::vector Y = { 2.7, -1.2, -0.5, 1.5, -1.0, -0.7, 1.2, -1.3}; @@ -224,7 +271,7 @@ BOOST_AUTO_TEST_CASE( ThirdDeriv ) // 2nd derivative is piecewise linear, so we can use a larger // step size for finite differences for 3rd derivative const double h = 1e-5; // step size for finite differences - const double tol = 1e-10; // error tolerance in 3rd derivative + const double tol = 2e-10; // error tolerance in 3rd derivative std::vector X = {-10.2, 0.1, 0.8, 2.4, 3.1, 3.2, 4.7, 19.1}; std::vector Y = { 2.7, -1.2, -0.5, 1.5, -1.0, -0.7, 1.2, -1.3}; @@ -298,45 +345,89 @@ BOOST_AUTO_TEST_CASE( Monotonicity ) BOOST_AUTO_TEST_CASE( BoundaryTypes ) { - const double dx = 1e-5; // step size to check continuity + const double dx = 1e-10; // step size to check continuity + const double tol = 1e-14; // error tolerance std::vector X = {-10.2, 0.1, 0.8, 2.4, 3.1, 3.2, 4.7, 19.1}; std::vector Y = { 2.7, -1.2, -0.5, 1.5, -1.0, -0.7, 1.2, -1.3}; - { + std::vector types = {tk::spline::cspline, tk::spline::cspline_hermite}; + + for(tk::spline::spline_type type : types) { // natural boundary conditions, i.e. 2nd derivative is zero tk::spline s; - s.set_points(X,Y); - BOOST_CHECK_CLOSE( s.deriv(2,X[0]-dx), 0.0, 0.0 ); // 2nd deriv = 0 - BOOST_CHECK_CLOSE( s.deriv(2,X.back()+dx), 0.0, 0.0 ); // 2nd deriv = 0 + s.set_points(X,Y,type); + BOOST_CHECK_SMALL( s.deriv(2,X[0]-dx), tol ); // 2nd deriv = 0 + BOOST_CHECK_SMALL( s.deriv(2,X.back()+dx), tol ); // 2nd deriv = 0 } - { + for(tk::spline::spline_type type : types) { // left and right with given 2nd and 1st derivative, respectively tk::spline s; double deriv1=1.3; double deriv2=-0.7; s.set_boundary(tk::spline::second_deriv, deriv2, tk::spline::first_deriv, deriv1); - s.set_points(X,Y); + s.set_points(X,Y,type); BOOST_CHECK_CLOSE( s.deriv(2,X[0]), deriv2, 0.0 ); BOOST_CHECK_CLOSE( s.deriv(2,X[0]-1.0), deriv2, 0.0 ); - BOOST_CHECK_CLOSE( s.deriv(1,X.back()), deriv1, 1e-12); - BOOST_CHECK_CLOSE( s.deriv(1,X.back()+1.0), deriv1, 1e-12); + BOOST_CHECK_CLOSE( s.deriv(1,X.back()), deriv1, tol*100.0); + BOOST_CHECK_CLOSE( s.deriv(1,X.back()+1.0), deriv1, tol*100.0); BOOST_CHECK_CLOSE( s.deriv(2,X.back()+1.0), 0.0, 0.0); } - { + for(tk::spline::spline_type type : types) { // left and right with given 1st and 2nd derivative, respectively tk::spline s; double deriv1=-4.7; double deriv2=-1.2; s.set_boundary(tk::spline::first_deriv, deriv1, tk::spline::second_deriv, deriv2); - s.set_points(X,Y); + s.set_points(X,Y,type); BOOST_CHECK_CLOSE( s.deriv(2,X.back()), deriv2, 0.0 ); BOOST_CHECK_CLOSE( s.deriv(2,X.back()+1.0), deriv2, 0.0 ); - BOOST_CHECK_CLOSE( s.deriv(1,X[0]), deriv1, 1e-12); - BOOST_CHECK_CLOSE( s.deriv(1,X[0]-1.0), deriv1, 1e-12); + BOOST_CHECK_CLOSE( s.deriv(1,X[0]), deriv1, tol*100.0); + BOOST_CHECK_CLOSE( s.deriv(1,X[0]-1.0), deriv1, tol*100.0); BOOST_CHECK_CLOSE( s.deriv(2,X[0]-1.0), 0.0, 0.0); } + + for(tk::spline::spline_type type : types) { + // left and right with not-a-knot condition + assert(X.size()>=4); + tk::spline s; + s.set_boundary(tk::spline::not_a_knot, 0.0, + tk::spline::not_a_knot, 0.0); + s.set_points(X,Y,type); + // f''' is continuous around X[1] and X[n-2] + int n=(int) X.size(); + BOOST_CHECK_SMALL( s.deriv(3,X[1]-dx)-s.deriv(3,X[1]+dx), tol ); + BOOST_CHECK_SMALL( s.deriv(3,X[n-2]-dx)-s.deriv(3,X[n-2]+dx), tol ); + // f'' continuous around X[0] and X[n-1] + BOOST_CHECK_SMALL( s.deriv(2,X[0]-dx)-s.deriv(2,X[0]+dx), 10.0*dx ); + BOOST_CHECK_SMALL( s.deriv(2,X[n-1]-dx)-s.deriv(2,X[n-1]+dx), 10.0*dx ); + } + { + // C^2 not-a-knot splines generated from cubic functions + // will be identical to the cubic function + double a=2.0, b=-2.7, c=-0.3, d=0.1; + for(size_t size : {4, 5, 6, 7, 20, 55, 57, 93, 121, 283}) { + std::vector XX(size), YY(size); + XX[0]=rand_laplace(); + YY[0]=f(a,b,c,d,XX[0]); + for(size_t i=1; i& x, spline_type type) { assert(x.size()==y.size()); - assert(x.size()>2); + assert(x.size()>=3); + // not-a-knot with 3 points has many solutions + if(m_left==not_a_knot || m_right==not_a_knot) + assert(x.size()>=4); m_type=type; m_made_monotonic=false; m_x=x; @@ -273,7 +277,9 @@ void spline::set_points(const std::vector& x, // setting up the matrix and right hand side of the equation system // for the parameters b[] - internal::band_matrix A(n,1,1); + int n_upper = (m_left == spline::not_a_knot) ? 2 : 1; + int n_lower = (m_right == spline::not_a_knot) ? 2 : 1; + internal::band_matrix A(n,n_upper,n_lower); std::vector rhs(n); for(int i=1; i& x, A(0,0)=2.0*(x[1]-x[0]); A(0,1)=1.0*(x[1]-x[0]); rhs[0]=3.0*((y[1]-y[0])/(x[1]-x[0])-m_left_value); + } else if(m_left == spline::not_a_knot) { + // f'''(x[1]) exists, i.e. d[0]=d[1], or re-expressed in c: + // -h1*c[0] + (h0+h1)*c[1] - h0*c[2] = 0 + A(0,0) = -(x[2]-x[1]); + A(0,1) = x[2]-x[0]; + A(0,2) = -(x[1]-x[0]); + rhs[0] = 0.0; } else { assert(false); } @@ -308,6 +321,13 @@ void spline::set_points(const std::vector& x, A(n-1,n-1)=2.0*(x[n-1]-x[n-2]); A(n-1,n-2)=1.0*(x[n-1]-x[n-2]); rhs[n-1]=3.0*(m_right_value-(y[n-1]-y[n-2])/(x[n-1]-x[n-2])); + } else if(m_right == spline::not_a_knot) { + // f'''(x[n-2]) exists, i.e. d[n-3]=d[n-2], or re-expressed in c: + // -h_{n-2}*c[n-3] + (h_{n-3}+h_{n-2})*c[n-2] - h_{n-3}*c[n-1] = 0 + A(n-1,n-3) = -(x[n-1]-x[n-2]); + A(n-1,n-2) = x[n-1]-x[n-3]; + A(n-1,n-1) = -(x[n-2]-x[n-3]); + rhs[0] = 0.0; } else { assert(false); } @@ -352,6 +372,12 @@ void spline::set_points(const std::vector& x, } else if(m_left==second_deriv) { const double h = m_x[1]-m_x[0]; m_b[0]=0.5*(-m_b[1]-0.5*m_left_value*h+3.0*(m_y[1]-m_y[0])/h); + } else if(m_left == not_a_knot) { + // f''' continuous at x[1] + const double h0 = m_x[1]-m_x[0]; + const double h1 = m_x[2]-m_x[1]; + m_b[0]= -m_b[1] + 2.0*(m_y[1]-m_y[0])/h0 + + h0*h0/(h1*h1)*(m_b[1]+m_b[2]-2.0*(m_y[2]-m_y[1])/h1); } else { assert(false); } @@ -362,6 +388,14 @@ void spline::set_points(const std::vector& x, const double h = m_x[n-1]-m_x[n-2]; m_b[n-1]=0.5*(-m_b[n-2]+0.5*m_right_value*h+3.0*(m_y[n-1]-m_y[n-2])/h); m_c[n-1]=0.5*m_right_value; + } else if(m_right == not_a_knot) { + // f''' continuous at x[n-2] + const double h0 = m_x[n-2]-m_x[n-3]; + const double h1 = m_x[n-1]-m_x[n-2]; + m_b[n-1]= -m_b[n-2] + 2.0*(m_y[n-1]-m_y[n-2])/h1 + h1*h1/(h0*h0) + *(m_b[n-3]+m_b[n-2]-2.0*(m_y[n-2]-m_y[n-3])/h0); + // f'' continuous at x[n-1]: c[n-1] = 3*d[n-2]*h[n-2] + c[n-1] + m_c[n-1]=(m_b[n-2]+2.0*m_b[n-1])/h1-3.0*(m_y[n-1]-m_y[n-2])/(h1*h1); } else { assert(false); }