X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Forthogonal_polynomials.py;h=15c695611fb3e287c357075f45e55eeddf550f7f;hb=1bbade9f41ffbfe366b15d0db657f666bc1f025d;hp=0da947cd6f0bf21705f501ec4523868135ff658f;hpb=c649291bc271a68a94e172f1946abffb6e5d9fdd;p=sage.d.git

diff --git a/mjo/orthogonal_polynomials.py b/mjo/orthogonal_polynomials.py
index 0da947c..15c6956 100644
--- a/mjo/orthogonal_polynomials.py
+++ b/mjo/orthogonal_polynomials.py
@@ -1,7 +1,7 @@
 from sage.all import *
 
 def legendre_p(n, x, a = -1, b = 1):
-    """
+    r"""
     Returns the ``n``th Legendre polynomial of the first kind over the
     interval [a, b] with respect to ``x``.
 
@@ -28,6 +28,10 @@ def legendre_p(n, x, a = -1, b = 1):
     returned. Otherwise, the value of the ``n``th polynomial at ``x``
     will be returned.
 
+    SETUP::
+
+        sage: from mjo.orthogonal_polynomials import legendre_p
+
     EXAMPLES:
 
     Create the standard Legendre polynomials in `x`::
@@ -60,15 +64,32 @@ def legendre_p(n, x, a = -1, b = 1):
 
     And finite field elements::
 
-        sage: legendre_P(3, GF(11)(5))
+        sage: legendre_p(3, GF(11)(5))
         8
 
+    Solve a simple least squares problem over `[-\pi, \pi]`::
+
+        sage: a = -pi
+        sage: b = pi
+        sage: def inner_product(v1, v2):
+        ....:     return integrate(v1*v2, x, a, b)
+        sage: def norm(v):
+        ....:     return sqrt(inner_product(v,v))
+        sage: def project(basis, v):
+        ....:     return sum( inner_product(v, b)*b/norm(b)**2
+        ....:                 for b in basis)
+        sage: f = sin(x)
+        sage: legendre_basis = [ legendre_p(k, x, a, b) for k in range(4) ]
+        sage: proj = project(legendre_basis, f)
+        sage: proj.simplify_trig()
+        5/2*(7*(pi^2 - 15)*x^3 - 3*(pi^4 - 21*pi^2)*x)/pi^6
+
     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
+        sage: all( eq(k) for k in range(20) ) # long time
         True
 
     We can evaluate the result of the zeroth polynomial::
@@ -172,6 +193,6 @@ def legendre_p(n, x, a = -1, b = 1):
 
     # From Abramowitz & Stegun, (22.3.2) with alpha = beta = 0.
     # Also massaged to support finite field elements.
-    P = sum([ c(m)*g(m) for m in range(0,n+1) ])/(2**n)
+    P = sum( c(m)*g(m) for m in range(n+1) )/(2**n)
 
     return P