Add the Integration.Gaussian module where we begin to implement Gaussian quadrature.

[numerical-analysis.git] / src / Integration / Gaussian.hs

diff --git a/src/Integration/Gaussian.hs b/src/Integration/Gaussian.hs
new file mode 100644 (file)
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