Source code for sympy.functions.special.polynomials

"""
This module mainly implements special orthogonal polynomials.

See also functions.combinatorial.numbers which contains some
combinatorial polynomials.

"""

from __future__ import print_function, division

from sympy.core.basic import C
from sympy.core.singleton import S
from sympy.core import Rational
from sympy.core.function import Function, ArgumentIndexError
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.functions.special.gamma_functions import gamma
from sympy.functions.combinatorial.factorials import factorial

from sympy.polys.orthopolys import (
    jacobi_poly,
    gegenbauer_poly,
    chebyshevt_poly,
    chebyshevu_poly,
    laguerre_poly,
    hermite_poly,
    legendre_poly
)

_x = C.Dummy('x')


class OrthogonalPolynomial(Function):
    """Base class for orthogonal polynomials.
    """

    @classmethod
    def _eval_at_order(cls, n, x):
        if n.is_integer and n >= 0:
            return cls._ortho_poly(int(n), _x).subs(_x, x)

    def _eval_conjugate(self):
        return self.func(self.args[0], self.args[1].conjugate())

#----------------------------------------------------------------------------
# Jacobi polynomials
#


[docs]class jacobi(OrthogonalPolynomial): r""" Jacobi polynomial :math:`P_n^{\left(\alpha, \beta\right)}(x)` jacobi(n, alpha, beta, x) gives the nth Jacobi polynomial in x, :math:`P_n^{\left(\alpha, \beta\right)}(x)`. The Jacobi polynomials are orthogonal on :math:`[-1, 1]` with respect to the weight :math:`\left(1-x\right)^\alpha \left(1+x\right)^\beta`. Examples ======== >>> from sympy import jacobi, S, conjugate, diff >>> from sympy.abc import n,a,b,x >>> jacobi(0, a, b, x) 1 >>> jacobi(1, a, b, x) a/2 - b/2 + x*(a/2 + b/2 + 1) >>> jacobi(2, a, b, x) # doctest:+SKIP (a**2/8 - a*b/4 - a/8 + b**2/8 - b/8 + x**2*(a**2/8 + a*b/4 + 7*a/8 + b**2/8 + 7*b/8 + 3/2) + x*(a**2/4 + 3*a/4 - b**2/4 - 3*b/4) - 1/2) >>> jacobi(n, a, b, x) jacobi(n, a, b, x) >>> jacobi(n, a, a, x) RisingFactorial(a + 1, n)*gegenbauer(n, a + 1/2, x)/RisingFactorial(2*a + 1, n) >>> jacobi(n, 0, 0, x) legendre(n, x) >>> jacobi(n, S(1)/2, S(1)/2, x) RisingFactorial(3/2, n)*chebyshevu(n, x)/factorial(n + 1) >>> jacobi(n, -S(1)/2, -S(1)/2, x) RisingFactorial(1/2, n)*chebyshevt(n, x)/factorial(n) >>> jacobi(n, a, b, -x) (-1)**n*jacobi(n, b, a, x) >>> jacobi(n, a, b, 0) 2**(-n)*gamma(a + n + 1)*hyper((-b - n, -n), (a + 1,), -1)/(factorial(n)*gamma(a + 1)) >>> jacobi(n, a, b, 1) RisingFactorial(a + 1, n)/factorial(n) >>> conjugate(jacobi(n, a, b, x)) jacobi(n, conjugate(a), conjugate(b), conjugate(x)) >>> diff(jacobi(n,a,b,x), x) (a/2 + b/2 + n/2 + 1/2)*jacobi(n - 1, a + 1, b + 1, x) See Also ======== gegenbauer, chebyshevt_root, chebyshevu, chebyshevu_root, legendre, assoc_legendre, hermite, laguerre, assoc_laguerre, sympy.polys.orthopolys.jacobi_poly, sympy.polys.orthopolys.gegenbauer_poly sympy.polys.orthopolys.chebyshevt_poly sympy.polys.orthopolys.chebyshevu_poly sympy.polys.orthopolys.hermite_poly sympy.polys.orthopolys.legendre_poly sympy.polys.orthopolys.laguerre_poly References ========== .. [1] http://en.wikipedia.org/wiki/Jacobi_polynomials .. [2] http://mathworld.wolfram.com/JacobiPolynomial.html .. [3] http://functions.wolfram.com/Polynomials/JacobiP/ """ @classmethod def eval(cls, n, a, b, x): # Simplify to other polynomials # P^{a, a}_n(x) if a == b: if a == -S.Half: return C.RisingFactorial(S.Half, n) / C.factorial(n) * chebyshevt(n, x) elif a == S.Zero: return legendre(n, x) elif a == S.Half: return C.RisingFactorial(3*S.Half, n) / C.factorial(n + 1) * chebyshevu(n, x) else: return C.RisingFactorial(a + 1, n) / C.RisingFactorial(2*a + 1, n) * gegenbauer(n, a + S.Half, x) elif b == -a: # P^{a, -a}_n(x) return C.gamma(n + a + 1) / C.gamma(n + 1) * (1 + x)**(a/2) / (1 - x)**(a/2) * assoc_legendre(n, -a, x) elif a == -b: # P^{-b, b}_n(x) return C.gamma(n - b + 1) / C.gamma(n + 1) * (1 - x)**(b/2) / (1 + x)**(b/2) * assoc_legendre(n, b, x) if not n.is_Number: # Symbolic result P^{a,b}_n(x) # P^{a,b}_n(-x) ---> (-1)**n * P^{b,a}_n(-x) if x.could_extract_minus_sign(): return S.NegativeOne**n * jacobi(n, b, a, -x) # We can evaluate for some special values of x if x == S.Zero: return (2**(-n) * C.gamma(a + n + 1) / (C.gamma(a + 1) * C.factorial(n)) * C.hyper([-b - n, -n], [a + 1], -1)) if x == S.One: return C.RisingFactorial(a + 1, n) / C.factorial(n) elif x == S.Infinity: if n.is_positive: # TODO: Make sure a+b+2*n \notin Z return C.RisingFactorial(a + b + n + 1, n) * S.Infinity else: # n is a given fixed integer, evaluate into polynomial return jacobi_poly(n, a, b, x) def fdiff(self, argindex=4): if argindex == 1: # Diff wrt n raise ArgumentIndexError(self, argindex) elif argindex == 2: # Diff wrt a n, a, b, x = self.args k = C.Dummy("k") f1 = 1 / (a + b + n + k + 1) f2 = ((a + b + 2*k + 1) * C.RisingFactorial(b + k + 1, n - k) / ((n - k) * C.RisingFactorial(a + b + k + 1, n - k))) return C.Sum(f1 * (jacobi(n, a, b, x) + f2*jacobi(k, a, b, x)), (k, 0, n - 1)) elif argindex == 3: # Diff wrt b n, a, b, x = self.args k = C.Dummy("k") f1 = 1 / (a + b + n + k + 1) f2 = (-1)**(n - k) * ((a + b + 2*k + 1) * C.RisingFactorial(a + k + 1, n - k) / ((n - k) * C.RisingFactorial(a + b + k + 1, n - k))) return C.Sum(f1 * (jacobi(n, a, b, x) + f2*jacobi(k, a, b, x)), (k, 0, n - 1)) elif argindex == 4: # Diff wrt x n, a, b, x = self.args return S.Half * (a + b + n + 1) * jacobi(n - 1, a + 1, b + 1, x) else: raise ArgumentIndexError(self, argindex) def _eval_rewrite_as_polynomial(self, n, a, b, x): # TODO: Make sure n \in N k = C.Dummy("k") kern = (C.RisingFactorial(-n, k) * C.RisingFactorial(a + b + n + 1, k) * C.RisingFactorial(a + k + 1, n - k) / C.factorial(k) * ((1 - x)/2)**k) return 1 / C.factorial(n) * C.Sum(kern, (k, 0, n)) def _eval_conjugate(self): n, a, b, x = self.args return self.func(n, a.conjugate(), b.conjugate(), x.conjugate())
[docs]def jacobi_normalized(n, a, b, x): r""" Jacobi polynomial :math:`P_n^{\left(\alpha, \beta\right)}(x)` jacobi_normalized(n, alpha, beta, x) gives the nth Jacobi polynomial in x, :math:`P_n^{\left(\alpha, \beta\right)}(x)`. The Jacobi polynomials are orthogonal on :math:`[-1, 1]` with respect to the weight :math:`\left(1-x\right)^\alpha \left(1+x\right)^\beta`. This functions returns the polynomials normilzed: .. math:: \int_{-1}^{1} P_m^{\left(\alpha, \beta\right)}(x) P_n^{\left(\alpha, \beta\right)}(x) (1-x)^{\alpha} (1+x)^{\beta} \mathrm{d}x = \delta_{m,n} Examples ======== >>> from sympy import jacobi_normalized >>> from sympy.abc import n,a,b,x >>> jacobi_normalized(n, a, b, x) jacobi(n, a, b, x)/sqrt(2**(a + b + 1)*gamma(a + n + 1)*gamma(b + n + 1)/((a + b + 2*n + 1)*factorial(n)*gamma(a + b + n + 1))) See Also ======== gegenbauer, chebyshevt_root, chebyshevu, chebyshevu_root, legendre, assoc_legendre, hermite, laguerre, assoc_laguerre, sympy.polys.orthopolys.jacobi_poly, sympy.polys.orthopolys.gegenbauer_poly sympy.polys.orthopolys.chebyshevt_poly sympy.polys.orthopolys.chebyshevu_poly sympy.polys.orthopolys.hermite_poly sympy.polys.orthopolys.legendre_poly sympy.polys.orthopolys.laguerre_poly References ========== .. [1] http://en.wikipedia.org/wiki/Jacobi_polynomials .. [2] http://mathworld.wolfram.com/JacobiPolynomial.html .. [3] http://functions.wolfram.com/Polynomials/JacobiP/ """ nfactor = (S(2)**(a + b + 1) * (gamma(n + a + 1) * gamma(n + b + 1)) / (2*n + a + b + 1) / (factorial(n) * gamma(n + a + b + 1))) return jacobi(n, a, b, x) / sqrt(nfactor) #---------------------------------------------------------------------------- # Gegenbauer polynomials #
[docs]class gegenbauer(OrthogonalPolynomial): r""" Gegenbauer polynomial :math:`C_n^{\left(\alpha\right)}(x)` gegenbauer(n, alpha, x) gives the nth Gegenbauer polynomial in x, :math:`C_n^{\left(\alpha\right)}(x)`. The Gegenbauer polynomials are orthogonal on :math:`[-1, 1]` with respect to the weight :math:`\left(1-x^2\right)^{\alpha-\frac{1}{2}}`. Examples ======== >>> from sympy import gegenbauer, conjugate, diff >>> from sympy.abc import n,a,x >>> gegenbauer(0, a, x) 1 >>> gegenbauer(1, a, x) 2*a*x >>> gegenbauer(2, a, x) -a + x**2*(2*a**2 + 2*a) >>> gegenbauer(3, a, x) x**3*(4*a**3/3 + 4*a**2 + 8*a/3) + x*(-2*a**2 - 2*a) >>> gegenbauer(n, a, x) gegenbauer(n, a, x) >>> gegenbauer(n, a, -x) (-1)**n*gegenbauer(n, a, x) >>> gegenbauer(n, a, 0) 2**n*sqrt(pi)*gamma(a + n/2)/(gamma(a)*gamma(-n/2 + 1/2)*gamma(n + 1)) >>> gegenbauer(n, a, 1) gamma(2*a + n)/(gamma(2*a)*gamma(n + 1)) >>> conjugate(gegenbauer(n, a, x)) gegenbauer(n, conjugate(a), conjugate(x)) >>> diff(gegenbauer(n, a, x), x) 2*a*gegenbauer(n - 1, a + 1, x) See Also ======== jacobi, chebyshevt_root, chebyshevu, chebyshevu_root, legendre, assoc_legendre, hermite, laguerre, assoc_laguerre, sympy.polys.orthopolys.jacobi_poly sympy.polys.orthopolys.gegenbauer_poly sympy.polys.orthopolys.chebyshevt_poly sympy.polys.orthopolys.chebyshevu_poly sympy.polys.orthopolys.hermite_poly sympy.polys.orthopolys.legendre_poly sympy.polys.orthopolys.laguerre_poly References ========== .. [1] http://en.wikipedia.org/wiki/Gegenbauer_polynomials .. [2] http://mathworld.wolfram.com/GegenbauerPolynomial.html .. [3] http://functions.wolfram.com/Polynomials/GegenbauerC3/ """ @classmethod def eval(cls, n, a, x): # For negative n the polynomials vanish # See http://functions.wolfram.com/Polynomials/GegenbauerC3/03/01/03/0012/ if n.is_negative: return S.Zero # Some special values for fixed a if a == S.Half: return legendre(n, x) elif a == S.One: return chebyshevu(n, x) elif a == S.NegativeOne: return S.Zero if not n.is_Number: # Handle this before the general sign extraction rule if x == S.NegativeOne: if (C.re(a) > S.Half) == True: return S.ComplexInfinity else: # No sec function available yet #return (C.cos(S.Pi*(a+n)) * C.sec(S.Pi*a) * C.gamma(2*a+n) / # (C.gamma(2*a) * C.gamma(n+1))) return None # Symbolic result C^a_n(x) # C^a_n(-x) ---> (-1)**n * C^a_n(x) if x.could_extract_minus_sign(): return S.NegativeOne**n * gegenbauer(n, a, -x) # We can evaluate for some special values of x if x == S.Zero: return (2**n * sqrt(S.Pi) * C.gamma(a + S.Half*n) / (C.gamma((1 - n)/2) * C.gamma(n + 1) * C.gamma(a)) ) if x == S.One: return C.gamma(2*a + n) / (C.gamma(2*a) * C.gamma(n + 1)) elif x == S.Infinity: if n.is_positive: return C.RisingFactorial(a, n) * S.Infinity else: # n is a given fixed integer, evaluate into polynomial return gegenbauer_poly(n, a, x) def fdiff(self, argindex=3): if argindex == 1: # Diff wrt n raise ArgumentIndexError(self, argindex) elif argindex == 2: # Diff wrt a n, a, x = self.args k = C.Dummy("k") factor1 = 2 * (1 + (-1)**(n - k)) * (k + a) / ((k + n + 2*a) * (n - k)) factor2 = 2*(k + 1) / ((k + 2*a) * (2*k + 2*a + 1)) + \ 2 / (k + n + 2*a) kern = factor1*gegenbauer(k, a, x) + factor2*gegenbauer(n, a, x) return C.Sum(kern, (k, 0, n - 1)) elif argindex == 3: # Diff wrt x n, a, x = self.args return 2*a*gegenbauer(n - 1, a + 1, x) else: raise ArgumentIndexError(self, argindex) def _eval_rewrite_as_polynomial(self, n, a, x): k = C.Dummy("k") kern = ((-1)**k * C.RisingFactorial(a, n - k) * (2*x)**(n - 2*k) / (C.factorial(k) * C.factorial(n - 2*k))) return C.Sum(kern, (k, 0, C.floor(n/2))) def _eval_conjugate(self): n, a, x = self.args return self.func(n, a.conjugate(), x.conjugate()) #---------------------------------------------------------------------------- # Chebyshev polynomials of first and second kind #
[docs]class chebyshevt(OrthogonalPolynomial): r""" Chebyshev polynomial of the first kind, :math:`T_n(x)` chebyshevt(n, x) gives the nth Chebyshev polynomial (of the first kind) in x, :math:`T_n(x)`. The Chebyshev polynomials of the first kind are orthogonal on :math:`[-1, 1]` with respect to the weight :math:`\frac{1}{\sqrt{1-x^2}}`. Examples ======== >>> from sympy import chebyshevt, chebyshevu, diff >>> from sympy.abc import n,x >>> chebyshevt(0, x) 1 >>> chebyshevt(1, x) x >>> chebyshevt(2, x) 2*x**2 - 1 >>> chebyshevt(n, x) chebyshevt(n, x) >>> chebyshevt(n, -x) (-1)**n*chebyshevt(n, x) >>> chebyshevt(-n, x) chebyshevt(n, x) >>> chebyshevt(n, 0) cos(pi*n/2) >>> chebyshevt(n, -1) (-1)**n >>> diff(chebyshevt(n, x), x) n*chebyshevu(n - 1, x) See Also ======== jacobi, gegenbauer, chebyshevt_root, chebyshevu, chebyshevu_root, legendre, assoc_legendre, hermite, laguerre, assoc_laguerre, sympy.polys.orthopolys.jacobi_poly sympy.polys.orthopolys.gegenbauer_poly sympy.polys.orthopolys.chebyshevt_poly sympy.polys.orthopolys.chebyshevu_poly sympy.polys.orthopolys.hermite_poly sympy.polys.orthopolys.legendre_poly sympy.polys.orthopolys.laguerre_poly References ========== .. [1] http://en.wikipedia.org/wiki/Chebyshev_polynomial .. [2] http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html .. [3] http://mathworld.wolfram.com/ChebyshevPolynomialoftheSecondKind.html .. [4] http://functions.wolfram.com/Polynomials/ChebyshevT/ .. [5] http://functions.wolfram.com/Polynomials/ChebyshevU/ """ _ortho_poly = staticmethod(chebyshevt_poly) @classmethod def eval(cls, n, x): if not n.is_Number: # Symbolic result T_n(x) # T_n(-x) ---> (-1)**n * T_n(x) if x.could_extract_minus_sign(): return S.NegativeOne**n * chebyshevt(n, -x) # T_{-n}(x) ---> T_n(x) if n.could_extract_minus_sign(): return chebyshevt(-n, x) # We can evaluate for some special values of x if x == S.Zero: return C.cos(S.Half * S.Pi * n) if x == S.One: return S.One elif x == S.Infinity: return S.Infinity else: # n is a given fixed integer, evaluate into polynomial if n.is_negative: # T_{-n}(x) == T_n(x) return cls._eval_at_order(-n, x) else: return cls._eval_at_order(n, x) def fdiff(self, argindex=2): if argindex == 1: # Diff wrt n raise ArgumentIndexError(self, argindex) elif argindex == 2: # Diff wrt x n, x = self.args return n * chebyshevu(n - 1, x) else: raise ArgumentIndexError(self, argindex) def _eval_rewrite_as_polynomial(self, n, x): k = C.Dummy("k") kern = C.binomial(n, 2*k) * (x**2 - 1)**k * x**(n - 2*k) return C.Sum(kern, (k, 0, C.floor(n/2)))
[docs]class chebyshevu(OrthogonalPolynomial): r""" Chebyshev polynomial of the second kind, :math:`U_n(x)` chebyshevu(n, x) gives the nth Chebyshev polynomial of the second kind in x, :math:`U_n(x)`. The Chebyshev polynomials of the second kind are orthogonal on :math:`[-1, 1]` with respect to the weight :math:`\sqrt{1-x^2}`. Examples ======== >>> from sympy import chebyshevt, chebyshevu, diff >>> from sympy.abc import n,x >>> chebyshevu(0, x) 1 >>> chebyshevu(1, x) 2*x >>> chebyshevu(2, x) 4*x**2 - 1 >>> chebyshevu(n, x) chebyshevu(n, x) >>> chebyshevu(n, -x) (-1)**n*chebyshevu(n, x) >>> chebyshevu(-n, x) -chebyshevu(n - 2, x) >>> chebyshevu(n, 0) cos(pi*n/2) >>> chebyshevu(n, 1) n + 1 >>> diff(chebyshevu(n, x), x) (-x*chebyshevu(n, x) + (n + 1)*chebyshevt(n + 1, x))/(x**2 - 1) See Also ======== jacobi, gegenbauer, chebyshevt, chebyshevt_root, chebyshevu_root, legendre, assoc_legendre, hermite, laguerre, assoc_laguerre, sympy.polys.orthopolys.jacobi_poly sympy.polys.orthopolys.gegenbauer_poly sympy.polys.orthopolys.chebyshevt_poly sympy.polys.orthopolys.chebyshevu_poly sympy.polys.orthopolys.hermite_poly sympy.polys.orthopolys.legendre_poly sympy.polys.orthopolys.laguerre_poly References ========== .. [1] http://en.wikipedia.org/wiki/Chebyshev_polynomial .. [2] http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html .. [3] http://mathworld.wolfram.com/ChebyshevPolynomialoftheSecondKind.html .. [4] http://functions.wolfram.com/Polynomials/ChebyshevT/ .. [5] http://functions.wolfram.com/Polynomials/ChebyshevU/ """ _ortho_poly = staticmethod(chebyshevu_poly) @classmethod def eval(cls, n, x): if not n.is_Number: # Symbolic result U_n(x) # U_n(-x) ---> (-1)**n * U_n(x) if x.could_extract_minus_sign(): return S.NegativeOne**n * chebyshevu(n, -x) # U_{-n}(x) ---> -U_{n-2}(x) if n.could_extract_minus_sign(): if n == S.NegativeOne: return S.Zero else: return -chebyshevu(-n - 2, x) # We can evaluate for some special values of x if x == S.Zero: return C.cos(S.Half * S.Pi * n) if x == S.One: return S.One + n elif x == S.Infinity: return S.Infinity else: # n is a given fixed integer, evaluate into polynomial if n.is_negative: # U_{-n}(x) ---> -U_{n-2}(x) if n == S.NegativeOne: return S.Zero else: return -cls._eval_at_order(-n - 2, x) else: return cls._eval_at_order(n, x) def fdiff(self, argindex=2): if argindex == 1: # Diff wrt n raise ArgumentIndexError(self, argindex) elif argindex == 2: # Diff wrt x n, x = self.args return ((n + 1) * chebyshevt(n + 1, x) - x * chebyshevu(n, x)) / (x**2 - 1) else: raise ArgumentIndexError(self, argindex) def _eval_rewrite_as_polynomial(self, n, x): k = C.Dummy("k") kern = S.NegativeOne**k * C.factorial( n - k) * (2*x)**(n - 2*k) / (C.factorial(k) * C.factorial(n - 2*k)) return C.Sum(kern, (k, 0, C.floor(n/2)))
[docs]class chebyshevt_root(Function): r""" chebyshev_root(n, k) returns the kth root (indexed from zero) of the nth Chebyshev polynomial of the first kind; that is, if 0 <= k < n, chebyshevt(n, chebyshevt_root(n, k)) == 0. Examples ======== >>> from sympy import chebyshevt, chebyshevt_root >>> chebyshevt_root(3, 2) -sqrt(3)/2 >>> chebyshevt(3, chebyshevt_root(3, 2)) 0 See Also ======== jacobi, gegenbauer, chebyshevt, chebyshevu, chebyshevu_root, legendre, assoc_legendre, hermite, laguerre, assoc_laguerre, sympy.polys.orthopolys.jacobi_poly sympy.polys.orthopolys.gegenbauer_poly sympy.polys.orthopolys.chebyshevt_poly sympy.polys.orthopolys.chebyshevu_poly sympy.polys.orthopolys.hermite_poly sympy.polys.orthopolys.legendre_poly sympy.polys.orthopolys.laguerre_poly """ @classmethod def eval(cls, n, k): if not ((0 <= k) and (k < n)): raise ValueError("must have 0 <= k < n, " "got k = %s and n = %s" % (k, n)) return C.cos(S.Pi*(2*k + 1)/(2*n))
[docs]class chebyshevu_root(Function): r""" chebyshevu_root(n, k) returns the kth root (indexed from zero) of the nth Chebyshev polynomial of the second kind; that is, if 0 <= k < n, chebyshevu(n, chebyshevu_root(n, k)) == 0. Examples ======== >>> from sympy import chebyshevu, chebyshevu_root >>> chebyshevu_root(3, 2) -sqrt(2)/2 >>> chebyshevu(3, chebyshevu_root(3, 2)) 0 See Also ======== chebyshevt, chebyshevt_root, chebyshevu, legendre, assoc_legendre, hermite, laguerre, assoc_laguerre, sympy.polys.orthopolys.jacobi_poly sympy.polys.orthopolys.gegenbauer_poly sympy.polys.orthopolys.chebyshevt_poly sympy.polys.orthopolys.chebyshevu_poly sympy.polys.orthopolys.hermite_poly sympy.polys.orthopolys.legendre_poly sympy.polys.orthopolys.laguerre_poly """ @classmethod def eval(cls, n, k): if not ((0 <= k) and (k < n)): raise ValueError("must have 0 <= k < n, " "got k = %s and n = %s" % (k, n)) return C.cos(S.Pi*(k + 1)/(n + 1)) #---------------------------------------------------------------------------- # Legendre polynomials and Associated Legendre polynomials #
[docs]class legendre(OrthogonalPolynomial): r""" legendre(n, x) gives the nth Legendre polynomial of x, :math:`P_n(x)` The Legendre polynomials are orthogonal on [-1, 1] with respect to the constant weight 1. They satisfy :math:`P_n(1) = 1` for all n; further, :math:`P_n` is odd for odd n and even for even n. Examples ======== >>> from sympy import legendre, diff >>> from sympy.abc import x, n >>> legendre(0, x) 1 >>> legendre(1, x) x >>> legendre(2, x) 3*x**2/2 - 1/2 >>> legendre(n, x) legendre(n, x) >>> diff(legendre(n,x), x) n*(x*legendre(n, x) - legendre(n - 1, x))/(x**2 - 1) See Also ======== jacobi, gegenbauer, chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root, assoc_legendre, hermite, laguerre, assoc_laguerre, sympy.polys.orthopolys.jacobi_poly sympy.polys.orthopolys.gegenbauer_poly sympy.polys.orthopolys.chebyshevt_poly sympy.polys.orthopolys.chebyshevu_poly sympy.polys.orthopolys.hermite_poly sympy.polys.orthopolys.legendre_poly sympy.polys.orthopolys.laguerre_poly References ========== .. [1] http://en.wikipedia.org/wiki/Legendre_polynomial .. [2] http://mathworld.wolfram.com/LegendrePolynomial.html .. [3] http://functions.wolfram.com/Polynomials/LegendreP/ .. [4] http://functions.wolfram.com/Polynomials/LegendreP2/ """ _ortho_poly = staticmethod(legendre_poly) @classmethod def eval(cls, n, x): if not n.is_Number: # Symbolic result L_n(x) # L_n(-x) ---> (-1)**n * L_n(x) if x.could_extract_minus_sign(): return S.NegativeOne**n * legendre(n, -x) # L_{-n}(x) ---> L_{n-1}(x) if n.could_extract_minus_sign(): return legendre(-n - S.One, x) # We can evaluate for some special values of x if x == S.Zero: return sqrt(S.Pi)/(C.gamma(S.Half - n/2)*C.gamma(S.One + n/2)) elif x == S.One: return S.One elif x == S.Infinity: return S.Infinity else: # n is a given fixed integer, evaluate into polynomial if n.is_negative: raise ValueError( "The index n must be nonnegative integer (got %r)" % n) else: return cls._eval_at_order(n, x) def fdiff(self, argindex=2): if argindex == 1: # Diff wrt n raise ArgumentIndexError(self, argindex) elif argindex == 2: # Diff wrt x # Find better formula, this is unsuitable for x = 1 n, x = self.args return n/(x**2 - 1)*(x*legendre(n, x) - legendre(n - 1, x)) else: raise ArgumentIndexError(self, argindex) def _eval_rewrite_as_polynomial(self, n, x): k = C.Dummy("k") kern = (-1)**k*C.binomial(n, k)**2*((1 + x)/2)**(n - k)*((1 - x)/2)**k return C.Sum(kern, (k, 0, n))
[docs]class assoc_legendre(Function): r""" assoc_legendre(n,m, x) gives :math:`P_n^m(x)`, where n and m are the degree and order or an expression which is related to the nth order Legendre polynomial, :math:`P_n(x)` in the following manner: .. math:: P_n^m(x) = (-1)^m (1 - x^2)^{\frac{m}{2}} \frac{\mathrm{d}^m P_n(x)}{\mathrm{d} x^m} Associated Legendre polynomial are orthogonal on [-1, 1] with: - weight = 1 for the same m, and different n. - weight = 1/(1-x**2) for the same n, and different m. Examples ======== >>> from sympy import assoc_legendre >>> from sympy.abc import x, m, n >>> assoc_legendre(0,0, x) 1 >>> assoc_legendre(1,0, x) x >>> assoc_legendre(1,1, x) -sqrt(-x**2 + 1) >>> assoc_legendre(n,m,x) assoc_legendre(n, m, x) See Also ======== jacobi, gegenbauer, chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root, legendre, hermite, laguerre, assoc_laguerre, sympy.polys.orthopolys.jacobi_poly sympy.polys.orthopolys.gegenbauer_poly sympy.polys.orthopolys.chebyshevt_poly sympy.polys.orthopolys.chebyshevu_poly sympy.polys.orthopolys.hermite_poly sympy.polys.orthopolys.legendre_poly sympy.polys.orthopolys.laguerre_poly References ========== .. [1] http://en.wikipedia.org/wiki/Associated_Legendre_polynomials .. [2] http://mathworld.wolfram.com/LegendrePolynomial.html .. [3] http://functions.wolfram.com/Polynomials/LegendreP/ .. [4] http://functions.wolfram.com/Polynomials/LegendreP2/ """ @classmethod def _eval_at_order(cls, n, m): P = legendre_poly(n, _x, polys=True).diff((_x, m)) return (-1)**m * (1 - _x**2)**Rational(m, 2) * P.as_expr() @classmethod def eval(cls, n, m, x): if m.could_extract_minus_sign(): # P^{-m}_n ---> F * P^m_n return S.NegativeOne**(-m) * (C.factorial(m + n)/C.factorial(n - m)) * assoc_legendre(n, -m, x) if m == 0: # P^0_n ---> L_n return legendre(n, x) if x == 0: return 2**m*sqrt(S.Pi) / (C.gamma((1 - m - n)/2)*C.gamma(1 - (m - n)/2)) if n.is_Number and m.is_Number and n.is_integer and m.is_integer: if n.is_negative: raise ValueError("%s : 1st index must be nonnegative integer (got %r)" % (cls, n)) if abs(m) > n: raise ValueError("%s : abs('2nd index') must be <= '1st index' (got %r, %r)" % (cls, n, m)) return cls._eval_at_order(int(n), abs(int(m))).subs(_x, x) def fdiff(self, argindex=3): if argindex == 1: # Diff wrt n raise ArgumentIndexError(self, argindex) elif argindex == 2: # Diff wrt m raise ArgumentIndexError(self, argindex) elif argindex == 3: # Diff wrt x # Find better formula, this is unsuitable for x = 1 n, m, x = self.args return 1/(x**2 - 1)*(x*n*assoc_legendre(n, m, x) - (m + n)*assoc_legendre(n - 1, m, x)) else: raise ArgumentIndexError(self, argindex) def _eval_rewrite_as_polynomial(self, n, m, x): k = C.Dummy("k") kern = C.factorial(2*n - 2*k)/(2**n*C.factorial(n - k)*C.factorial( k)*C.factorial(n - 2*k - m))*(-1)**k*x**(n - m - 2*k) return (1 - x**2)**(m/2) * C.Sum(kern, (k, 0, C.floor((n - m)*S.Half))) def _eval_conjugate(self): n, m, x = self.args return self.func(n, m.conjugate(), x.conjugate()) #---------------------------------------------------------------------------- # Hermite polynomials #
[docs]class hermite(OrthogonalPolynomial): r""" hermite(n, x) gives the nth Hermite polynomial in x, :math:`H_n(x)` The Hermite polynomials are orthogonal on :math:`(-\infty, \infty)` with respect to the weight :math:`\exp\left(-\frac{x^2}{2}\right)`. Examples ======== >>> from sympy import hermite, diff >>> from sympy.abc import x, n >>> hermite(0, x) 1 >>> hermite(1, x) 2*x >>> hermite(2, x) 4*x**2 - 2 >>> hermite(n, x) hermite(n, x) >>> diff(hermite(n,x), x) 2*n*hermite(n - 1, x) >>> diff(hermite(n,x), x) 2*n*hermite(n - 1, x) >>> hermite(n, -x) (-1)**n*hermite(n, x) See Also ======== jacobi, gegenbauer, chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root, legendre, assoc_legendre, laguerre, assoc_laguerre, sympy.polys.orthopolys.jacobi_poly sympy.polys.orthopolys.gegenbauer_poly sympy.polys.orthopolys.chebyshevt_poly sympy.polys.orthopolys.chebyshevu_poly sympy.polys.orthopolys.hermite_poly sympy.polys.orthopolys.legendre_poly sympy.polys.orthopolys.laguerre_poly References ========== .. [1] http://en.wikipedia.org/wiki/Hermite_polynomial .. [2] http://mathworld.wolfram.com/HermitePolynomial.html .. [3] http://functions.wolfram.com/Polynomials/HermiteH/ """ _ortho_poly = staticmethod(hermite_poly) @classmethod def eval(cls, n, x): if not n.is_Number: # Symbolic result H_n(x) # H_n(-x) ---> (-1)**n * H_n(x) if x.could_extract_minus_sign(): return S.NegativeOne**n * hermite(n, -x) # We can evaluate for some special values of x if x == S.Zero: return 2**n * sqrt(S.Pi) / C.gamma((S.One - n)/2) elif x == S.Infinity: return S.Infinity else: # n is a given fixed integer, evaluate into polynomial if n.is_negative: raise ValueError( "The index n must be nonnegative integer (got %r)" % n) else: return cls._eval_at_order(n, x) def fdiff(self, argindex=2): if argindex == 1: # Diff wrt n raise ArgumentIndexError(self, argindex) elif argindex == 2: # Diff wrt x n, x = self.args return 2*n*hermite(n - 1, x) else: raise ArgumentIndexError(self, argindex) def _eval_rewrite_as_polynomial(self, n, x): k = C.Dummy("k") kern = (-1)**k / (C.factorial(k)*C.factorial(n - 2*k)) * (2*x)**(n - 2*k) return C.factorial(n)*C.Sum(kern, (k, 0, C.floor(n/2))) #---------------------------------------------------------------------------- # Laguerre polynomials #
[docs]class laguerre(OrthogonalPolynomial): r""" Returns the nth Laguerre polynomial in x, :math:`L_n(x)`. Parameters ========== n : int Degree of Laguerre polynomial. Must be ``n >= 0``. Examples ======== >>> from sympy import laguerre, diff >>> from sympy.abc import x, n >>> laguerre(0, x) 1 >>> laguerre(1, x) -x + 1 >>> laguerre(2, x) x**2/2 - 2*x + 1 >>> laguerre(3, x) -x**3/6 + 3*x**2/2 - 3*x + 1 >>> laguerre(n, x) laguerre(n, x) >>> diff(laguerre(n, x), x) -assoc_laguerre(n - 1, 1, x) See Also ======== jacobi, gegenbauer, chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root, legendre, assoc_legendre, hermite, assoc_laguerre, sympy.polys.orthopolys.jacobi_poly sympy.polys.orthopolys.gegenbauer_poly sympy.polys.orthopolys.chebyshevt_poly sympy.polys.orthopolys.chebyshevu_poly sympy.polys.orthopolys.hermite_poly sympy.polys.orthopolys.legendre_poly sympy.polys.orthopolys.laguerre_poly References ========== .. [1] http://en.wikipedia.org/wiki/Laguerre_polynomial .. [2] http://mathworld.wolfram.com/LaguerrePolynomial.html .. [3] http://functions.wolfram.com/Polynomials/LaguerreL/ .. [4] http://functions.wolfram.com/Polynomials/LaguerreL3/ """ @classmethod def eval(cls, n, x): if not n.is_Number: # Symbolic result L_n(x) # L_{n}(-x) ---> exp(-x) * L_{-n-1}(x) # L_{-n}(x) ---> exp(x) * L_{n-1}(-x) if n.could_extract_minus_sign(): return C.exp(x) * laguerre(n - 1, -x) # We can evaluate for some special values of x if x == S.Zero: return S.One elif x == S.NegativeInfinity: return S.Infinity elif x == S.Infinity: return S.NegativeOne**n * S.Infinity else: # n is a given fixed integer, evaluate into polynomial if n.is_negative: raise ValueError( "The index n must be nonnegative integer (got %r)" % n) else: return laguerre_poly(n, x, 0) def fdiff(self, argindex=2): if argindex == 1: # Diff wrt n raise ArgumentIndexError(self, argindex) elif argindex == 2: # Diff wrt x n, x = self.args return -assoc_laguerre(n - 1, 1, x) else: raise ArgumentIndexError(self, argindex) def _eval_rewrite_as_polynomial(self, n, x): # TODO: Should make sure n is in N_0 k = C.Dummy("k") kern = C.RisingFactorial(-n, k) / C.factorial(k)**2 * x**k return C.Sum(kern, (k, 0, n))
[docs]class assoc_laguerre(OrthogonalPolynomial): r""" Returns the nth generalized Laguerre polynomial in x, :math:`L_n(x)`. Parameters ========== n : int Degree of Laguerre polynomial. Must be ``n >= 0``. alpha : Expr Arbitrary expression. For ``alpha=0`` regular Laguerre polynomials will be generated. Examples ======== >>> from sympy import laguerre, assoc_laguerre, diff >>> from sympy.abc import x, n, a >>> assoc_laguerre(0, a, x) 1 >>> assoc_laguerre(1, a, x) a - x + 1 >>> assoc_laguerre(2, a, x) a**2/2 + 3*a/2 + x**2/2 + x*(-a - 2) + 1 >>> assoc_laguerre(3, a, x) a**3/6 + a**2 + 11*a/6 - x**3/6 + x**2*(a/2 + 3/2) + x*(-a**2/2 - 5*a/2 - 3) + 1 >>> assoc_laguerre(n, a, 0) binomial(a + n, a) >>> assoc_laguerre(n, a, x) assoc_laguerre(n, a, x) >>> assoc_laguerre(n, 0, x) laguerre(n, x) >>> diff(assoc_laguerre(n, a, x), x) -assoc_laguerre(n - 1, a + 1, x) >>> diff(assoc_laguerre(n, a, x), a) Sum(assoc_laguerre(_k, a, x)/(-a + n), (_k, 0, n - 1)) See Also ======== jacobi, gegenbauer, chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root, legendre, assoc_legendre, hermite, laguerre, sympy.polys.orthopolys.jacobi_poly sympy.polys.orthopolys.gegenbauer_poly sympy.polys.orthopolys.chebyshevt_poly sympy.polys.orthopolys.chebyshevu_poly sympy.polys.orthopolys.hermite_poly sympy.polys.orthopolys.legendre_poly sympy.polys.orthopolys.laguerre_poly References ========== .. [1] http://en.wikipedia.org/wiki/Laguerre_polynomial#Assoc_laguerre_polynomials .. [2] http://mathworld.wolfram.com/AssociatedLaguerrePolynomial.html .. [3] http://functions.wolfram.com/Polynomials/LaguerreL/ .. [4] http://functions.wolfram.com/Polynomials/LaguerreL3/ """ @classmethod def eval(cls, n, alpha, x): # L_{n}^{0}(x) ---> L_{n}(x) if alpha == S.Zero: return laguerre(n, x) if not n.is_Number: # We can evaluate for some special values of x if x == S.Zero: return C.binomial(n + alpha, alpha) elif x == S.Infinity and n > S.Zero: return S.NegativeOne**n * S.Infinity elif x == S.NegativeInfinity and n > S.Zero: return S.Infinity else: # n is a given fixed integer, evaluate into polynomial if n.is_negative: raise ValueError( "The index n must be nonnegative integer (got %r)" % n) else: return laguerre_poly(n, x, alpha) def fdiff(self, argindex=3): if argindex == 1: # Diff wrt n raise ArgumentIndexError(self, argindex) elif argindex == 2: # Diff wrt alpha n, alpha, x = self.args k = C.Dummy("k") return C.Sum(assoc_laguerre(k, alpha, x) / (n - alpha), (k, 0, n - 1)) elif argindex == 3: # Diff wrt x n, alpha, x = self.args return -assoc_laguerre(n - 1, alpha + 1, x) else: raise ArgumentIndexError(self, argindex) def _eval_rewrite_as_polynomial(self, n, x): # TODO: Should make sure n is in N_0 k = C.Dummy("k") kern = C.RisingFactorial( -n, k) / (C.gamma(k + alpha + 1) * C.factorial(k)) * x**k return C.gamma(n + alpha + 1) / C.factorial(n) * C.Sum(kern, (k, 0, n)) def _eval_conjugate(self): n, alpha, x = self.args return self.func(n, alpha.conjugate(), x.conjugate())