| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
| 65.1 Introduction to orthogonal polynomials | ||
| 65.2 Functions and Variables for orthogonal polynomials |
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
orthopoly is a package for symbolic and numerical evaluation of several
kinds of orthogonal polynomials, including Chebyshev, Laguerre, Hermite, Jacobi,
Legendre, and ultraspherical (Gegenbauer) polynomials. Additionally,
orthopoly includes support for the spherical Bessel, spherical Hankel,
and spherical harmonic functions.
For the most part, orthopoly follows the conventions of Abramowitz and
Stegun Handbook of Mathematical Functions, Chapter 22 (10th printing,
December 1972); additionally, we use Gradshteyn and Ryzhik,
Table of Integrals, Series, and Products (1980 corrected and enlarged
edition), and Eugen Merzbacher Quantum Mechanics (2nd edition, 1970).
Barton Willis of the University of Nebraska at Kearney (UNK) wrote the
orthopoly package and its documentation. The package is released under
the GNU General Public License (GPL).
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
load (orthopoly) loads the orthopoly package.
To find the third-order Legendre polynomial,
(%i1) legendre_p (3, x);
3 2
5 (1 - x) 15 (1 - x)
(%o1) - ---------- + ----------- - 6 (1 - x) + 1
2 2
To express this as a sum of powers of x, apply ratsimp
or
rat
to the result.
(%i2) [ratsimp (%), rat (%)];
3 3
5 x - 3 x 5 x - 3 x
(%o2)/R/ [----------, ----------]
2 2
Alternatively, make the second argument to legendre_p (its "main"
variable) a canonical rational expression (CRE).
(%i1) legendre_p (3, rat (x));
3
5 x - 3 x
(%o1)/R/ ----------
2
For floating point evaluation, orthopoly uses a running error analysis
to estimate an upper bound for the error. For example,
(%i1) jacobi_p (150, 2, 3, 0.2); (%o1) interval(- 0.062017037936715, 1.533267919277521E-11)
Intervals have the form interval (c, r), where c is the
center and r is the radius of the interval. Since Maxima does not support
arithmetic on intervals, in some situations, such as graphics, you want to
suppress the error and output only the center of the interval. To do this, set
the option variable orthopoly_returns_intervals
to false.
(%i1) orthopoly_returns_intervals : false; (%o1) false (%i2) jacobi_p (150, 2, 3, 0.2); (%o2) - 0.062017037936715
Refer to the section Floating point Evaluation for more information.
Most functions in orthopoly have a gradef
property; thus
(%i1) diff (hermite (n, x), x);
(%o1) 2 n H (x)
n - 1
(%i2) diff (gen_laguerre (n, a, x), x);
(a) (a)
n L (x) - (n + a) L (x) unit_step(n)
n n - 1
(%o2) ------------------------------------------
x
The unit step function in the second example prevents an error that would otherwise arise by evaluating with n equal to 0.
(%i3) ev (%, n = 0); (%o3) 0
The gradef property only applies to the "main" variable; derivatives
with respect other arguments usually result in an error message; for example
(%i1) diff (hermite (n, x), x);
(%o1) 2 n H (x)
n - 1
(%i2) diff (hermite (n, x), n);
Maxima doesn't know the derivative of hermite with respect the first
argument
-- an error. Quitting. To debug this try debugmode(true);
Generally, functions in orthopoly map over lists and matrices. For
the mapping to fully evaluate, the option variables doallmxops
and
listarith
must both be true (the defaults).
To illustrate the mapping over matrices, consider
(%i1) hermite (2, x);
2
(%o1) - 2 (1 - 2 x )
(%i2) m : matrix ([0, x], [y, 0]);
[ 0 x ]
(%o2) [ ]
[ y 0 ]
(%i3) hermite (2, m);
[ 2 ]
[ - 2 - 2 (1 - 2 x ) ]
(%o3) [ ]
[ 2 ]
[ - 2 (1 - 2 y ) - 2 ]
In the second example, the i, j element of the value
is hermite (2, m[i,j]); this is not the same as computing
-2 + 4 m . m, as seen in the next example.
(%i4) -2 * matrix ([1, 0], [0, 1]) + 4 * m . m;
[ 4 x y - 2 0 ]
(%o4) [ ]
[ 0 4 x y - 2 ]
If you evaluate a function at a point outside its domain, generally
orthopoly returns the function unevaluated. For example,
(%i1) legendre_p (2/3, x);
(%o1) P (x)
2/3
orthopoly supports translation into TeX; it also does two-dimensional
output on a terminal.
(%i1) spherical_harmonic (l, m, theta, phi);
m
(%o1) Y (theta, phi)
l
(%i2) tex (%);
$$Y_{l}^{m}\left(\vartheta,\varphi\right)$$
(%o2) false
(%i3) jacobi_p (n, a, a - b, x/2);
(a, a - b) x
(%o3) P (-)
n 2
(%i4) tex (%);
$$P_{n}^{\left(a,a-b\right)}\left({{x}\over{2}}\right)$$
(%o4) false
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
When an expression involves several orthogonal polynomials with symbolic orders, it's possible that the expression actually vanishes, yet Maxima is unable to simplify it to zero. If you divide by such a quantity, you'll be in trouble. For example, the following expression vanishes for integers n greater than 1, yet Maxima is unable to simplify it to zero.
(%i1) (2*n - 1) * legendre_p (n - 1, x) * x - n * legendre_p (n, x)
+ (1 - n) * legendre_p (n - 2, x);
(%o1) (2 n - 1) P (x) x - n P (x) + (1 - n) P (x)
n - 1 n n - 2
For a specific n, we can reduce the expression to zero.
(%i2) ev (% ,n = 10, ratsimp); (%o2) 0
Generally, the polynomial form of an orthogonal polynomial is ill-suited for floating point evaluation. Here's an example.
(%i1) p : jacobi_p (100, 2, 3, x)$ (%i2) subst (0.2, x, p); (%o2) 3.4442767023833592E+35 (%i3) jacobi_p (100, 2, 3, 0.2); (%o3) interval(0.18413609135169, 6.8990300925815987E-12) (%i4) float(jacobi_p (100, 2, 3, 2/10)); (%o4) 0.18413609135169
The true value is about 0.184; this calculation suffers from extreme subtractive cancellation error. Expanding the polynomial and then evaluating, gives a better result.
(%i5) p : expand(p)$ (%i6) subst (0.2, x, p); (%o6) 0.18413609766122982
This isn't a general rule; expanding the polynomial does not always result in an expression that is better suited for numerical evaluation. By far, the best way to do numerical evaluation is to make one or more of the function arguments floating point numbers. By doing that, specialized floating point algorithms are used for evaluation.
Maxima's float
function is somewhat indiscriminate; if you apply
float to an expression involving an orthogonal polynomial with a
symbolic degree or order parameter, these parameters may be
converted into floats; after that, the expression will not evaluate
fully. Consider
(%i1) assoc_legendre_p (n, 1, x);
1
(%o1) P (x)
n
(%i2) float (%);
1.0
(%o2) P (x)
n
(%i3) ev (%, n=2, x=0.9);
1.0
(%o3) P (0.9)
2
The expression in (%o3) will not evaluate to a float; orthopoly doesn't
recognize floating point values where it requires an integer. Similarly,
numerical evaluation of the pochhammer
function for orders that
exceed pochhammer_max_index
can be troublesome; consider
(%i1) x : pochhammer (1, 10), pochhammer_max_index : 5;
(%o1) (1)
10
Applying float doesn't evaluate x to a float
(%i2) float (x);
(%o2) (1.0)
10.0
To evaluate x to a float, you'll need to bind pochhammer_max_index
to 11 or greater and apply float to x.
(%i3) float (x), pochhammer_max_index : 11; (%o3) 3628800.0
The default value of pochhammer_max_index is 100;
change its value after loading orthopoly.
Finally, be aware that reference books vary on the definitions of the orthogonal polynomials; we've generally used the conventions of conventions of Abramowitz and Stegun.
Before you suspect a bug in orthopoly, check some special cases
to determine if your definitions match those used by orthopoly.
Definitions often differ by a normalization; occasionally, authors
use "shifted" versions of the functions that makes the family
orthogonal on an interval other than (-1, 1). To define, for example,
a Legendre polynomial that is orthogonal on (0, 1), define
(%i1) shifted_legendre_p (n, x) := legendre_p (n, 2*x - 1)$
(%i2) shifted_legendre_p (2, rat (x));
2
(%o2)/R/ 6 x - 6 x + 1
(%i3) legendre_p (2, rat (x));
2
3 x - 1
(%o3)/R/ --------
2
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
Most functions in orthopoly use a running error analysis to
estimate the error in floating point evaluation; the
exceptions are the spherical Bessel functions and the associated Legendre
polynomials of the second kind. For numerical evaluation, the spherical
Bessel functions call SLATEC functions. No specialized method is used
for numerical evaluation of the associated Legendre polynomials of the
second kind.
The running error analysis ignores errors that are second or higher order in the machine epsilon (also known as unit roundoff). It also ignores a few other errors. It's possible (although unlikely) that the actual error exceeds the estimate.
Intervals have the form interval (c, r), where c is
the center of the interval and r is its radius. The center of an interval
can be a complex number, and the radius is always a positive real number.
Here is an example.
(%i1) fpprec : 50$ (%i2) y0 : jacobi_p (100, 2, 3, 0.2); (%o2) interval(0.1841360913516871, 6.8990300925815987E-12) (%i3) y1 : bfloat (jacobi_p (100, 2, 3, 1/5)); (%o3) 1.8413609135168563091370224958913493690868904463668b-1
Let's test that the actual error is smaller than the error estimate
(%i4) is (abs (part (y0, 1) - y1) < part (y0, 2)); (%o4) true
Indeed, for this example the error estimate is an upper bound for the true error.
Maxima does not support arithmetic on intervals.
(%i1) legendre_p (7, 0.1) + legendre_p (8, 0.1);
(%o1) interval(0.18032072148437508, 3.1477135311021797E-15)
+ interval(- 0.19949294375000004, 3.3769353084291579E-15)
A user could define arithmetic operators that do interval math. To define interval addition, we can define
(%i1) infix ("@+")$
(%i2) "@+"(x,y) := interval (part (x, 1) + part (y, 1), part (x, 2)
+ part (y, 2))$
(%i3) legendre_p (7, 0.1) @+ legendre_p (8, 0.1);
(%o3) interval(- 0.019172222265624955, 6.5246488395313372E-15)
The special floating point routines get called when the arguments are complex. For example,
(%i1) legendre_p (10, 2 + 3.0*%i);
(%o1) interval(- 3.876378825E+7 %i - 6.0787748E+7,
1.2089173052721777E-6)
Let's compare this to the true value.
(%i1) float (expand (legendre_p (10, 2 + 3*%i))); (%o1) - 3.876378825E+7 %i - 6.0787748E+7
Additionally, when the arguments are big floats, the special floating point routines get called; however, the big floats are converted into double floats and the final result is a double.
(%i1) ultraspherical (150, 0.5b0, 0.9b0); (%o1) interval(- 0.043009481257265, 3.3750051301228864E-14)
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
orthopoly To plot expressions that involve the orthogonal polynomials, you must do two things:
orthopoly_returns_intervals
to false,
orthopoly functions.
If function calls aren't quoted, Maxima evaluates them to polynomials before plotting; consequently, the specialized floating point code doesn't get called. Here is an example of how to plot an expression that involves a Legendre polynomial.
(%i1) plot2d ('(legendre_p (5, x)), [x, 0, 1]),
orthopoly_returns_intervals : false;
(%o1)
The entire expression legendre_p (5, x) is quoted; this is different
than just quoting the function name using 'legendre_p (5, x).
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
The orthopoly package defines the Pochhammer symbol and a unit step
function. orthopoly uses the Kronecker delta function and the unit step
function in gradef
statements.
To convert Pochhammer symbols into quotients of gamma functions,
use makegamma.
(%i1) makegamma (pochhammer (x, n));
gamma(x + n)
(%o1) ------------
gamma(x)
(%i2) makegamma (pochhammer (1/2, 1/2));
1
(%o2) ---------
sqrt(%pi)
Derivatives of the Pochhammer symbol are given in terms of the psi
function.
(%i1) diff (pochhammer (x, n), x);
(%o1) (x) (psi (x + n) - psi (x))
n 0 0
(%i2) diff (pochhammer (x, n), n);
(%o2) (x) psi (x + n)
n 0
You need to be careful with the expression in (%o1); the difference of the
psi functions has polynomials when x = -1, -2, .., -n.
These polynomials cancel with factors in pochhammer (x, n)
making the derivative a degree n - 1 polynomial when n is a
positive integer.
The Pochhammer symbol is defined for negative orders through its representation as a quotient of gamma functions. Consider
(%i1) q : makegamma (pochhammer (x, n));
gamma(x + n)
(%o1) ------------
gamma(x)
(%i2) sublis ([x=11/3, n= -6], q);
729
(%o2) - ----
2240
Alternatively, we can get this result directly.
(%i1) pochhammer (11/3, -6);
729
(%o1) - ----
2240
The unit step function is left-continuous; thus
(%i1) [unit_step (-1/10), unit_step (0), unit_step (1/10)]; (%o1) [0, 0, 1]
If you need a unit step function that is neither left or right continuous
at zero, define your own using signum,
for example,
(%i1) xunit_step (x) := (1 + signum (x))/2$
(%i2) [xunit_step (-1/10), xunit_step (0), xunit_step (1/10)];
1
(%o2) [0, -, 1]
2
Do not redefine unit_step
itself; some code in orthopoly
requires that the unit step function be left-continuous.
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
Generally, orthopoly does symbolic evaluation by using a hypergeometic
representation of the orthogonal polynomials. The hypergeometic functions are
evaluated using the (undocumented) functions hypergeo11 and
hypergeo21. The exceptions are the half-integer Bessel functions and the
associated Legendre function of the second kind. The half-integer Bessel
functions are evaluated using an explicit representation, and the associated
Legendre function of the second kind is evaluated using recursion.
For floating point evaluation, we again convert most functions into a hypergeometic form; we evaluate the hypergeometic functions using forward recursion. Again, the exceptions are the half-integer Bessel functions and the associated Legendre function of the second kind. Numerically, the half-integer Bessel functions are evaluated using the SLATEC code.
| [ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
The associated Legendre function of the first kind of degree n and order m.
Reference: Abramowitz and Stegun, equations 22.5.37, page 779, 8.6.6 (second equation), page 334, and 8.2.5, page 333.
The associated Legendre function of the second kind of degree n and order m.
Reference: Abramowitz and Stegun, equation 8.5.3 and 8.1.8.
The Chebyshev function of the first kind.
Reference: Abramowitz and Stegun, equation 22.5.47, page 779.
The Chebyshev function of the second kind.
Reference: Abramowitz and Stegun, equation 22.5.48, page 779.
The generalized Laguerre polynomial of degree n.
Reference: Abramowitz and Stegun, equation 22.5.54, page 780.
The Hermite polynomial.
Reference: Abramowitz and Stegun, equation 22.5.55, page 780.
Return true if the input is an interval and return false if it isn't.
The Jacobi polynomial.
The Jacobi polynomials are actually defined for all a and b;
however, the Jacobi polynomial weight
(1 - x)^a (1 + x)^b isn't integrable for
a <= -1 or b <= -1.
Reference: Abramowitz and Stegun, equation 22.5.42, page 779.
The Laguerre polynomial.
Reference: Abramowitz and Stegun, equations 22.5.16 and 22.5.54, page 780.
The Legendre polynomial of the first kind.
Reference: Abramowitz and Stegun, equations 22.5.50 and 22.5.51, page 779.
The Legendre polynomial of the first kind.
Reference: Abramowitz and Stegun, equations 8.5.3 and 8.1.8.
Returns a recursion relation for the orthogonal function family f with arguments args. The recursion is with respect to the polynomial degree.
(%i1) orthopoly_recur (legendre_p, [n, x]);
(2 n - 1) P (x) x + (1 - n) P (x)
n - 1 n - 2
(%o1) P (x) = -----------------------------------------
n n
The second argument to orthopoly_recur must be a list with the
correct number of arguments for the function f; if it isn't,
Maxima signals an error.
(%i1) orthopoly_recur (jacobi_p, [n, x]); Function jacobi_p needs 4 arguments, instead it received 2 -- an error. Quitting. To debug this try debugmode(true);
Additionally, when f isn't the name of one of the families of orthogonal polynomials, an error is signalled.
(%i1) orthopoly_recur (foo, [n, x]); A recursion relation for foo isn't known to Maxima -- an error. Quitting. To debug this try debugmode(true);
Default value: true
When orthopoly_returns_intervals is true, floating point results
are returned in the form interval (c, r), where c is
the center of an interval and r is its radius. The center can be a
complex number; in that case, the interval is a disk in the complex plane.
Returns a three element list; the first element is the formula of the weight for the orthogonal polynomial family f with arguments given by the list args; the second and third elements give the lower and upper endpoints of the interval of orthogonality. For example,
(%i1) w : orthopoly_weight (hermite, [n, x]);
2
- x
(%o1) [%e , - inf, inf]
(%i2) integrate(w[1]*hermite(3, x)*hermite(2, x), x, w[2], w[3]);
(%o2) 0
The main variable of f must be a symbol; if it isn't, Maxima signals an error.
The Pochhammer symbol. For nonnegative integers n with
n <= pochhammer_max_index, the expression
pochhammer (x, n) evaluates to the product
x (x + 1) (x + 2) ... (x + n - 1) when
n > 0 and to 1 when n = 0. For negative n,
pochhammer (x, n) is defined as
(-1)^n / pochhammer (1 - x, -n). Thus
(%i1) pochhammer (x, 3);
(%o1) x (x + 1) (x + 2)
(%i2) pochhammer (x, -3);
1
(%o2) - -----------------------
(1 - x) (2 - x) (3 - x)
To convert a Pochhammer symbol into a quotient of gamma functions,
(see Abramowitz and Stegun, equation 6.1.22) use makegamma,
for
example
(%i1) makegamma (pochhammer (x, n));
gamma(x + n)
(%o1) ------------
gamma(x)
When n exceeds pochhammer_max_index
or when n is symbolic,
pochhammer returns a noun form.
(%i1) pochhammer (x, n);
(%o1) (x)
n
Default value: 100
pochhammer (n, x) expands to a product if and only if
n <= pochhammer_max_index.
Examples:
(%i1) pochhammer (x, 3), pochhammer_max_index : 3;
(%o1) x (x + 1) (x + 2)
(%i2) pochhammer (x, 4), pochhammer_max_index : 3;
(%o2) (x)
4
Reference: Abramowitz and Stegun, equation 6.1.16, page 256.
The spherical Bessel function of the first kind.
Reference: Abramowitz and Stegun, equations 10.1.8, page 437 and 10.1.15, page 439.
The spherical Bessel function of the second kind.
Reference: Abramowitz and Stegun, equations 10.1.9, page 437 and 10.1.15, page 439.
The spherical Hankel function of the first kind.
Reference: Abramowitz and Stegun, equation 10.1.36, page 439.
The spherical Hankel function of the second kind.
Reference: Abramowitz and Stegun, equation 10.1.17, page 439.
The spherical harmonic function.
Reference: Merzbacher 9.64.
The left-continuous unit step function; thus unit_step (x) vanishes
for x <= 0 and equals 1 for x > 0.
If you want a unit step function that takes on the value 1/2 at zero, use
(1 + signum (x))/2.
The ultraspherical polynomial (also known as the Gegenbauer polynomial).
Reference: Abramowitz and Stegun, equation 22.5.46, page 779.
| [ << ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
This document was generated by Crategus on Dezember, 12 2012 using texi2html 1.76.