## `orthpoly::legendre` -- the Legendre polynomials

#### Introduction

`orthpoly::legendre``(n,x)` computes the value of the n-th degree Legendre polynomial at the point x.

#### Call(s)

`orthpoly::legendre(n, x)`

#### Parameters

`n` | - | a nonnegative integer: the degree of the polynomial.
`x` | - | an indeterminate or an arithmetical expression. An indeterminate is either an identifier (of domain type `DOM_IDENT`) or an indexed identifier (of type `"_index"`).

#### Returns

If `x` is an indeterminate, then a polynomial of domain type `DOM_POLY` is returned. If `x` is an arithmetical expression, then the value of the Legendre polynomial at this point is returned as an arithmetical expression. If `n` is not a nonnegative integer, then `orthpoly::legendre` returns itself symbolically.

#### Related Functions

`numeric::gldata`, `orthpoly::gegenbauer`, `orthpoly::jacobi`

#### Details

- These polynomials have rational coefficients.
- Evaluation for real floating point values x from the interval [-1.0,1.0] is numerically stable. Cf. example 2.
- Use `numeric::gldata` to compute the roots of the Legendre polynomials. Cf. example 3.

#### Example 1

Polynomials of domain type `DOM_POLY` are returned, if identifiers or indexed identifiers are specified:

    >> orthpoly::legendre(2, x)

                                        2
                              poly(3/2 x  - 1/2, [x])

    >> orthpoly::legendre(3, x[1])

                                     3
                        poly(5/2 x[1]  - 3/2 x[1], [x[1]])

However, using arithmetical expressions as input the ``values'' of these polynomials are returned:

    >> orthpoly::legendre(2, 6*x)

                                        2
                                    54 x  - 1/2

    >> orthpoly::legendre(3, x[1] + 2)

                             /             2       \
                             | 3 (x[1] + 2)        |
                5 (x[1] + 2) | ------------- - 1/2 |
                             \       2             /   2 x[1]
                ------------------------------------ - ------ - 4/3
                                 3                       3

``Arithmetical expressions'' include numbers:

    >> orthpoly::legendre(2, sqrt(2)), orthpoly::legendre(3, 8 + I),
       orthpoly::legendre(1000, 0.3)

                         5/2, 1208 + 476 I, -0.0256691675

If no integer degree is specified, then `orthpoly::legendre` returns itself symbolically:

    >> orthpoly::legendre(n, x), orthpoly::legendre(1/2, x)

               orthpoly::legendre(n, x), orthpoly::legendre(1/2, x)

#### Example 2

If a floating point value is desired, then a direct call such as

    >> orthpoly::legendre(100, 0.9)

                                   0.1022658206

is appropriate and yields a correct result. One should not evaluate the symbolic polynomial at a floating point value, because this may be numerically unstable:

    >> P100 := orthpoly::legendre(100, x):

    >> evalp(P100, x = 0.9)

                                  4.222688586e13

This result is caused by numerical round-off. Also with increased `DIGITS` only a few leading digits are correct:

    >> DIGITS := 30: evalp(P100, x = 0.9)

                         0.102266645970798303097499109014

    >> delete P100, DIGITS:

#### Example 3

We recommend to use `numeric::gldata` for computing roots of the Legendre polynomial P(n,x). This routine provides all roots of the function Q(n,y)=P(n,2*y-1):

    >> QRoots := numeric::gldata(5, DIGITS)[2]

           [0.04691007703, 0.2307653449, 1/2, 0.769234655, 0.9530899229]

These values are easily transformed to roots of P(n,x):

    >> PRoots := map(QRoots, y -> 2*y - 1)

           [-0.9061798459, -0.5384693101, 0, 0.5384693101, 0.9061798459]

    >> orthpoly::legendre(5, r) $ r in PRoots

          -1.08046818e-14, -1.084202173e-19, 0, 1.084202173e-19,

             1.08046818e-14

    >> delete QRoots, PRoots:

#### Background

- The Legendre polynomials are given by P(n,x)=1/2^n/n!*diff((x^2-1)^n, x$n)..
- They satisfy the recursion formula

      P(n,x)=(2*n-1)/n*x*P(n-1,x)-(n-1)/n*P(n-2,x)

  with P(0,x)=1 and P(1,x)=x.
- They are orthogonal on the interval [-1,1] with respect to the weight function w(x)=1.
- P(n,x) coincides with the Gegenbauer polynomial P(n,1/2,x).
- P(n,x) coincides with the Jacobi polynomial P(n,0,0,x).

#### Changes

- Indexed identifiers are now regarded as indeterminates. Non-integer degrees n are now accepted, leading to an unevaluated function call. Floating point evaluations are now numerically stable. Efficiency was improved.

Do you have questions or comments?
Copyright © SciFace Software GmbH & Co. KG 2000