From edfb80e76dece4059b22d8ee647a6ae0eb2b5403 Mon Sep 17 00:00:00 2001 From: NAThompson Date: Wed, 7 Aug 2019 14:10:08 -0400 Subject: [PATCH 1/3] Cubic B-spline second derivatives --- .../math/interpolators/cubic_b_spline.hpp | 9 ++++ .../detail/cubic_b_spline_detail.hpp | 43 +++++++++++++++++++ test/test_cubic_b_spline.cpp | 7 +++ 3 files changed, 59 insertions(+) diff --git a/include/boost/math/interpolators/cubic_b_spline.hpp b/include/boost/math/interpolators/cubic_b_spline.hpp index 73ac1d013..91a435aee 100644 --- a/include/boost/math/interpolators/cubic_b_spline.hpp +++ b/include/boost/math/interpolators/cubic_b_spline.hpp @@ -45,6 +45,8 @@ public: Real prime(Real x) const; + Real double_prime(Real x) const; + private: std::shared_ptr> m_imp; }; @@ -74,5 +76,12 @@ Real cubic_b_spline::prime(Real x) const return m_imp->prime(x); } +template +Real cubic_b_spline::double_prime(Real x) const +{ + return m_imp->double_prime(x); +} + + }} #endif diff --git a/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp b/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp index de013fabb..1347cbc30 100644 --- a/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp +++ b/include/boost/math/interpolators/detail/cubic_b_spline_detail.hpp @@ -32,6 +32,8 @@ public: Real prime(Real x) const; + Real double_prime(Real x) const; + private: std::vector m_beta; Real m_h_inv; @@ -79,6 +81,25 @@ Real b3_spline_prime(Real x) return (Real) 0; } +template +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()*x - 2) + x*(3*boost::math::constants::half()); + } + if (x < 2) + { + return (2 - x); + } + return (Real) 0; +} + template template @@ -283,5 +304,27 @@ Real cubic_b_spline_imp::prime(Real x) const return z*m_h_inv; } +template +Real cubic_b_spline_imp::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(0), boost::math::ltrunc(ceil(t - 2))); + size_t k_max = (size_t) (min)(static_cast(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 diff --git a/test/test_cubic_b_spline.cpp b/test/test_cubic_b_spline.cpp index 8655cabf2..8f4e66c65 100644 --- a/test/test_cubic_b_spline.cpp +++ b/test/test_cubic_b_spline.cpp @@ -30,6 +30,9 @@ void test_b3_spline() BOOST_CHECK_SMALL(boost::math::detail::b3_spline(-2.5), (Real) 0); BOOST_CHECK_SMALL(boost::math::detail::b3_spline_prime(2.5), (Real) 0); BOOST_CHECK_SMALL(boost::math::detail::b3_spline_prime(-2.5), (Real) 0); + BOOST_CHECK_SMALL(boost::math::detail::b3_spline_double_prime(2.5), (Real) 0); + BOOST_CHECK_SMALL(boost::math::detail::b3_spline_double_prime(-2.5), (Real) 0); + // On the boundary of support: BOOST_CHECK_SMALL(boost::math::detail::b3_spline(2), (Real) 0); @@ -52,6 +55,7 @@ void test_b3_spline() Real arg = i*0.01; BOOST_CHECK_CLOSE(boost::math::detail::b3_spline(arg), boost::math::detail::b3_spline(arg), eps); BOOST_CHECK_CLOSE(boost::math::detail::b3_spline_prime(-arg), -boost::math::detail::b3_spline_prime(arg), eps); + BOOST_CHECK_CLOSE(boost::math::detail::b3_spline_double_prime(-arg), boost::math::detail::b3_spline_double_prime(arg), eps); } } @@ -109,6 +113,9 @@ void test_constant_function() BOOST_CHECK_CLOSE(y, constant, 10*std::numeric_limits::epsilon()); Real y_prime = spline.prime(i*step + a + 0.002); BOOST_CHECK_SMALL(y_prime, 5000*std::numeric_limits::epsilon()); + Real y_double_prime = spline.double_prime(i*step + a + 0.002); + BOOST_CHECK_SMALL(y_double_prime, 5000*std::numeric_limits::epsilon()); + } // Test that correctly specified left and right-derivatives work properly: From b28d3f34ea8c2debba4f5cf9aa5de18a60aa1373 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Thu, 8 Aug 2019 11:31:26 +0100 Subject: [PATCH 2/3] hypergeometrics: indent equations. [CI SKIP] --- doc/sf/hypergeometric.qbk | 60 +++++++++++++++++++-------------------- 1 file changed, 30 insertions(+), 30 deletions(-) diff --git a/doc/sf/hypergeometric.qbk b/doc/sf/hypergeometric.qbk index 858584e9e..bc7462326 100644 --- a/doc/sf/hypergeometric.qbk +++ b/doc/sf/hypergeometric.qbk @@ -169,25 +169,25 @@ Otherwise we use the defining series. The function `hypergeometric_1F1(a, b, z)` returns the non-singular solution to [@https://en.wikipedia.org/wiki/Confluent_hypergeometric_function Kummer's equation] -[/\large $$z \frac{d^2 w}{d z^2} + (b-z) \frac{dw}{dz} - aw = 0 $$] -[$../equations/hypergeometric_1f1_2.svg] +[:[/\large $$z \frac{d^2 w}{d z^2} + (b-z) \frac{dw}{dz} - aw = 0 $$] +[$../equations/hypergeometric_1f1_2.svg]] which for |/z/| < 1 has the hypergeometric series expansion -[$../equations/hypergeometric_1f1_1.svg] +[:[$../equations/hypergeometric_1f1_1.svg]] where (/a/)[sub /n/] denotes rising factorial. This function has the same definition as Mathematica's `Hypergeometric1F1[a, b, z]` and Maple's `KummerM(a, b, z)`. The "regularized" versions of the function return: -[/ \Large $$ \textbf{M}(a, b; z) = \frac{{_1F_1}(a, b; z)}{\Gamma(b)} = \sum_{n=0}^{\infty} \frac{(a)_n z^n}{\Gamma(b+n) n!} $$] -[$../equations/hypergeometric_1f1_17.svg] +[:[/ \Large $$ \textbf{M}(a, b; z) = \frac{{_1F_1}(a, b; z)}{\Gamma(b)} = \sum_{n=0}^{\infty} \frac{(a)_n z^n}{\Gamma(b+n) n!} $$] +[$../equations/hypergeometric_1f1_17.svg]] The "log" versions of the function return: -[/ \Large $$ \ln(|_1F_1(a, b, z)|) $$ ] -[$../equations/hypergeometric_1f1_18.svg] +[:[/ \Large $$ \ln(|_1F_1(a, b, z)|) $$ ] +[$../equations/hypergeometric_1f1_18.svg]] When the optional `sign` argument is supplied it is set on exit to either +1 or -1 depending on the sign of [sub 1]/F/[sub 1](/a/, /b/, /z/). @@ -276,12 +276,12 @@ readers will have to refer to the code for the transitions between methods, as w In almost all cases where /z < 0/ we use [@https://en.wikipedia.org/wiki/Confluent_hypergeometric_function Kummer's relation] to make /z/ positive: -[/\large $$_1F_1(a, b; -z) = e^{-z} {_1F_1}(a, b; z)$$] -[$../equations/hypergeometric_1f1_12.svg] +[:[/\large $$_1F_1(a, b; -z) = e^{-z} {_1F_1}(a, b; z)$$] +[$../equations/hypergeometric_1f1_12.svg]] The main series representation -[$../equations/hypergeometric_1f1_1.svg] +[:[$../equations/hypergeometric_1f1_1.svg]] is used only when @@ -293,29 +293,29 @@ to life again as /b/ crosses the origin. A&S 13.3.6 gives -[/\large $$_1F_1(a, b; z) = e^{ \frac{z}{2} } \Gamma(b-a- \frac{1}{2} ) ( \frac{z}{4})^{a-b+ \frac{1}{2}} \sum_{n=0}^{\infty} \frac{(2b-2a-1)_n(b-2a)_n(b-a-\frac{1}{2}+n)}{n!(b)_n} (-1)^n I_{b-a-\frac{1}{2}+n}(\frac{z}{2})$$] -[$../equations/hypergeometric_1f1_3.svg] +[:[/\large $$_1F_1(a, b; z) = e^{ \frac{z}{2} } \Gamma(b-a- \frac{1}{2} ) ( \frac{z}{4})^{a-b+ \frac{1}{2}} \sum_{n=0}^{\infty} \frac{(2b-2a-1)_n(b-2a)_n(b-a-\frac{1}{2}+n)}{n!(b)_n} (-1)^n I_{b-a-\frac{1}{2}+n}(\frac{z}{2})$$] +[$../equations/hypergeometric_1f1_3.svg]] which is particularly useful for ['a [cong] b] and ['z > 0], or /a/ \u226a 1, and ['z < 0]. Then we have Tricomi's approximation (given in simplified form in A&S 13.3.7) [link math_toolkit.hypergeometric.hypergeometric_refs (7)] -[/\large $$_1F_1(a, b; z) = \Gamma(b) e^{ \frac{1}{2} z} \sum_{n=0}^{\infty} 2^{-n}z^n A_n(a,b) E_{b-1+n}(z(\frac{b}{2}-a)) $$] -[$../equations/hypergeometric_1f1_4.svg] +[:[/\large $$_1F_1(a, b; z) = \Gamma(b) e^{ \frac{1}{2} z} \sum_{n=0}^{\infty} 2^{-n}z^n A_n(a,b) E_{b-1+n}(z(\frac{b}{2}-a)) $$] +[$../equations/hypergeometric_1f1_4.svg]] with -[/\large $$A_0(a,b) = 1, A_1(a,b) = 0, A_2(a,b) = \frac{b}{2}, A_3(a,b)= - \frac{1}{3}(b-2a) $$] -[$../equations/hypergeometric_1f1_5.svg] +[:[/\large $$A_0(a,b) = 1, A_1(a,b) = 0, A_2(a,b) = \frac{b}{2}, A_3(a,b)= - \frac{1}{3}(b-2a) $$] +[$../equations/hypergeometric_1f1_5.svg]] and -[/\large $$(n+1)A_{n+1} = (n+b-1)A_{n-1} - (b-2a)A_{n-2} \quad;\quad n \geq 2 $$] -[$../equations/hypergeometric_1f1_6.svg] +[:[/\large $$(n+1)A_{n+1} = (n+b-1)A_{n-1} - (b-2a)A_{n-2} \quad;\quad n \geq 2 $$] +[$../equations/hypergeometric_1f1_6.svg]] With ['E[sub v]] defined as: -[/ +[:[/ \begin{equation*} \begin{split} E_v(z) & = z^{-\frac{1}{2}v} J_v(2 \sqrt{z})\\ @@ -323,19 +323,19 @@ E_v(-z) & = z^{-\frac{1}{2}v} I_v(2 \sqrt{z})\\ E_v(0) & = \frac{1}{\Gamma(v + 1)} \end{split} \end{equation*}] -[$../equations/hypergeometric_1f1_7.svg] +[$../equations/hypergeometric_1f1_7.svg]] This approximation is especially effective when ['a < 0]. For /a/ and /b/ both large and positive and somewhat similar in size then an expansion in terms of gamma functions can be used [link math_toolkit.hypergeometric.hypergeometric_refs (6)]: -[/\large $$_1F_1(a, b; x) = \frac{1}{B(a, b-a)} e^x \sum_{k=0}^{\infty} \frac{(1-a)_k}{k!} \frac{\gamma(b-a+k, x)}{x^{b-a+k}} $$] -[$../equations/hypergeometric_1f1_8.svg] +[:[/\large $$_1F_1(a, b; x) = \frac{1}{B(a, b-a)} e^x \sum_{k=0}^{\infty} \frac{(1-a)_k}{k!} \frac{\gamma(b-a+k, x)}{x^{b-a+k}} $$] +[$../equations/hypergeometric_1f1_8.svg]] For /z/ large we have the asymptotic expansion: -[/\large $$_1F_1(a, b; x) \approx \frac{e^x x^{a-b}}{\Gamma(a)} \sum_{k=0}^{\infty} \frac{(1-a)_k(b-a)_k}{k! x^k} $$] -[$../equations/hypergeometric_1f1_9.svg] +[:[/\large $$_1F_1(a, b; x) \approx \frac{e^x x^{a-b}}{\Gamma(a)} \sum_{k=0}^{\infty} \frac{(1-a)_k(b-a)_k}{k! x^k} $$] +[$../equations/hypergeometric_1f1_9.svg]] which is a special case of the gamma function expansion above. @@ -348,8 +348,8 @@ That effectively completes the "direct" methods used, the remaining areas are tr These require extreme care in their use, as they often change the direction of stability, or else are not stable at all. Sometimes this is a well defined and characterized change in direction: for example with /a,b/ and /z/ all positive recurrence on /a/ via -[/\large $$(b-a) _1F_1(a-1, b; z) + (2a-b+z) _1F_1(a, b; z) -a _1F_1(a+1, b; z) = 0 $$] -[$../equations/hypergeometric_1f1_10.svg] +[:[/\large $$(b-a) _1F_1(a-1, b; z) + (2a-b+z) _1F_1(a, b; z) -a _1F_1(a+1, b; z) = 0 $$] +[$../equations/hypergeometric_1f1_10.svg]] abruptly changes from stable forwards recursion to stable backwards if /2a-b+z/ changes sign. Thus recurrence on /a/, even when [sub 1]/F/[sub 1](/a/+/N/, /b/, /z/) is strictly increasing, needs careful handling as /a/ \u2192 0. @@ -367,8 +367,8 @@ the continued fractions associated with the recurrence relations to calculate [s and then normalize either by iterating forwards on /b/ until /b > 0/ and comparison with a reference value, or by using the Wronskian -[/\large $$_1F_1(a, b; z) \frac{d}{dz}z^{1-b}_1F_1(1+a-b, 2-b; z) - z^{1-b}_1F_1(1+a-b, 2-b; z)\frac{d}{dz}{_1F_1}(a, b; z) = (1-b)z^{-b}e^z,$$] -[$../equations/hypergeometric_1f1_11.svg] +[:[/\large $$_1F_1(a, b; z) \frac{d}{dz}z^{1-b}_1F_1(1+a-b, 2-b; z) - z^{1-b}_1F_1(1+a-b, 2-b; z)\frac{d}{dz}{_1F_1}(a, b; z) = (1-b)z^{-b}e^z,$$] +[$../equations/hypergeometric_1f1_11.svg]] which is of course the well known Miller's method [link math_toolkit.hypergeometric.hypergeometric_refs (12)]. @@ -455,8 +455,8 @@ are destroyed by cancellation. The function `hypergeometric_pFq` returns the result of: -[/\Large $$ _pF_q(\{a_1...a_n\}, \{b_1...b_n\}; z) = \sum_{k=0}^{\infty} \frac{\Pi_{j=1}^n(a_j)_n z^k}{\Pi_{j=1}^n(b_j)_n k!} $$] -[$../equations/hypergeometric_pfq_1.svg] +[:[/\Large $$ _pF_q(\{a_1...a_n\}, \{b_1...b_n\}; z) = \sum_{k=0}^{\infty} \frac{\Pi_{j=1}^n(a_j)_n z^k}{\Pi_{j=1}^n(b_j)_n k!} $$] +[$../equations/hypergeometric_pfq_1.svg]] It is most naturally used via initializer lists as in: From d7772d1560badc4b253c3e6c757fd396f273a18d Mon Sep 17 00:00:00 2001 From: NAThompson Date: Thu, 8 Aug 2019 11:26:57 -0400 Subject: [PATCH 3/3] Document second derivative of cubic B-spline. [CI SKIP] --- doc/interpolators/cubic_b_spline.qbk | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/doc/interpolators/cubic_b_spline.qbk b/doc/interpolators/cubic_b_spline.qbk index de1475ff8..77dbd8bc7 100644 --- a/doc/interpolators/cubic_b_spline.qbk +++ b/doc/interpolators/cubic_b_spline.qbk @@ -30,6 +30,8 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) Real operator()(Real x) const; Real prime(Real x) const; + + Real double_prime(Real x) const; }; }} // namespaces @@ -37,7 +39,7 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) [heading Cubic B-Spline Interpolation] -The cubic B-spline class provided by boost allows fast and accurate interpolation of a function which is known at equally spaced points. +The cubic B-spline class provided by Boost allows fast and accurate interpolation of a function which is known at equally spaced points. The cubic B-spline interpolation is numerically stable as it uses compactly supported basis functions constructed via iterative convolution. This is to be contrasted to traditional cubic spline interpolation which is ill-conditioned as the global support of cubic polynomials causes small changes far from the evaluation point to exert a large influence on the calculated value. @@ -79,6 +81,15 @@ and to evaluate the derivative of the interpolant we use Be aware that the accuracy guarantees on the derivative of the spline are an order lower than the guarantees on the original function, see [@http://www.springer.com/us/book/9780387984087 Numerical Analysis, Graduate Texts in Mathematics, 181, Rainer Kress] for details. +The last interesting member is the second derivative, evaluated via + + double ypp = spline.double_prime(x); + +The basis functions of the spline are cubic polynomials, so the second derivative is simply linear interpolation. +But the interpolation is not constrained by data (you don't pass in the second derivatives), and hence is less accurate than would be expected from assuming it is a linear interpolation. +The problem is especially pronounced at the boundaries, where the second derivative is totally unconstrained. +Use the second derivative of the cubic B-spline interpolator only in desperation. + Finally, note that this is an interpolator, not an extrapolator. Therefore, you should strenuously avoid evaluating the spline outside the endpoints. However, it is not an error if you do, as often you cannot control where (say) an ODE stepper will evaluate your function.