171 lines
7.9 KiB
Plaintext
171 lines
7.9 KiB
Plaintext
[section:recurrence Tools For 3-Term Recurrence Relations]
|
|
|
|
[h4 Synopsis]
|
|
|
|
``
|
|
#include <boost/math/tools/recurrence.hpp>
|
|
``
|
|
|
|
namespace boost{ namespace math{ namespace tools{
|
|
|
|
template <class Recurrence, class T>
|
|
T function_ratio_from_backwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter);
|
|
|
|
template <class Recurrence, class T>
|
|
T function_ratio_from_forwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter);
|
|
|
|
template <class NextCoefs, class T>
|
|
T apply_recurrence_relation_forward(const NextCoefs& get_coefs, unsigned number_of_steps, T first, T second, int* log_scaling = 0, T* previous = 0);
|
|
|
|
template <class T, class NextCoefs>
|
|
T apply_recurrence_relation_backward(const NextCoefs& get_coefs, unsigned number_of_steps, T first, T second, int* log_scaling = 0, T* previous = 0);
|
|
|
|
template <class Recurrence>
|
|
struct forward_recurrence_iterator;
|
|
|
|
template <class Recurrence>
|
|
struct backward_recurrence_iterator;
|
|
|
|
}}} // namespaces
|
|
|
|
[h4 Description]
|
|
|
|
All of the tools in this header require a description of the recurrence relation: this takes the form of
|
|
a functor that returns a tuple containing the 3 coefficients, specifically, given a recurrence relation:
|
|
|
|
[/\Large $$ a_nF_{n-1} + b_nF_n + c_nF_{n+1} = 0 $$]
|
|
[equation three_term_recurrence]
|
|
|
|
And a functor `F` then the expression:
|
|
|
|
[expression F(n);]
|
|
|
|
Returns a tuple containing [role serif_italic { a[sub n], b[sub n], c[sub n] }].
|
|
|
|
For example, the recurrence relation for the Bessel J and Y functions when written in this form is:
|
|
|
|
[/\Large $$ J_{v-1}(x) - \frac{2v}{x}J_v(x) + J_{v+1}(x)= 0 $$]
|
|
[$../equations/three_term_recurrence_bessel_jy.svg]
|
|
|
|
Therefore, given local variables /x/ and /v/ of type `double` the recurrence relation for Bessel J and Y can be encoded
|
|
in a lambda expression like this:
|
|
|
|
auto recurrence_functor_jy = [&](int n) { return std::make_tuple(1.0, -2 * (v + n) / x, 1.0); };
|
|
|
|
Similarly, the Bessel I and K recurrence relation differs just by the sign of the final term:
|
|
|
|
[/\Large $$ I_{v-1}(x) - \frac{2v}{x}I_v(x) - I_{v+1}(x)= 0 $$]
|
|
[$../equations/three_term_recurrence_bessel_ik.svg]
|
|
|
|
And this could be encoded as:
|
|
|
|
auto recurrence_functor_ik = [&](int n) { return std::make_tuple(1.0, -2 * (v + n) / x, -1.0); };
|
|
|
|
The tools are then as follows:
|
|
|
|
template <class Recurrence, class T>
|
|
T function_ratio_from_backwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter);
|
|
|
|
Given a functor `r` which encodes the recurrence relation for function `F` at some location /n/, then returns the ratio:
|
|
|
|
[/\Large $$ F_n / F_{n-1} $$]
|
|
[$../equations/three_term_recurrence_backwards_ratio.svg]
|
|
|
|
This calculation is stable only if recurrence is stable in the backwards direction. Further the ratio calculated
|
|
is for the dominant solution (in the backwards direction) of the recurrence relation, if there are multiple solutions,
|
|
then there is no guarantee that this will find the one you want or expect.
|
|
|
|
Argument /factor/ is the tolerance required for convergence of the continued fraction associated with
|
|
the recurrence relation, and should be no smaller than machine epsilon. Argument /max_iter/ sets
|
|
the maximum number of permitted iterations in the associated continued fraction.
|
|
|
|
template <class Recurrence, class T>
|
|
T function_ratio_from_forwards_recurrence(const Recurrence& r, const T& factor, boost::uintmax_t& max_iter);
|
|
|
|
Given a functor `r` which encodes the recurrence relation for function F at some location /n/, then returns the ratio:
|
|
|
|
[/\Large $$ F_n / F_{n+1} $$]
|
|
[$../equations/three_term_recurrence_forwards_ratio.svg]
|
|
|
|
This calculation is stable only if recurrence is stable in the forwards direction. Further the ratio calculated
|
|
is for the dominant solution (in the forwards direction) of the recurrence relation, if there are multiple solutions,
|
|
then there is no guarantee that this will find the one you want or expect.
|
|
|
|
Argument /factor/ is the tolerance required for convergence of the continued fraction associated with
|
|
the recurrence relation, and should be no smaller than machine epsilon. Argument /max_iter/ sets
|
|
the maximum number of permitted iterations in the associated continued fraction.
|
|
|
|
template <class NextCoefs, class T>
|
|
T apply_recurrence_relation_forward(const NextCoefs& get_coefs, unsigned number_of_steps, T first, T second, int* log_scaling = 0, T* previous = 0);
|
|
|
|
Applies a recurrence relation in a stable forward direction, starting with the values F[sub n-1] and F[sub n].
|
|
|
|
[variablelist
|
|
[[get_coefs] [Functor that returns the corefficients of the recurrence relation. The coefficients should be centered on position /second/.]]
|
|
[[number_of_steps][The number of steps to apply the recurrence relation onwards from /second/.]]
|
|
[[first] [The value of F[sub n-1]]]
|
|
[[second] [The value of F[sub n]]]
|
|
[[log_scaling][When provided, the recurrence relations may be rescaled internally to avoid over/underflow issues. The result should be multiplied by `exp(*log_scaling)` to get the true value of the result.]]
|
|
[[previous][When provided, is set to the value of F[sub n + number_of_steps - 1]]
|
|
]
|
|
]
|
|
|
|
Returns F[sub n + number_of_steps].
|
|
|
|
template <class NextCoefs, class T>
|
|
T apply_recurrence_relation_backward(const NextCoefs& get_coefs, unsigned number_of_steps, T first, T second, int* log_scaling = 0, T* previous = 0);
|
|
|
|
Applies a recurrence relation in a stable backward direction, starting with the values F[sub n+1] and F[sub n].
|
|
|
|
[variablelist
|
|
[[get_coefs] [Functor that returns the corefficients of the recurrence relation. The coefficients should be centered on position /second/.]]
|
|
[[number_of_steps][The number of steps to apply the recurrence relation backwards from /second/.]]
|
|
[[first] [The value of F[sub n+1]]]
|
|
[[second] [The value of F[sub n]]]
|
|
[[log_scaling][When provided, the recurrence relations may be rescaled internally to avoid over/underflow issues. The result should be multiplied by `exp(*log_scaling)` to get the true value of the result.]]
|
|
[[previous][When provided, is set to the value of F[sub n - number_of_steps + 1]]
|
|
]
|
|
]
|
|
|
|
Returns F[sub n - number_of_steps].
|
|
|
|
template <class Recurrence>
|
|
struct forward_recurrence_iterator
|
|
{
|
|
typedef typename boost::remove_reference<decltype(std::get<0>(std::declval<Recurrence&>()(0)))>::type value_type;
|
|
|
|
forward_recurrence_iterator(const Recurrence& r, value_type f_n_minus_1, value_type f_n);
|
|
forward_recurrence_iterator(const Recurrence& r, value_type f_n);
|
|
/* Operators omitted for clarity */
|
|
};
|
|
|
|
Type `forward_recurrence_iterator` defines a forward-iterator for a recurrence relation stable in the
|
|
forward direction. The constructors take the recurrence relation, plus either one or two values: if
|
|
only one value is provided, then the second is computed by using the recurrence relation to calculate the function ratio.
|
|
|
|
template <class Recurrence>
|
|
struct backward_recurrence_iterator
|
|
{
|
|
typedef typename boost::remove_reference<decltype(std::get<0>(std::declval<Recurrence&>()(0)))>::type value_type;
|
|
|
|
backward_recurrence_iterator(const Recurrence& r, value_type f_n_plus_1, value_type f_n);
|
|
backward_recurrence_iterator(const Recurrence& r, value_type f_n);
|
|
/* Operators omitted for clarity */
|
|
};
|
|
|
|
Type `backward_recurrence_iterator` defines a forward-iterator for a recurrence relation stable in the
|
|
backward direction. The constructors take the recurrence relation, plus either one or two values: if
|
|
only one value is provided, then the second is computed by using the recurrence relation to calculate the function ratio.
|
|
|
|
Note that /incrementing/ this iterator moves the value returned successively to F[sub n-1], F[sub n-2] etc.
|
|
|
|
[endsect] [/section:recurrence Tools For 3-Term Recurrence Relations]
|
|
|
|
[/
|
|
Copyright 2019 John Maddock.
|
|
Distributed under the Boost Software License, Version 1.0.
|
|
(See accompanying file LICENSE_1_0.txt or copy at
|
|
http://www.boost.org/LICENSE_1_0.txt).
|
|
]
|
|
|