mjo/orthogonal_polynomials.py

   1 from sage.all import *
   2
   3 def legendre_p(n, x, a = -1, b = 1):
   4     """
   5     Returns the ``n``th Legendre polynomial of the first kind over the
   6     interval [a, b] with respect to ``x``.
   7
   8     When [a,b] is not [-1,1], we scale the standard Legendre
   9     polynomial (which is defined over [-1,1]) via an affine map. The
  10     resulting polynomials are still orthogonal and possess the
  11     property that `P(a) = P(b) = 1`.
  12
  13     INPUT:
  14
  15       * ``n`` -- The index of the polynomial.
  16
  17       * ``x`` -- Either the variable to use as the independent
  18         variable in the polynomial, or a point at which to evaluate
  19         the polynomial.
  20
  21       * ``a`` -- The "left" endpoint of the interval. Must be a real number.
  22
  23       * ``b`` -- The "right" endpoint of the interval. Must be a real number.
  24
  25     OUTPUT:
  26
  27     If ``x`` is a variable, a polynomial (symbolic expression) will be
  28     returned. Otherwise, the value of the ``n``th polynomial at ``x``
  29     will be returned.
  30
  31     SETUP::
  32
  33         sage: from mjo.orthogonal_polynomials import legendre_p
  34
  35     EXAMPLES:
  36
  37     Create the standard Legendre polynomials in `x`::
  38
  39         sage: legendre_p(0,x)
  40         1
  41         sage: legendre_p(1,x)
  42         x
  43
  44     Reuse the variable from a polynomial ring::
  45         sage: P.<t> = QQ[]
  46         sage: legendre_p(2,t)
  47         3/2*t^2 - 1/2
  48
  49     If ``x`` is a real number, the result should be as well::
  50
  51         sage: legendre_p(3, 1.1)
  52         1.67750000000000
  53
  54     Similarly for complex numbers::
  55
  56         sage: legendre_p(3, 1 + I)
  57         7/2*I - 13/2
  58
  59     Even matrices work::
  60
  61         sage: legendre_p(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
  62         [-179  242]
  63         [-484  547]
  64
  65     And finite field elements::
  66
  67         sage: legendre_p(3, GF(11)(5))
  68         8
  69
  70     Solve a simple least squares problem over `[-\pi, \pi]`::
  71
  72         sage: a = -pi
  73         sage: b = pi
  74         sage: def inner_product(v1, v2):
  75         ...       return integrate(v1*v2, x, a, b)
  76         ...
  77         sage: def norm(v):
  78         ...       return sqrt(inner_product(v,v))
  79         ...
  80         sage: def project(basis, v):
  81         ...       return sum([ inner_product(v, b)*b/norm(b)**2
  82         ...                    for b in basis])
  83         ...
  84         sage: f = sin(x)
  85         sage: legendre_basis = [ legendre_p(k, x, a, b) for k in range(0,4) ]
  86         sage: proj = project(legendre_basis, f)
  87         sage: proj.simplify_trig()
  88         5/2*(7*(pi^2 - 15)*x^3 - 3*(pi^4 - 21*pi^2)*x)/pi^6
  89
  90     TESTS:
  91
  92     We should agree with Maxima for all `n`::
  93
  94         sage: eq = lambda k: bool(legendre_p(k,x) == legendre_P(k,x))
  95         sage: all([eq(k) for k in range(0,20) ]) # long time
  96         True
  97
  98     We can evaluate the result of the zeroth polynomial::
  99
 100         sage: f = legendre_p(0,x)
 101         sage: f(x=10)
 102         1
 103
 104     We should have |P(a)| = |P(b)| = 1 for all a,b::
 105
 106         sage: a = RR.random_element()
 107         sage: b = RR.random_element()
 108         sage: k = ZZ.random_element(20)
 109         sage: P = legendre_p(k, x, a, b)
 110         sage: abs(P(x=a)) # abs tol 1e-12
 111         1
 112         sage: abs(P(x=b)) # abs tol 1e-12
 113         1
 114
 115     Two different polynomials should be orthogonal with respect to the
 116     inner product over `[a,b]`. Note that this test can fail if QQ is
 117     replaced with RR, since integrate() can return NaN::
 118
 119         sage: a = QQ.random_element()
 120         sage: b = QQ.random_element()
 121         sage: j = ZZ.random_element(20)
 122         sage: k = j + 1
 123         sage: Pj = legendre_p(j, x, a, b)
 124         sage: Pk = legendre_p(k, x, a, b)
 125         sage: integrate(Pj*Pk, x, a, b) # abs tol 1e-12
 126         0
 127
 128     The first few polynomials shifted to [0,1] are known to be::
 129
 130         sage: p0 = 1
 131         sage: p1 = 2*x - 1
 132         sage: p2 = 6*x^2 - 6*x + 1
 133         sage: p3 = 20*x^3 - 30*x^2 + 12*x - 1
 134         sage: bool(legendre_p(0, x, 0, 1) == p0)
 135         True
 136         sage: bool(legendre_p(1, x, 0, 1) == p1)
 137         True
 138         sage: bool(legendre_p(2, x, 0, 1) == p2)
 139         True
 140         sage: bool(legendre_p(3, x, 0, 1) == p3)
 141         True
 142
 143     The zeroth polynomial is an element of the ring that we're working
 144     in::
 145
 146         sage: legendre_p(0, MatrixSpace(ZZ, 2)([1, 2, -4, 7]))
 147         [1 0]
 148         [0 1]
 149
 150     """
 151     if not n in ZZ:
 152         raise TypeError('n must be a natural number')
 153
 154     if n < 0:
 155         raise ValueError('n must be nonnegative')
 156
 157     if not (a in RR and b in RR):
 158         raise TypeError('both `a` and `b` must be real numbers')
 159
 160     if n == 0:
 161         # Easy base case, save time. Attempt to return a value in the
 162         # same field/ring as `x`.
 163         return x.parent()(1)
 164
 165     # Even though we know a,b are real we use a different ring. We
 166     # prefer ZZ so that we can support division of finite field
 167     # elements by (a-b). Eventually this should be supported for QQ as
 168     # well, although it does not work at the moment. The preference of
 169     # SR over RR is to return something attractive when e.g. a=pi.
 170     if a in ZZ:
 171         a = ZZ(a)
 172     else:
 173         a = SR(a)
 174
 175     if b in ZZ:
 176         b = ZZ(b)
 177     else:
 178         b = SR(b)
 179
 180     # Ensure that (2**n) is an element of ZZ. This is used later --
 181     # we can divide finite field elements by integers but we can't
 182     # multiply them by rationals at the moment.
 183     n = ZZ(n)
 184
 185     def phi(t):
 186         # This is an affine map from [a,b] into [-1,1] and so
 187         # preserves orthogonality.
 188         return (a + b - 2*t)/(a - b)
 189
 190     def c(m):
 191         return binomial(n,m)*binomial(n, n-m)
 192
 193     def g(m):
 194         # As given in A&S, but with `x` replaced by `phi(x)`.
 195         return ( ((phi(x) - 1)**(n-m)) * (phi(x) + 1)**m )
 196
 197     # From Abramowitz & Stegun, (22.3.2) with alpha = beta = 0.
 198     # Also massaged to support finite field elements.
 199     P = sum([ c(m)*g(m) for m in range(0,n+1) ])/(2**n)
 200
 201     return P