266 lines
8.8 KiB
Plaintext
266 lines
8.8 KiB
Plaintext
[section:legendre Legendre (and Associated) Polynomials]
|
|
|
|
[h4 Synopsis]
|
|
|
|
``
|
|
#include <boost/math/special_functions/legendre.hpp>
|
|
``
|
|
|
|
namespace boost{ namespace math{
|
|
|
|
template <class T>
|
|
``__sf_result`` legendre_p(int n, T x);
|
|
|
|
template <class T, class ``__Policy``>
|
|
``__sf_result`` legendre_p(int n, T x, const ``__Policy``&);
|
|
|
|
template <class T>
|
|
``__sf_result`` legendre_p_prime(int n, T x);
|
|
|
|
template <class T, class ``__Policy``>
|
|
``__sf_result`` legendre_p_prime(int n, T x, const ``__Policy``&);
|
|
|
|
template <class T, class ``__Policy``>
|
|
std::vector<T> legendre_p_zeros(int l, const ``__Policy``&);
|
|
|
|
template <class T>
|
|
std::vector<T> legendre_p_zeros(int l);
|
|
|
|
template <class T>
|
|
``__sf_result`` legendre_p(int n, int m, T x);
|
|
|
|
template <class T, class ``__Policy``>
|
|
``__sf_result`` legendre_p(int n, int m, T x, const ``__Policy``&);
|
|
|
|
template <class T>
|
|
``__sf_result`` legendre_q(unsigned n, T x);
|
|
|
|
template <class T, class ``__Policy``>
|
|
``__sf_result`` legendre_q(unsigned n, T x, const ``__Policy``&);
|
|
|
|
template <class T1, class T2, class T3>
|
|
``__sf_result`` legendre_next(unsigned l, T1 x, T2 Pl, T3 Plm1);
|
|
|
|
template <class T1, class T2, class T3>
|
|
``__sf_result`` legendre_next(unsigned l, unsigned m, T1 x, T2 Pl, T3 Plm1);
|
|
|
|
|
|
}} // namespaces
|
|
|
|
The return type of these functions is computed using the __arg_promotion_rules:
|
|
note than when there is a single template argument the result is the same type
|
|
as that argument or `double` if the template argument is an integer type.
|
|
|
|
[optional_policy]
|
|
|
|
[h4 Description]
|
|
|
|
template <class T>
|
|
``__sf_result`` legendre_p(int l, T x);
|
|
|
|
template <class T, class ``__Policy``>
|
|
``__sf_result`` legendre_p(int l, T x, const ``__Policy``&);
|
|
|
|
Returns the Legendre Polynomial of the first kind:
|
|
|
|
[equation legendre_0]
|
|
|
|
Requires -1 <= x <= 1, otherwise returns the result of __domain_error.
|
|
|
|
Negative orders are handled via the reflection formula:
|
|
|
|
[:P[sub -l-1](x) = P[sub l](x)]
|
|
|
|
The following graph illustrates the behaviour of the first few
|
|
Legendre Polynomials:
|
|
|
|
[graph legendre_p]
|
|
|
|
template <class T>
|
|
``__sf_result`` legendre_p_prime(int n, T x);
|
|
|
|
template <class T, class ``__Policy``>
|
|
``__sf_result`` legendre_p_prime(int n, T x, const ``__Policy``&);
|
|
|
|
Returns the derivatives of the Legendre polynomials.
|
|
|
|
template <class T, class ``__Policy``>
|
|
std::vector<T> legendre_p_zeros(int l, const ``__Policy``&);
|
|
|
|
template <class T>
|
|
std::vector<T> legendre_p_zeros(int l);
|
|
|
|
The zeros of the Legendre polynomials are calculated by Newton's method using an initial guess given by Tricomi with root bracketing provided by Szego.
|
|
|
|
Since the Legendre polynomials are alternatively even and odd, only the non-negative zeros are returned.
|
|
For the odd Legendre polynomials, the first zero is always zero.
|
|
The rest of the zeros are returned in increasing order.
|
|
|
|
Note that the argument to the routine is an integer, and the output is a floating-point type.
|
|
Hence the template argument is mandatory.
|
|
The time to extract a single root is linear in `l` (this is scaling to evaluate the Legendre polynomials), so recovering all roots is [bigo](`l`[super 2]).
|
|
Algorithms with linear scaling [@ https://doi.org/10.1137/06067016X exist] for recovering all roots, but requires tooling not currently built into boost.math.
|
|
This implementation proceeds under the assumption that calculating zeros of these functions will not be a bottleneck for any workflow.
|
|
|
|
template <class T>
|
|
``__sf_result`` legendre_p(int l, int m, T x);
|
|
|
|
template <class T, class ``__Policy``>
|
|
``__sf_result`` legendre_p(int l, int m, T x, const ``__Policy``&);
|
|
|
|
Returns the associated Legendre polynomial of the first kind:
|
|
|
|
[equation legendre_1]
|
|
|
|
Requires -1 <= x <= 1, otherwise returns the result of __domain_error.
|
|
|
|
Negative values of /l/ and /m/ are handled via the identity relations:
|
|
|
|
[equation legendre_3]
|
|
|
|
[caution The definition of the associated Legendre polynomial used here
|
|
includes a leading Condon-Shortley phase term of (-1)[super m]. This
|
|
matches the definition given by Abramowitz and Stegun (8.6.6) and that
|
|
used by [@http://mathworld.wolfram.com/LegendrePolynomial.html Mathworld]
|
|
and [@http://documents.wolfram.com/mathematica/functions/LegendreP
|
|
Mathematica's LegendreP function]. However, uses in the literature
|
|
do not always include this phase term, and strangely the specification
|
|
for the associated Legendre function in the C++ TR1 (assoc_legendre)
|
|
also omits it, in spite of stating that it uses Abramowitz and Stegun
|
|
as the final arbiter on these matters.
|
|
|
|
See:
|
|
|
|
[@http://mathworld.wolfram.com/LegendrePolynomial.html
|
|
Weisstein, Eric W. "Legendre Polynomial."
|
|
From MathWorld--A Wolfram Web Resource].
|
|
|
|
Abramowitz, M. and Stegun, I. A. (Eds.). "Legendre Functions" and
|
|
"Orthogonal Polynomials." Ch. 22 in Chs. 8 and 22 in Handbook of
|
|
Mathematical Functions with Formulas, Graphs, and Mathematical Tables,
|
|
9th printing. New York: Dover, pp. 331-339 and 771-802, 1972.
|
|
]
|
|
|
|
template <class T>
|
|
``__sf_result`` legendre_q(unsigned n, T x);
|
|
|
|
template <class T, class ``__Policy``>
|
|
``__sf_result`` legendre_q(unsigned n, T x, const ``__Policy``&);
|
|
|
|
Returns the value of the Legendre polynomial that is the second solution
|
|
to the Legendre differential equation, for example:
|
|
|
|
[equation legendre_2]
|
|
|
|
Requires -1 <= x <= 1, otherwise __domain_error is called.
|
|
|
|
The following graph illustrates the first few Legendre functions of the
|
|
second kind:
|
|
|
|
[graph legendre_q]
|
|
|
|
template <class T1, class T2, class T3>
|
|
``__sf_result`` legendre_next(unsigned l, T1 x, T2 Pl, T3 Plm1);
|
|
|
|
Implements the three term recurrence relation for the Legendre
|
|
polynomials, this function can be used to create a sequence of
|
|
values evaluated at the same /x/, and for rising /l/. This recurrence
|
|
relation holds for Legendre Polynomials of both the first and second kinds.
|
|
|
|
[equation legendre_4]
|
|
|
|
For example we could produce a vector of the first 10 polynomial
|
|
values using:
|
|
|
|
double x = 0.5; // Abscissa value
|
|
vector<double> v;
|
|
v.push_back(legendre_p(0, x));
|
|
v.push_back(legendre_p(1, x));
|
|
for(unsigned l = 1; l < 10; ++l)
|
|
v.push_back(legendre_next(l, x, v[l], v[l-1]));
|
|
// Double check values:
|
|
for(unsigned l = 1; l < 10; ++l)
|
|
assert(v[l] == legendre_p(l, x));
|
|
|
|
Formally the arguments are:
|
|
|
|
[variablelist
|
|
[[l][The degree of the last polynomial calculated.]]
|
|
[[x][The abscissa value]]
|
|
[[Pl][The value of the polynomial evaluated at degree /l/.]]
|
|
[[Plm1][The value of the polynomial evaluated at degree /l-1/.]]
|
|
]
|
|
|
|
template <class T1, class T2, class T3>
|
|
``__sf_result`` legendre_next(unsigned l, unsigned m, T1 x, T2 Pl, T3 Plm1);
|
|
|
|
Implements the three term recurrence relation for the Associated Legendre
|
|
polynomials, this function can be used to create a sequence of
|
|
values evaluated at the same /x/, and for rising /l/.
|
|
|
|
[equation legendre_5]
|
|
|
|
For example we could produce a vector of the first m+10 polynomial
|
|
values using:
|
|
|
|
double x = 0.5; // Abscissa value
|
|
int m = 10; // order
|
|
vector<double> v;
|
|
v.push_back(legendre_p(m, m, x));
|
|
v.push_back(legendre_p(1 + m, m, x));
|
|
for(unsigned l = 1; l < 10; ++l)
|
|
v.push_back(legendre_next(l + 10, m, x, v[l], v[l-1]));
|
|
// Double check values:
|
|
for(unsigned l = 1; l < 10; ++l)
|
|
assert(v[l] == legendre_p(10 + l, m, x));
|
|
|
|
Formally the arguments are:
|
|
|
|
[variablelist
|
|
[[l][The degree of the last polynomial calculated.]]
|
|
[[m][The order of the Associated Polynomial.]]
|
|
[[x][The abscissa value]]
|
|
[[Pl][The value of the polynomial evaluated at degree /l/.]]
|
|
[[Plm1][The value of the polynomial evaluated at degree /l-1/.]]
|
|
]
|
|
|
|
[h4 Accuracy]
|
|
|
|
The following table shows peak errors (in units of epsilon)
|
|
for various domains of input arguments.
|
|
Note that only results for the widest floating point type on the system are
|
|
given as narrower types have __zero_error.
|
|
|
|
[table_legendre_p]
|
|
|
|
[table_legendre_q]
|
|
|
|
[table_legendre_p_associated_]
|
|
|
|
Note that the worst errors occur when the order increases, values greater than
|
|
~120 are very unlikely to produce sensible results, especially in the associated
|
|
polynomial case when the degree is also large. Further the relative errors
|
|
are likely to grow arbitrarily large when the function is very close to a root.
|
|
|
|
[h4 Testing]
|
|
|
|
A mixture of spot tests of values calculated using functions.wolfram.com,
|
|
and randomly generated test data are
|
|
used: the test data was computed using
|
|
[@http://shoup.net/ntl/doc/RR.txt NTL::RR] at 1000-bit precision.
|
|
|
|
[h4 Implementation]
|
|
|
|
These functions are implemented using the stable three term
|
|
recurrence relations. These relations guarantee low absolute error
|
|
but cannot guarantee low relative error near one of the roots of the
|
|
polynomials.
|
|
|
|
[endsect] [/section:beta_function The Beta Function]
|
|
[/
|
|
Copyright 2006 John Maddock and Paul A. Bristow.
|
|
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).
|
|
]
|