From cf94f2c074a13ca9fd011fae5c0a668afefd3a18 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Tue, 4 Feb 2014 17:41:19 -0500
Subject: [PATCH] Add the Integration.Gaussian module where we begin to
 implement Gaussian quadrature.

---
 .ghci                       |  2 ++
 numerical-analysis.cabal    |  1 +
 src/Integration/Gaussian.hs | 59 +++++++++++++++++++++++++++++++++++++
 3 files changed, 62 insertions(+)
 create mode 100644 src/Integration/Gaussian.hs

diff --git a/.ghci b/.ghci
index 942d73d..e6f7087 100644
--- a/.ghci
+++ b/.ghci
@@ -4,6 +4,7 @@
 -- Load everything.
 :{
 :load src/BigFloat.hs
+  src/Integration/Gaussian.hs
   src/Integration/Simpson.hs
   src/Integration/Trapezoid.hs
   src/Linear/Iteration.hs
@@ -19,6 +20,7 @@
 :}
 
 import BigFloat
+import Integration.Gaussian
 import Integration.Simpson
 import Integration.Trapezoid
 import Linear.Iteration
diff --git a/numerical-analysis.cabal b/numerical-analysis.cabal
index fe47ead..92b716e 100644
--- a/numerical-analysis.cabal
+++ b/numerical-analysis.cabal
@@ -18,6 +18,7 @@ data-files: makefile
 
 library
   exposed-modules:
+    Integration.Gaussian,
     Integration.Simpson,
     Integration.Trapezoid,
     Linear.Iteration,
diff --git a/src/Integration/Gaussian.hs b/src/Integration/Gaussian.hs
new file mode 100644
index 0000000..3169b20
--- /dev/null
+++ b/src/Integration/Gaussian.hs
@@ -0,0 +1,59 @@
+{-# LANGUAGE NoImplicitPrelude #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+-- | Gaussian quadrature over [-1, 1] using the Legendre polynomials.
+--
+module Integration.Gaussian (
+  jacobi_matrix )
+where
+
+import Algebra.Absolute ( abs )
+import qualified Algebra.Algebraic as Algebraic ( C )
+import qualified Algebra.Ring as Ring ( C )
+import Data.Vector.Fixed ( Arity )
+import NumericPrelude hiding ( abs )
+
+import Linear.Matrix ( Mat(..), construct )
+import Linear.QR ( eigenvectors_symmetric )
+
+
+-- | Compute the Jacobi matrix for the Legendre polynomials over
+--   [-1,1]. This matrix is derived from the three-term recurrence
+--   relation which is expressable as a matrix equation. When the
+--   polynomials are ortho/normal/, its Jacobi matrix should be
+--   symmetric. But since the Legendre polynomials are not
+--   ortho/normal/, we have to make a few modifications. The procedure
+--   is described in Amparo, Gil, and Segura, \"Numerical Methods for
+--   Special Functions\".
+--
+--   Examples:
+--
+--   >>> import Linear.Matrix ( Mat2, frobenius_norm, fromList )
+--
+--   >>> let c = 1/(sqrt 3)
+--   >>> let expected = fromList [[0,c],[c, 0]] :: Mat2 Double
+--   >>> let actual = jacobi_matrix :: Mat2 Double
+--   >>> frobenius_norm (actual - expected) < 1e-14
+--   True
+--
+jacobi_matrix :: forall m a. (Arity m, Algebraic.C a, Ring.C a)
+              => Mat m m a
+jacobi_matrix =
+  construct lambda
+  where
+    coeff_a :: Int -> a
+    coeff_a n = fromRational' $ (toRational (n + 1))/(toRational (2*n + 1))
+
+    coeff_c :: Int -> a
+    coeff_c n = fromRational' $ (toRational n)/(toRational (2*n + 1))
+
+    max2 :: Ord b => b -> b -> b
+    max2 x y = if x >= y then x else y
+
+    lambda i j
+      | j == i = fromInteger 0
+      | abs (i-j) == 1 =
+          let k = max2 i j
+          in  sqrt $ (coeff_a (k-1)) * (coeff_c k)
+
+      | otherwise = fromInteger 0
-- 
2.43.2