From: Michael Orlitzky <michael@orlitzky.com>
Date: Thu, 30 Jan 2014 21:02:38 +0000 (-0500)
Subject: Add the Polynomials.Orthogonal module.
X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=edbf6900d50865b4b6278f13cd75d55df1289c90;p=numerical-analysis.git

Add the Polynomials.Orthogonal module.
---

diff --git a/.ghci b/.ghci
index 38ef02d..b8019ba 100644
--- a/.ghci
+++ b/.ghci
@@ -13,6 +13,7 @@
   src/Misc.hs
   src/Normed.hs
   src/ODE/IVP.hs
+  src/Polynomials/Orthogonal.hs
   src/Roots/Simple.hs
 :}
 
@@ -26,6 +27,7 @@ import Linear.Vector
 import Misc
 import Normed
 import ODE.IVP
+import Polynomials.Orthogonal
 import Roots.Simple
 
 -- Use a calmer prompt.
diff --git a/makefile b/makefile
index 0ffe8d3..a333c58 100644
--- a/makefile
+++ b/makefile
@@ -7,23 +7,19 @@ SRCS = $(shell find src/ -name '*.hs')
 
 .PHONY : test publish_doc doc dist hlint
 
-$(BIN): $(SRCS)
-	runghc Setup.hs clean
-	runghc Setup.hs configure --user --flags=${FLAGS}
+$(BIN): $(PN).cabal $(SRCS)
+	runghc Setup.hs configure --user
 	runghc Setup.hs build
 
 $(DOCTESTS_BIN): $(SRCS) test/Doctests.hs
-	runghc Setup.hs configure --user --flags=${FLAGS} --enable-tests
+	runghc Setup.hs configure --user --enable-tests
 	runghc Setup.hs build
 
 
-profile: $(SRCS)
+profile: $(PN).cabal $(SRCS)
 	runghc Setup.hs configure --user --enable-executable-profiling
 	runghc Setup.hs build
 
-hpc: $(SRCS)
-	FLAGS="hpc" make
-
 clean:
 	runghc Setup.hs clean
 	rm -f dist/
@@ -34,17 +30,12 @@ clean:
 test: $(BIN) $(DOCTESTS_BIN)
 	runghc Setup.hs test
 
-dist:
-	runghc Setup.hs configure
-	runghc Setup.hs sdist
-
 
 doc: $(SRCS)
 	runghc Setup.hs configure --user --flags=${FLAGS}
-	runghc Setup.hs hscolour --executables
-	runghc Setup.hs haddock --internal    \
-				--executables \
-				--hyperlink-source
+	runghc Setup.hs hscolour --all
+	runghc Setup.hs haddock  --all \
+                                 --haddock-options="--ignore-all-exports"
 
 hlint:
 	hlint --ignore="Use camelCase"     \
diff --git a/numerical-analysis.cabal b/numerical-analysis.cabal
index 29dd006..2bc7468 100644
--- a/numerical-analysis.cabal
+++ b/numerical-analysis.cabal
@@ -19,7 +19,8 @@ data-files: makefile
 library
   exposed-modules: Integration.Simpson, Integration.Trapezoid,
     Linear.Iteration, Linear.Matrix, Linear.System, Linear.Vector,
-    Misc, Normed, ODE.IVP, Roots.Simple, Roots.Fast
+    Misc, Normed, ODE.IVP, Polynomials.Orthogonal, Roots.Simple,
+    Roots.Fast
 
   build-depends:
     base              >= 3 && < 5,
diff --git a/src/Polynomials/Orthogonal.hs b/src/Polynomials/Orthogonal.hs
new file mode 100644
index 0000000..28cf41d
--- /dev/null
+++ b/src/Polynomials/Orthogonal.hs
@@ -0,0 +1,75 @@
+-- | The polynomials over [a,b] form a basis for L_2[a,b]. But often
+--   the \"obvious\" choice of a basis {1, x, x^2,...} is not
+--   convenient, because its elements are not orthogonal.
+--
+--   There are several sets of polynomials that are known to be
+--   orthogonal over various intervals. Here we provide formulae for a
+--   few popular sets (bases).
+--
+--   First we present the simple formulae that generates polynomials
+--   orthogonal over a specific interval; say, [-1, 1] (these are what
+--   you'll find in a textbook). Then, we present a \"prime\" version
+--   that produces polynomials orthogonal over an arbitrary interval
+--   [a, b]. This is accomplished via an affine transform which shifts
+--   and scales the variable but preserves orthogonality.
+--
+--   The choice of basis depends not only on the interval, but on the
+--   weight function. The inner product used is not necessarily the
+--   standard one in L_2, rather,
+--
+--     \<f, g\> = \int_{a}^{b} w(x)*f(x)*g(x) dx
+--
+--   where w(x) is some non-negative (or non-positive) weight function.
+--
+module Polynomials.Orthogonal
+where
+
+import NumericPrelude
+import qualified Algebra.RealField as RealField
+import qualified Prelude as P
+
+
+-- | The @n@th Legendre polynomial in @x@ over [-1,1]. These are
+--   orthogonal over [-1,1] with the weight function w(x)=1.
+--
+--   Examples:
+--
+--   >>> let points = [0, 1, 2, 3, 4, 5] :: [Double]
+--
+--   The zeroth polynomial is the constant function p(x)=1.
+--
+--   >>> map (legendre 0) points
+--   [1.0,1.0,1.0,1.0,1.0,1.0]
+--
+--   The first polynomial is the identity:
+--
+--   >>> map (legendre 1) points
+--   [0.0,1.0,2.0,3.0,4.0,5.0]
+--
+--   The second polynomial should be (3x^2 - 1)/2:
+--
+--   >>> let actual = map (legendre 2) points
+--   >>> let f x = (3*x^2 - 1)/2
+--   >>> let expected = map f points
+--   >>> actual == expected
+--   True
+--
+--   The third polynomial should be (5x^3 - 3x)/2:
+--
+--   >>> let actual = map (legendre 3) points
+--   >>> let f x = (5*x^3 - 3*x)/2
+--   >>> let expected = map f points
+--   >>> actual == expected
+--   True
+--
+legendre :: (RealField.C a)
+         => Integer -- ^ The degree (i.e. the number of) the polynomial you want.
+         -> a -- ^ The dependent variable @x@.
+         -> a
+legendre 0 _ = fromInteger 1
+legendre 1 x = x
+legendre n x =
+  (c1*x * (legendre (n-1) x) - c2*(legendre (n-2) x)) / (fromInteger n)
+  where
+    c1 = fromInteger $ 2*n - 1
+    c2 = fromInteger $ n - 1