Orthogonal Polynomials#

Chebyshev polynomials#

The Chebyshev polynomial of the first kind arises as a solution to the differential equation

(1x2)yxy+n2y=0

and those of the second kind as a solution to

(1x2)y3xy+n(n+2)y=0.

The Chebyshev polynomials of the first kind are defined by the recurrence relation

T0(x)=1,T1(x)=x,Tn+1(x)=2xTn(x)Tn1(x).

The Chebyshev polynomials of the second kind are defined by the recurrence relation

U0(x)=1,U1(x)=2x,Un+1(x)=2xUn(x)Un1(x).

For integers m,n, they satisfy the orthogonality relations

11Tn(x)Tm(x)dx1x2={0if nm,πif n=m=0,π/2if n=m0,

and

11Un(x)Um(x)1x2dx=π2δm,n.

They are named after Pafnuty Chebyshev (1821-1894, alternative transliterations: Tchebyshef or Tschebyscheff).

Hermite polynomials#

The Hermite polynomials are defined either by

Hn(x)=(1)nex2/2dndxnex2/2

(the “probabilists’ Hermite polynomials”), or by

Hn(x)=(1)nex2dndxnex2

(the “physicists’ Hermite polynomials”). Sage (via Maxima) implements the latter flavor. These satisfy the orthogonality relation

Hn(x)Hm(x)ex2dx=πn!2nδnm.

They are named in honor of Charles Hermite (1822-1901), but were first introduced by Laplace in 1810 and also studied by Chebyshev in 1859.

Legendre polynomials#

Each Legendre polynomial Pn(x) is an n-th degree polynomial. It may be expressed using Rodrigues’ formula:

Pn(x)=(2nn!)1dndxn[(x21)n].

These are solutions to Legendre’s differential equation:

ddx[(1x2)ddxP(x)]+n(n+1)P(x)=0

and satisfy the orthogonality relation

11Pm(x)Pn(x)dx=22n+1δmn.

The Legendre function of the second kind Qn(x) is another (linearly independent) solution to the Legendre differential equation. It is not an “orthogonal polynomial” however.

The associated Legendre functions of the first kind Pm(x) can be given in terms of the “usual” Legendre polynomials by

Pm(x)=(1)m(1x2)m/2dmdxmP(x)=(1)m2!(1x2)m/2d+mdx+m(x21).

Assuming 0m, they satisfy the orthogonality relation:

11Pk(m)P(m)dx=2(+m)!(2+1)(m)! δk,,

where δk, is the Kronecker delta.

The associated Legendre functions of the second kind Qm(x) can be given in terms of the “usual” Legendre polynomials by

Qm(x)=(1)m(1x2)m/2dmdxmQ(x).

They are named after Adrien-Marie Legendre (1752-1833).

Laguerre polynomials#

Laguerre polynomials may be defined by the Rodrigues formula

Ln(x)=exn!dndxn(exxn).

They are solutions of Laguerre’s equation:

xy+(1x)y+ny=0

and satisfy the orthogonality relation

0Lm(x)Ln(x)exdx=δmn.

The generalized Laguerre polynomials may be defined by the Rodrigues formula:

Ln(α)(x)=xαexn!dndxn(exxn+α).

(These are also sometimes called the associated Laguerre polynomials.) The simple Laguerre polynomials are recovered from the generalized polynomials by setting α=0.

They are named after Edmond Laguerre (1834-1886).

Jacobi polynomials#

Jacobi polynomials are a class of orthogonal polynomials. They are obtained from hypergeometric series in cases where the series is in fact finite:

Pn(α,β)(z)=(α+1)nn!2F1(n,1+α+β+n;α+1;1z2),

where ()n is Pochhammer’s symbol (for the rising factorial), (Abramowitz and Stegun p561.) and thus have the explicit expression

Pn(α,β)(z)=Γ(α+n+1)n!Γ(α+β+n+1)m=0n(nm)Γ(α+β+n+m+1)Γ(α+m+1)(z12)m.

They are named after Carl Gustav Jaboc Jacobi (1804-1851).

Gegenbauer polynomials#

Ultraspherical or Gegenbauer polynomials are given in terms of the Jacobi polynomials Pn(α,β)(x) with α=β=a1/2 by

Cn(a)(x)=Γ(a+1/2)Γ(2a)Γ(n+2a)Γ(n+a+1/2)Pn(a1/2,a1/2)(x).

They satisfy the orthogonality relation

11(1x2)a1/2Cm(a)(x)Cn(a)(x)dx=δmn212aπΓ(n+2a)(n+a)Γ2(a)Γ(n+1),

for a>1/2. They are obtained from hypergeometric series in cases where the series is in fact finite:

Cn(a)(z)=(2a)nn!2F1(n,2a+n;a+12;1z2)

where n is the falling factorial. (See Abramowitz and Stegun p561.)

They are named for Leopold Gegenbauer (1849-1903).

Krawtchouk polynomials#

The Krawtchouk polynomials are discrete orthogonal polynomials that are given by the hypergeometric series

Kj(x;n,p)=(1)j(nj)pj2F1(j,x;n;p1).

Since they are discrete orthogonal polynomials, they satisfy an orthogonality relation defined on a discrete (in this case finite) set of points:

m=0nKi(m;n,p)Kj(m;n,p)(nm)pmqnm=(nj)(pq)jδij,

where q=1p. They can also be described by the recurrence relation

jKj(x;n,p)=(x(nj+1)p(j1)q)Kj1(x;n,p)pq(nj+2)Kj2(x;n,p),

where K0(x;n,p)=1 and K1(x;n,p)=xnp.

They are named for Mykhailo Krawtchouk (1892-1942).

Meixner polynomials#

The Meixner polynomials are discrete orthogonal polynomials that are given by the hypergeometric series

Mn(x;n,p)=(1)j(nj)pj2F1(j,x;n;p1).

They satisfy an orthogonality relation:

k=0M~n(k;b,c)M~m(k;b,c)(b)kk!ck=cnn!(b)n(1c)bδmn,

where M~n(x;b,c)=Mn(x;b,c)/(b)x, for b>0. They can also be described by the recurrence relation

c(n1+b)Mn(x;b,c)=((c1)x+n1+c(n1+b))(b+n1)Mn1(x;b,c)(b+n1)(b+n2)(n1)Mn2(x;b,c),

where M0(x;b,c)=0 and M1(x;b,c)=(1c1)x+b.

They are named for Josef Meixner (1908-1994).

Hahn polynomials#

The Hahn polynomials are discrete orthogonal polynomials that are given by the hypergeometric series

Qk(x;a,b,n)=3F2(k,k+a+b+1,x;a+1,n;1).

They satisfy an orthogonality relation:

k=0n1Qi(k;a,b,n)Qj(k;a,b,n)ρ(k)=δijπi,

where

ρ(k)=(a+kk)(b+nknk),πi=δij(1)ii!(b+1)i(i+a+b+1)n+1n!(2i+a+b+1)(n)i(a+1)i.

They can also be described by the recurrence relation

AQk(x;a,b,n)=(x+A+C)Qk1(x;a,b,n)CQk2(x;a,b,n),

where Q0(x;a,b,n)=1 and Q1(x;a,b,n)=1a+b+2(a+1)nx and

A=(k+a+b)(k+a)(nk+1)(2k+a+b1)(2k+a+b),C=(k1)(k+b1)(k+a+b+n)(2k+a+b2)(2k+a+b1).

They are named for Wolfgang Hahn (1911-1998), although they were first introduced by Chebyshev in 1875.

Pochhammer symbol#

For completeness, the Pochhammer symbol, introduced by Leo August Pochhammer, (x)n, is used in the theory of special functions to represent the “rising factorial” or “upper factorial”

(x)n=x(x+1)(x+2)(x+n1)=(x+n1)!(x1)!.

On the other hand, the falling factorial or lower factorial is

xn=x!(xn)!,

in the notation of Ronald L. Graham, Donald E. Knuth and Oren Patashnik in their book Concrete Mathematics.

Todo

Implement Zernike polynomials. Wikipedia article Zernike_polynomials

REFERENCES:

AUTHORS:

  • David Joyner (2006-06)

  • Stefan Reiterer (2010-)

  • Ralf Stephan (2015-)

The original module wrapped some of the orthogonal/special functions in the Maxima package “orthopoly” and was written by Barton Willis of the University of Nebraska at Kearney.

class sage.functions.orthogonal_polys.ChebyshevFunction(name, nargs=2, latex_name=None, conversions=None)#

Bases: sage.functions.orthogonal_polys.OrthogonalFunction

Abstract base class for Chebyshev polynomials of the first and second kind.

EXAMPLES:

sage: chebyshev_T(3,x)
4*x^3 - 3*x
class sage.functions.orthogonal_polys.Func_assoc_legendre_P#

Bases: sage.symbolic.function.BuiltinFunction

Return the Ferrers function Pnm(x) of first kind for x(1,1) with general order m and general degree n.

Ferrers functions of first kind are one of two linearly independent solutions of the associated Legendre differential equation

(1x2)d2wdx22xdwdx+(n(n+1)m21x2)w=0

on the interval x(1,1) and are usually denoted by Pnm(x).

See also

The other linearly independent solution is called Ferrers function of second kind and denoted by Qnm(x), see Func_assoc_legendre_Q.

Warning

Ferrers functions must be carefully distinguished from associated Legendre functions which are defined on C(,1] and have not yet been implemented.

EXAMPLES:

We give the first Ferrers functions for non-negative integers n and m in the interval 1<x<1:

sage: for n in range(4):
....:     for m in range(n+1):
....:         print(f"P_{n}^{m}({x}) = {gen_legendre_P(n, m, x)}")
P_0^0(x) = 1
P_1^0(x) = x
P_1^1(x) = -sqrt(-x^2 + 1)
P_2^0(x) = 3/2*x^2 - 1/2
P_2^1(x) = -3*sqrt(-x^2 + 1)*x
P_2^2(x) = -3*x^2 + 3
P_3^0(x) = 5/2*x^3 - 3/2*x
P_3^1(x) = -3/2*(5*x^2 - 1)*sqrt(-x^2 + 1)
P_3^2(x) = -15*(x^2 - 1)*x
P_3^3(x) = -15*(-x^2 + 1)^(3/2)

These expressions for non-negative integers are computed by the Rodrigues-type given in eval_gen_poly(). Negative values for n are obtained by the following identity:

Pnm(x)=Pn1m(x).

For n being a non-negative integer, negative values for m are obtained by

Pn|m|(x)=(1)|m|(n|m|)!(n+|m|)!Pn|m|(x),

where |m|n.

Here are some specific values with negative integers:

sage: gen_legendre_P(-2, -1, x)
1/2*sqrt(-x^2 + 1)
sage: gen_legendre_P(2, -2, x)
-1/8*x^2 + 1/8
sage: gen_legendre_P(3, -2, x)
-1/8*(x^2 - 1)*x
sage: gen_legendre_P(1, -2, x)
0

Here are some other random values with floating numbers:

sage: m = var('m'); assume(m, 'integer')
sage: gen_legendre_P(m, m, .2)
0.960000000000000^(1/2*m)*(-1)^m*factorial(2*m)/(2^m*factorial(m))
sage: gen_legendre_P(.2, m, 0)
sqrt(pi)*2^m/(gamma(-1/2*m + 1.10000000000000)*gamma(-1/2*m + 0.400000000000000))
sage: gen_legendre_P(.2, .2, .2)
0.757714892929573

REFERENCES:

deprecated_function_alias(trac_number, func)#

Create an aliased version of a function or a method which raises a deprecation warning message.

If f is a function or a method, write g = deprecated_function_alias(trac_number, f) to make a deprecated aliased version of f.

INPUT:

  • trac_number – integer. The trac ticket number where the deprecation is introduced.

  • func – the function or method to be aliased

EXAMPLES:

sage: from sage.misc.superseded import deprecated_function_alias
sage: g = deprecated_function_alias(13109, number_of_partitions)
sage: g(5)
doctest:...: DeprecationWarning: g is deprecated. Please use sage.combinat.partition.number_of_partitions instead.
See http://trac.sagemath.org/13109 for details.
7

This also works for methods:

sage: class cls():
....:    def new_meth(self): return 42
....:    old_meth = deprecated_function_alias(13109, new_meth)
sage: cls().old_meth()
doctest:...: DeprecationWarning: old_meth is deprecated. Please use new_meth instead.
See http://trac.sagemath.org/13109 for details.
42

trac ticket #11585:

sage: def a(): pass
sage: b = deprecated_function_alias(13109, a)
sage: b()
doctest:...: DeprecationWarning: b is deprecated. Please use a instead.
See http://trac.sagemath.org/13109 for details.

AUTHORS:

  • Florent Hivert (2009-11-23), with the help of Mike Hansen.

  • Luca De Feo (2011-07-11), printing the full module path when different from old path

eval_gen_poly(n, m, arg, **kwds)#

Return the Ferrers function of first kind Pnm(x) for integers n>1,m>1 given by the following Rodrigues-type formula:

Pnm(x)=(1)m+n(1x2)m/22nn!dm+ndxm+n(1x2)n.

INPUT:

  • n – an integer degree

  • m – an integer order

  • x – either an integer or a non-numerical symbolic expression

EXAMPLES:

sage: gen_legendre_P(7,4,x)
3465/2*(13*x^3 - 3*x)*(x^2 - 1)^2
sage: gen_legendre_P(3,1,sqrt(x))
-3/2*(5*x - 1)*sqrt(-x + 1)

REFERENCE:

eval_poly(*args, **kwds)#

Deprecated: Use eval_gen_poly() instead. See trac ticket #25034 for details.

class sage.functions.orthogonal_polys.Func_assoc_legendre_Q#

Bases: sage.symbolic.function.BuiltinFunction

EXAMPLES:

sage: loads(dumps(gen_legendre_Q))
gen_legendre_Q
sage: maxima(gen_legendre_Q(2,1,3, hold=True))._sage_().simplify_full()
1/4*sqrt(2)*(36*pi - 36*I*log(2) + 25*I)
eval_recursive(n, m, x, **kwds)#

Return the associated Legendre Q(n, m, arg) function for integers n>1,m>1.

EXAMPLES:

sage: gen_legendre_Q(3,4,x)
48/(x^2 - 1)^2
sage: gen_legendre_Q(4,5,x)
-384/((x^2 - 1)^2*sqrt(-x^2 + 1))
sage: gen_legendre_Q(0,1,x)
-1/sqrt(-x^2 + 1)
sage: gen_legendre_Q(0,2,x)
-1/2*((x + 1)^2 - (x - 1)^2)/(x^2 - 1)
sage: gen_legendre_Q(2,2,x).subs(x=2).expand()
9/2*I*pi - 9/2*log(3) + 14/3
class sage.functions.orthogonal_polys.Func_chebyshev_T#

Bases: sage.functions.orthogonal_polys.ChebyshevFunction

Chebyshev polynomials of the first kind.

REFERENCE:

  • [AS1964] 22.5.31 page 778 and 6.1.22 page 256.

EXAMPLES:

sage: chebyshev_T(5,x)
16*x^5 - 20*x^3 + 5*x
sage: var('k')
k
sage: test = chebyshev_T(k,x)
sage: test
chebyshev_T(k, x)
eval_algebraic(n, x)#

Evaluate chebyshev_T as polynomial, using a recursive formula.

INPUT:

  • n – an integer

  • x – a value to evaluate the polynomial at (this can be any ring element)

EXAMPLES:

sage: chebyshev_T.eval_algebraic(5, x)
2*(2*(2*x^2 - 1)*x - x)*(2*x^2 - 1) - x
sage: chebyshev_T(-7, x) - chebyshev_T(7,x)
0
sage: R.<t> = ZZ[]
sage: chebyshev_T.eval_algebraic(-1, t)
t
sage: chebyshev_T.eval_algebraic(0, t)
1
sage: chebyshev_T.eval_algebraic(1, t)
t
sage: chebyshev_T(7^100, 1/2)
1/2
sage: chebyshev_T(7^100, Mod(2,3))
2
sage: n = 97; x = RIF(pi/2/n)
sage: chebyshev_T(n, cos(x)).contains_zero()
True
sage: R.<t> = Zp(2, 8, 'capped-abs')[]
sage: chebyshev_T(10^6+1, t)
(2^7 + O(2^8))*t^5 + O(2^8)*t^4 + (2^6 + O(2^8))*t^3 + O(2^8)*t^2 + (1 + 2^6 + O(2^8))*t + O(2^8)
eval_formula(n, x)#

Evaluate chebyshev_T using an explicit formula. See [AS1964] 227 (p. 782) for details for the recursions. See also [Koe1999] for fast evaluation techniques.

INPUT:

  • n – an integer

  • x – a value to evaluate the polynomial at (this can be any ring element)

EXAMPLES:

sage: chebyshev_T.eval_formula(-1,x)
x
sage: chebyshev_T.eval_formula(0,x)
1
sage: chebyshev_T.eval_formula(1,x)
x
sage: chebyshev_T.eval_formula(2,0.1) == chebyshev_T._evalf_(2,0.1)
True
sage: chebyshev_T.eval_formula(10,x)
512*x^10 - 1280*x^8 + 1120*x^6 - 400*x^4 + 50*x^2 - 1
sage: chebyshev_T.eval_algebraic(10,x).expand()
512*x^10 - 1280*x^8 + 1120*x^6 - 400*x^4 + 50*x^2 - 1
class sage.functions.orthogonal_polys.Func_chebyshev_U#

Bases: sage.functions.orthogonal_polys.ChebyshevFunction

Class for the Chebyshev polynomial of the second kind.

REFERENCE:

  • [AS1964] 22.8.3 page 783 and 6.1.22 page 256.

EXAMPLES:

sage: R.<t> = QQ[]
sage: chebyshev_U(2,t)
4*t^2 - 1
sage: chebyshev_U(3,t)
8*t^3 - 4*t
eval_algebraic(n, x)#

Evaluate chebyshev_U as polynomial, using a recursive formula.

INPUT:

  • n – an integer

  • x – a value to evaluate the polynomial at (this can be any ring element)

EXAMPLES:

sage: chebyshev_U.eval_algebraic(5,x)
-2*((2*x + 1)*(2*x - 1)*x - 4*(2*x^2 - 1)*x)*(2*x + 1)*(2*x - 1)
sage: parent(chebyshev_U(3, Mod(8,9)))
Ring of integers modulo 9
sage: parent(chebyshev_U(3, Mod(1,9)))
Ring of integers modulo 9
sage: chebyshev_U(-3,x) + chebyshev_U(1,x)
0
sage: chebyshev_U(-1,Mod(5,8))
0
sage: parent(chebyshev_U(-1,Mod(5,8)))
Ring of integers modulo 8
sage: R.<t> = ZZ[]
sage: chebyshev_U.eval_algebraic(-2, t)
-1
sage: chebyshev_U.eval_algebraic(-1, t)
0
sage: chebyshev_U.eval_algebraic(0, t)
1
sage: chebyshev_U.eval_algebraic(1, t)
2*t
sage: n = 97; x = RIF(pi/n)
sage: chebyshev_U(n-1, cos(x)).contains_zero()
True
sage: R.<t> = Zp(2, 6, 'capped-abs')[]
sage: chebyshev_U(10^6+1, t)
(2 + O(2^6))*t + O(2^6)
eval_formula(n, x)#

Evaluate chebyshev_U using an explicit formula.

See [AS1964] 227 (p. 782) for details on the recursions. See also [Koe1999] for the recursion formulas.

INPUT:

  • n – an integer

  • x – a value to evaluate the polynomial at (this can be any ring element)

EXAMPLES:

sage: chebyshev_U.eval_formula(10, x)
1024*x^10 - 2304*x^8 + 1792*x^6 - 560*x^4 + 60*x^2 - 1
sage: chebyshev_U.eval_formula(-2, x)
-1
sage: chebyshev_U.eval_formula(-1, x)
0
sage: chebyshev_U.eval_formula(0, x)
1
sage: chebyshev_U.eval_formula(1, x)
2*x
sage: chebyshev_U.eval_formula(2,0.1) == chebyshev_U._evalf_(2,0.1)
True
class sage.functions.orthogonal_polys.Func_gen_laguerre#

Bases: sage.functions.orthogonal_polys.OrthogonalFunction

REFERENCE:

  • [AS1964] 22.5.16, page 778 and page 789.

class sage.functions.orthogonal_polys.Func_hahn#

Bases: sage.functions.orthogonal_polys.OrthogonalFunction

Hahn polynomials Qk(x;a,b,n).

INPUT:

  • k – the degree

  • x – the independent variable x

  • a, b – the parameters a,b

  • n – the number of discrete points

EXAMPLES:

We verify the orthogonality for n=3:

sage: n = 2
sage: a,b = SR.var('a,b')
sage: def rho(k,a,b,n):
....:     return binomial(a+k,k) * binomial(b+n-k,n-k)
sage: M = matrix([[sum(rho(k,a,b,n)
....:                  * hahn(i,k,a,b,n) * hahn(j,k,a,b,n)
....:                  for k in range(n+1)).expand().factor()
....:              for i in range(n+1)] for j in range(n+1)])
sage: M = M.factor()
sage: P = rising_factorial
sage: def diag(i,a,b,n):
....:    return ((-1)^i * factorial(i) * P(b+1,i) * P(i+a+b+1,n+1)
....:            / (factorial(n) * (2*i+a+b+1) * P(-n,i) * P(a+1,i)))
sage: all(M[i,i] == diag(i,a,b,n) for i in range(3))
True
sage: all(M[i,j] == 0 for i in range(3) for j in range(3) if i != j)
True
eval_formula(k, x, a, b, n)#

Evaluate self using an explicit formula.

EXAMPLES:

sage: k,x,a,b,n = var('k,x,a,b,n')
sage: Q2 = hahn.eval_formula(2, x, a, b, n).simplify_full()
sage: Q2.coefficient(x^2).factor()
(a + b + 4)*(a + b + 3)/((a + 2)*(a + 1)*(n - 1)*n)
sage: Q2.coefficient(x).factor()
-(2*a*n - a + b + 4*n)*(a + b + 3)/((a + 2)*(a + 1)*(n - 1)*n)
sage: Q2(x=0)
1
eval_recursive(k, x, a, b, n, *args, **kwds)#

Return the Hahn polynomial Qk(x;a,b,n) using the recursive formula.

EXAMPLES:

sage: x,a,b,n = var('x,a,b,n')
sage: hahn.eval_recursive(0,x,a,b,n)
1
sage: hahn.eval_recursive(1,x,a,b,n)
-(a + b + 2)*x/((a + 1)*n) + 1
sage: bool(hahn(2,x,a,b,n) == hahn.eval_recursive(2,x,a,b,n))
True
sage: bool(hahn(3,x,a,b,n) == hahn.eval_recursive(3,x,a,b,n))
True
sage: bool(hahn(4,x,a,b,n) == hahn.eval_recursive(4,x,a,b,n))
True
sage: M = matrix([[-1/2,-1],[1,0]])
sage: ret = hahn.eval_recursive(2, M, 1, 2, n).simplify_full().factor()
sage: ret
[1/4*(4*n^2 + 8*n - 19)/((n - 1)*n)          3/2*(4*n + 3)/((n - 1)*n)]
[        -3/2*(4*n + 3)/((n - 1)*n)          (n^2 - n - 7)/((n - 1)*n)]
class sage.functions.orthogonal_polys.Func_hermite#

Bases: sage.symbolic.function.GinacFunction

Return the Hermite polynomial for integers n>1.

REFERENCE:

  • [AS1964] 22.5.40 and 22.5.41, page 779.

EXAMPLES:

sage: x = PolynomialRing(QQ, 'x').gen()
sage: hermite(2,x)
4*x^2 - 2
sage: hermite(3,x)
8*x^3 - 12*x
sage: hermite(3,2)
40
sage: S.<y> = PolynomialRing(RR)
sage: hermite(3,y)
8.00000000000000*y^3 - 12.0000000000000*y
sage: R.<x,y> = QQ[]
sage: hermite(3,y^2)
8*y^6 - 12*y^2
sage: w = var('w')
sage: hermite(3,2*w)
64*w^3 - 24*w
sage: hermite(5,3.1416)
5208.69733891963
sage: hermite(5,RealField(100)(pi))
5208.6167627118104649470287166

Check that trac ticket #17192 is fixed:

sage: x = PolynomialRing(QQ, 'x').gen()
sage: hermite(0,x)
1

sage: hermite(-1,x)
Traceback (most recent call last):
...
RuntimeError: hermite_eval: The index n must be a nonnegative integer

sage: hermite(-7,x)
Traceback (most recent call last):
...
RuntimeError: hermite_eval: The index n must be a nonnegative integer

sage: m,x = SR.var('m,x')
sage: hermite(m, x).diff(m)
Traceback (most recent call last):
...
RuntimeError: derivative w.r.t. to the index is not supported yet
class sage.functions.orthogonal_polys.Func_jacobi_P#

Bases: sage.functions.orthogonal_polys.OrthogonalFunction

Return the Jacobi polynomial Pn(a,b)(x) for integers n>1 and a and b symbolic or a>1 and b>1.

The Jacobi polynomials are actually defined for all a and b. However, the Jacobi polynomial weight (1x)a(1+x)b is not integrable for a1 or b1.

REFERENCE:

EXAMPLES:

sage: x = PolynomialRing(QQ, 'x').gen()
sage: jacobi_P(2,0,0,x)
3/2*x^2 - 1/2
sage: jacobi_P(2,1,2,1.2)
5.01000000000000
class sage.functions.orthogonal_polys.Func_krawtchouk#

Bases: sage.functions.orthogonal_polys.OrthogonalFunction

Krawtchouk polynomials Kj(x;n,p).

INPUT:

  • j – the degree

  • x – the independent variable x

  • n – the number of discrete points

  • p – the parameter p

See also

sage.coding.delsarte_bounds.krawtchouk() K¯ln,q(x), which are related by

(q)jK¯jn,q1(x)=Kj(x;n,1q).

EXAMPLES:

We verify the orthogonality for n=4:

sage: n = 4
sage: p = SR.var('p')
sage: matrix([[sum(binomial(n,m) * p**m * (1-p)**(n-m)
....:              * krawtchouk(i,m,n,p) * krawtchouk(j,m,n,p)
....:              for m in range(n+1)).expand().factor()
....:          for i in range(n+1)] for j in range(n+1)])
[               1                0                0                0                0]
[               0     -4*(p - 1)*p                0                0                0]
[               0                0  6*(p - 1)^2*p^2                0                0]
[               0                0                0 -4*(p - 1)^3*p^3                0]
[               0                0                0                0    (p - 1)^4*p^4]

We verify the relationship between the Krawtchouk implementations:

sage: q = SR.var('q')
sage: all(codes.bounds.krawtchouk(n, 1/q, j, x)*(-q)^j
....:     == krawtchouk(j, x, n, 1-q) for j in range(n+1))
True
eval_formula(k, x, n, p)#

Evaluate self using an explicit formula.

EXAMPLES:

sage: x,n,p = var('x,n,p')
sage: krawtchouk.eval_formula(3, x, n, p).expand().collect(x)
-1/6*n^3*p^3 + 1/2*n^2*p^3 - 1/3*n*p^3 - 1/2*(n*p - 2*p + 1)*x^2
 + 1/6*x^3 + 1/6*(3*n^2*p^2 - 9*n*p^2 + 3*n*p + 6*p^2 - 6*p + 2)*x
eval_recursive(j, x, n, p, *args, **kwds)#

Return the Krawtchouk polynomial Kj(x;n,p) using the recursive formula.

EXAMPLES:

sage: x,n,p = var('x,n,p')
sage: krawtchouk.eval_recursive(0,x,n,p)
1
sage: krawtchouk.eval_recursive(1,x,n,p)
-n*p + x
sage: krawtchouk.eval_recursive(2,x,n,p).collect(x)
1/2*n^2*p^2 + 1/2*n*(p - 1)*p - n*p^2 + 1/2*n*p - 1/2*(2*n*p - 2*p + 1)*x + 1/2*x^2
sage: bool(krawtchouk.eval_recursive(2,x,n,p) == krawtchouk(2,x,n,p))
True
sage: bool(krawtchouk.eval_recursive(3,x,n,p) == krawtchouk(3,x,n,p))
True
sage: bool(krawtchouk.eval_recursive(4,x,n,p) == krawtchouk(4,x,n,p))
True
sage: M = matrix([[-1/2,-1],[1,0]])
sage: krawtchouk.eval_recursive(2, M, 3, 1/2)
[ 9/8  7/4]
[-7/4  1/4]
class sage.functions.orthogonal_polys.Func_laguerre#

Bases: sage.functions.orthogonal_polys.OrthogonalFunction

REFERENCE:

  • [AS1964] 22.5.16, page 778 and page 789.

class sage.functions.orthogonal_polys.Func_legendre_P#

Bases: sage.symbolic.function.GinacFunction

EXAMPLES:

sage: legendre_P(4, 2.0)
55.3750000000000
sage: legendre_P(1, x)
x
sage: legendre_P(4, x+1)
35/8*(x + 1)^4 - 15/4*(x + 1)^2 + 3/8
sage: legendre_P(1/2, I+1.)
1.05338240025858 + 0.359890322109665*I
sage: legendre_P(0, SR(1)).parent()
Symbolic Ring

sage: legendre_P(0, 0)
1
sage: legendre_P(1, x)
x

sage: legendre_P(4, 2.)
55.3750000000000
sage: legendre_P(5.5,1.00001)
1.00017875754114
sage: legendre_P(1/2, I+1).n()
1.05338240025858 + 0.359890322109665*I
sage: legendre_P(1/2, I+1).n(59)
1.0533824002585801 + 0.35989032210966539*I
sage: legendre_P(42, RR(12345678))
2.66314881466753e309
sage: legendre_P(42, Reals(20)(12345678))
2.6632e309
sage: legendre_P(201/2, 0).n()
0.0561386178630179
sage: legendre_P(201/2, 0).n(100)
0.056138617863017877699963095883

sage: R.<x> = QQ[]
sage: legendre_P(4,x)
35/8*x^4 - 15/4*x^2 + 3/8
sage: legendre_P(10000,x).coefficient(x,1)
0
sage: var('t,x')
(t, x)
sage: legendre_P(-5,t)
35/8*t^4 - 15/4*t^2 + 3/8
sage: legendre_P(4, x+1)
35/8*(x + 1)^4 - 15/4*(x + 1)^2 + 3/8
sage: legendre_P(4, sqrt(2))
83/8
sage: legendre_P(4, I*e)
35/8*e^4 + 15/4*e^2 + 3/8

sage: n = var('n')
sage: derivative(legendre_P(n,x), x)
(n*x*legendre_P(n, x) - n*legendre_P(n - 1, x))/(x^2 - 1)
sage: derivative(legendre_P(3,x), x)
15/2*x^2 - 3/2
sage: derivative(legendre_P(n,x), n)
Traceback (most recent call last):
...
RuntimeError: derivative w.r.t. to the index is not supported yet
class sage.functions.orthogonal_polys.Func_legendre_Q#

Bases: sage.symbolic.function.BuiltinFunction

EXAMPLES:

sage: loads(dumps(legendre_Q))
legendre_Q
sage: maxima(legendre_Q(20,x, hold=True))._sage_().coefficient(x,10)
-29113619535/131072*log(-(x + 1)/(x - 1))
eval_formula(n, arg, **kwds)#

Return expanded Legendre Q(n, arg) function expression.

REFERENCE:

EXAMPLES:

sage: legendre_Q.eval_formula(1, x)
1/2*x*(log(x + 1) - log(-x + 1)) - 1
sage: legendre_Q.eval_formula(2,x).expand().collect(log(1+x)).collect(log(1-x))
1/4*(3*x^2 - 1)*log(x + 1) - 1/4*(3*x^2 - 1)*log(-x + 1) - 3/2*x
sage: legendre_Q.eval_formula(20,x).coefficient(x,10)
-29113619535/131072*log(x + 1) + 29113619535/131072*log(-x + 1)
sage: legendre_Q(0, 2)
-1/2*I*pi + 1/2*log(3)
sage: legendre_Q(0, 2.)
0.549306144334055 - 1.57079632679490*I
eval_recursive(n, arg, **kwds)#

Return expanded Legendre Q(n, arg) function expression.

EXAMPLES:

sage: legendre_Q.eval_recursive(2,x)
3/4*x^2*(log(x + 1) - log(-x + 1)) - 3/2*x - 1/4*log(x + 1) + 1/4*log(-x + 1)
sage: legendre_Q.eval_recursive(20,x).expand().coefficient(x,10)
-29113619535/131072*log(x + 1) + 29113619535/131072*log(-x + 1)
class sage.functions.orthogonal_polys.Func_meixner#

Bases: sage.functions.orthogonal_polys.OrthogonalFunction

Meixner polynomials Mn(x;b,c).

INPUT:

  • n – the degree

  • x – the independent variable x

  • b, c – the parameters b,c

eval_formula(n, x, b, c)#

Evaluate self using an explicit formula.

EXAMPLES:

sage: x,b,c = var('x,b,c')
sage: meixner.eval_formula(3, x, b, c).expand().collect(x)
-x^3*(3/c - 3/c^2 + 1/c^3 - 1) + b^3
 + 3*(b - 2*b/c + b/c^2 - 1/c - 1/c^2 + 1/c^3 + 1)*x^2 + 3*b^2
 + (3*b^2 + 6*b - 3*b^2/c - 3*b/c - 3*b/c^2 - 2/c^3 + 2)*x + 2*b
eval_recursive(n, x, b, c, *args, **kwds)#

Return the Meixner polynomial Mn(x;b,c) using the recursive formula.

EXAMPLES:

sage: x,b,c = var('x,b,c')
sage: meixner.eval_recursive(0,x,b,c)
1
sage: meixner.eval_recursive(1,x,b,c)
-x*(1/c - 1) + b
sage: meixner.eval_recursive(2,x,b,c).simplify_full().collect(x)
-x^2*(2/c - 1/c^2 - 1) + b^2 + (2*b - 2*b/c - 1/c^2 + 1)*x + b
sage: bool(meixner(2,x,b,c) == meixner.eval_recursive(2,x,b,c))
True
sage: bool(meixner(3,x,b,c) == meixner.eval_recursive(3,x,b,c))
True
sage: bool(meixner(4,x,b,c) == meixner.eval_recursive(4,x,b,c))
True
sage: M = matrix([[-1/2,-1],[1,0]])
sage: ret = meixner.eval_recursive(2, M, b, c).simplify_full().factor()
sage: for i in range(2):  # make the output polynomials in 1/c
....:     for j in range(2):
....:         ret[i,j] = ret[i,j].collect(c)
sage: ret
[b^2 + 1/2*(2*b + 3)/c - 1/4/c^2 - 5/4    -2*b + (2*b - 1)/c + 3/2/c^2 - 1/2]
[    2*b - (2*b - 1)/c - 3/2/c^2 + 1/2             b^2 + b + 2/c - 1/c^2 - 1]
class sage.functions.orthogonal_polys.Func_ultraspherical#

Bases: sage.symbolic.function.GinacFunction

Return the ultraspherical (or Gegenbauer) polynomial gegenbauer(n,a,x),

Cna(x)=k=0n/2(1)kΓ(nk+a)Γ(a)k!(n2k)!(2x)n2k.

When n is a nonnegative integer, this formula gives a polynomial in z of degree n, but all parameters are permitted to be complex numbers. When a=1/2, the Gegenbauer polynomial reduces to a Legendre polynomial.

Computed using Pynac.

For numerical evaluation, consider using the mpmath library, as it also allows complex numbers (and negative n as well); see the examples below.

REFERENCE:

EXAMPLES:

sage: ultraspherical(8, 101/11, x)
795972057547264/214358881*x^8 - 62604543852032/19487171*x^6...
sage: x = PolynomialRing(QQ, 'x').gen()
sage: ultraspherical(2,3/2,x)
15/2*x^2 - 3/2
sage: ultraspherical(1,1,x)
2*x
sage: t = PolynomialRing(RationalField(),"t").gen()
sage: gegenbauer(3,2,t)
32*t^3 - 12*t
sage: x = SR.var('x')
sage: n = ZZ.random_element(5, 5001)
sage: a = QQ.random_element().abs() + 5
sage: s = (  (n+1)*ultraspherical(n+1,a,x)
....:      - 2*x*(n+a)*ultraspherical(n,a,x)
....:      + (n+2*a-1)*ultraspherical(n-1,a,x) )
sage: s.expand().is_zero()
True
sage: ultraspherical(5,9/10,3.1416)
6949.55439044240
sage: ultraspherical(5,9/10,RealField(100)(pi))
6949.4695419382702451843080687

sage: a,n = SR.var('a,n')
sage: gegenbauer(2,a,x)
2*(a + 1)*a*x^2 - a
sage: gegenbauer(3,a,x)
4/3*(a + 2)*(a + 1)*a*x^3 - 2*(a + 1)*a*x
sage: gegenbauer(3,a,x).expand()
4/3*a^3*x^3 + 4*a^2*x^3 + 8/3*a*x^3 - 2*a^2*x - 2*a*x
sage: gegenbauer(10,a,x).expand().coefficient(x,2)
1/12*a^6 + 5/4*a^5 + 85/12*a^4 + 75/4*a^3 + 137/6*a^2 + 10*a
sage: ex = gegenbauer(100,a,x)
sage: (ex.subs(a==55/98) - gegenbauer(100,55/98,x)).is_trivial_zero()
True

sage: gegenbauer(2,-3,x)
12*x^2 + 3
sage: gegenbauer(120,-99/2,3)
1654502372608570682112687530178328494861923493372493824
sage: gegenbauer(5,9/2,x)
21879/8*x^5 - 6435/4*x^3 + 1287/8*x
sage: gegenbauer(15,3/2,5)
3903412392243800

sage: derivative(gegenbauer(n,a,x),x)
2*a*gegenbauer(n - 1, a + 1, x)
sage: derivative(gegenbauer(3,a,x),x)
4*(a + 2)*(a + 1)*a*x^2 - 2*(a + 1)*a
sage: derivative(gegenbauer(n,a,x),a)
Traceback (most recent call last):
...
RuntimeError: derivative w.r.t. to the second index is not supported yet

Numerical evaluation with the mpmath library:

sage: from mpmath import gegenbauer as gegenbauer_mp
sage: from mpmath import mp
sage: mp.pretty = True; mp.dps=25
sage: gegenbauer_mp(-7,0.5,0.3)
0.1291811875
sage: gegenbauer_mp(2+3j, -0.75, -1000j)
(-5038991.358609026523401901 + 9414549.285447104177860806j)
class sage.functions.orthogonal_polys.OrthogonalFunction(name, nargs=2, latex_name=None, conversions=None)#

Bases: sage.symbolic.function.BuiltinFunction

Base class for orthogonal polynomials.

This class is an abstract base class for all orthogonal polynomials since they share similar properties. The evaluation as a polynomial is either done via maxima, or with pynac.

Convention: The first argument is always the order of the polynomial, the others are other values or parameters where the polynomial is evaluated.

eval_formula(*args)#

Evaluate this polynomial using an explicit formula.

EXAMPLES:

sage: from sage.functions.orthogonal_polys import OrthogonalFunction
sage: P = OrthogonalFunction('testo_P')
sage: P.eval_formula(1,2.0)
Traceback (most recent call last):
...
NotImplementedError: no explicit calculation of values implemented