from sage.all import * from sage.symbolic.expression import is_Expression def legendre_p(n, x, a = -1, b = 1): """ Returns the ``n``th Legendre polynomial of the first kind over the interval [a, b] with respect to ``x``. When [a,b] is not [-1,1], we scale the standard Legendre polynomial (which is defined over [-1,1]) via an affine map. The resulting polynomials are still orthogonal and possess the property that `P(a) = P(b) = 1`. INPUT: * ``n`` -- The index of the polynomial. * ``x`` -- Either the variable to use as the independent variable in the polynomial, or a point at which to evaluate the polynomial. * ``a`` -- The "left" endpoint of the interval. Must be a real number. * ``b`` -- The "right" endpoint of the interval. Must be a real number. OUTPUT: If ``x`` is a variable, a polynomial (symbolic expression) will be returned. Otherwise, the value of the ``n``th polynomial at ``x`` will be returned. TESTS: We should agree with Maxima for all `n`:: sage: eq = lambda k: bool(legendre_p(k,x) == legendre_P(k,x)) sage: all([eq(k) for k in range(0,20) ]) # long time True We can evaluate the result of the zeroth polynomial:: sage: f = legendre_p(0,x) sage: f(x=10) 1 We should have |P(a)| = |P(b)| = 1 for all a,b:: sage: a = RR.random_element() sage: b = RR.random_element() sage: k = ZZ.random_element(20) sage: P = legendre_p(k, x, a, b) sage: abs(P(x=a)) # abs tol 1e-12 1 sage: abs(P(x=b)) # abs tol 1e-12 1 Two different polynomials should be orthogonal with respect to the inner product over [a,b]:: sage: a = RR.random_element() sage: b = RR.random_element() sage: j = ZZ.random_element(20) sage: k = j + 1 sage: Pj = legendre_p(j, x, a, b) sage: Pk = legendre_p(k, x, a, b) sage: integrate(Pj*Pk, x, a, b) # abs tol 1e-12 0 The first few polynomials shifted to [0,1] are known to be:: sage: p0 = 1 sage: p1 = 2*x - 1 sage: p2 = 6*x^2 - 6*x + 1 sage: p3 = 20*x^3 - 30*x^2 + 12*x - 1 sage: bool(legendre_p(0, x, 0, 1) == p0) True sage: bool(legendre_p(1, x, 0, 1) == p1) True sage: bool(legendre_p(2, x, 0, 1) == p2) True sage: bool(legendre_p(3, x, 0, 1) == p3) True """ if not n in ZZ: raise TypeError('n must be a natural number') if n < 0: raise ValueError('n must be nonnegative') if not (a in RR and b in RR): raise TypeError('both `a` and `b` must be real numbers') a = RR(a) b = RR(b) n = ZZ(n) # Ensure that 1/(2**n) is not integer division. dn = 1/(2**n) def phi(t): # This is an affine map from [a,b] into [-1,1] and so preserves # orthogonality. return (2 / (b-a))*t + 1 - (2*b)/(b-a) def c(m): return binomial(n,m)*binomial(n, n-m) def g(m): # As given in A&S, but with `x` replaced by `phi(x)`. return ( ((phi(x) - 1)**(n-m)) * (phi(x) + 1)**m ) # From Abramowitz & Stegun, (22.3.2) with alpha = beta = 0. P = dn * sum([ c(m)*g(m) for m in range(0,n+1) ]) # If `x` is a symbolic expression, we want to return a symbolic # expression (even if that expression is e.g. `1`). if is_Expression(x): P = SR(P) return P