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

diff --git a/.ghci b/.ghci
index 942d73d9664e999e85e28246f843e969bb221468..e6f7087deaf56dbde83cdd76a9ec7d74e475b2e0 100644 (file)
--- 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 fe47eadb0c46d0b13b44e96653384f8132773a5c..92b716ef46b809f2f8fd3dc8b3b36c11e241ca69 100644 (file)
--- 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 (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