## `orthpoly::jacobi` -- the Jacobi polynomials

#### Introduction

`orthpoly::jacobi``(n,a,b,x)` computes the value of the n-th degree Jacobi polynomial with parameters a and b at the point x.

#### Call(s)

`orthpoly::jacobi(n, a, b, x)`

#### Parameters

`n` | - | a nonnegative integer: the degree of the polynomial.
`a, b` | - | arithmetical expressions.
`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 Jacobi polynomial at this point is returned as an arithmetical expression. If `n` is not a nonnegative integer, then `orthpoly::jacobi` returns itself symbolically.

#### Related Functions

`orthpoly::chebyshev1`, `orthpoly::chebyshev2`, `orthpoly::gegenbauer`, `orthpoly::legendre`

#### Details

- Evaluation for real floating point values x from the interval [-1.0,1.0] is numerically stable. Cf. example 2.

#### Example 1

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

    >> orthpoly::jacobi(2, a, b, x)

              / /                    2    2       \
              | | 7 a   7 b   a b   a    b        |  2
          poly| | --- + --- + --- + -- + -- + 3/2 | x  +
              \ \  8     8     4    8    8        /

             /              2    2 \     /  2                  2       \
             | 3 a   3 b   a    b  |     | a    b   a b   a   b        |
             | --- - --- + -- - -- | x + | -- - - - --- - - + -- - 1/2 |,
             \  4     4    4    4  /     \ 8    8    4    8   8        /

                 \
                 |
             [x] |
                 /

    >> orthpoly::jacobi(3, 4, 5, x[1])

                            3            2
             poly(455/8 x[1]  - 91/8 x[1]  - 91/8 x[1] + 7/8, [x[1]])

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

    >> orthpoly::jacobi(2, 4, b, 6*x)

          (b + 1) (6 x - 1)
          ----------------- +
                  4

                                   /     / b       \   b       \
             (6 x (b + 8) - b + 2) | 6 x | - + 7/2 | - - + 5/2 |
                                   \     \ 2       /   2       /
             ---------------------------------------------------
                                      4

    >> orthpoly::jacobi(2, 0, I, x[1] + 2)

          (1/4 + 1/4 I) x[1] + (((4 + I) x[1] + (6 + I))

             ((3/2 + 1/2 I) x[1] + (7/2 + 1/2 I))) / 4 + (1/4 + 1/4 I)

``Arithmetical expressions'' include numbers:

    >> orthpoly::jacobi(2, 1/2, -1/2, sqrt(2)),
       orthpoly::jacobi(3, 2, 5, 8 + I),
       orthpoly::jacobi(1000, 1, 2, 0.3);

             1/2   1/2
          3 2    (2    + 1/2)
          ------------------- - 3/8, 31733/2 + 12859/2 I, -0.06546648097
                   2

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

    >> orthpoly::jacobi(n, a, b, x), orthpoly::jacobi(1/2, 0, 1, 1)

           orthpoly::jacobi(n, a, b, x), orthpoly::jacobi(1/2, 0, 1, 1)

#### Example 2

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

    >> orthpoly::jacobi(100, 1/2, 3/2, 0.9)

                                   0.2560339406

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::jacobi(100, 1/2, 3/2, x):

    >> evalp(P100, x = 0.9)

                                  2.139740624e14

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.256005789994057173724575383078

    >> delete P100, DIGITS:

#### Background

- The Jacobi polynomials are given by the recursion formula 2*n*c(n)*c(2*n-2)*P(n,a,b,x) = c(2*n-1)*(c(2*n-2)*c(2*n)*x+a^2-b^2)*P(n-1,a,b,x) -2*(n-1+a)*(n-1+b)*c(2*n)*P(n-2,a,b,x) with c(i)= i+a+b and

      P(0,a,b,x)=1, P(1,a,b,x)=(a-b)/2+(1+(a+b)/2)*x.

- For fixed real a>-1, b>-1 the Jacobi polynomials are orthogonal on the interval [-1,1] with respect to the weight function w(x)=(1-x)^a*(1+x)^b.
- For special values of the parameters a, b the Jacobi polynomials are related to the Legendre polynomials

      P(n,x)=P(n,0,0,x),

  to the Chebyshev polynomials of the first kind

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

  to the Chebyshev polynomials of the second kind

      U(n,x)=2^(2*n)*n!*(n+1)!/(2*n+1)!*P(n,1/2,1/2,x),

  and to the Gegenbauer polynomials, respectively:

              G(n,a,x)= gamma(a+1/2)*gamma(n+2*a)/gamma(2*a)/gamma(n+a+1/2)*P(n,a-1/2,a-1/2,x).

#### Changes

- Indexed identifiers are now regarded as indeterminates. Now arbitrary expressions are accepted for the parameters a, b. Non-integer degrees n are now accepted, leading to an unevaluated function call. Efficiency was improved.

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