src/Integration/Gaussian.hs

   1 {-# LANGUAGE NoImplicitPrelude #-}
   2 {-# LANGUAGE ScopedTypeVariables #-}
   3
   4 -- | Gaussian quadrature over [-1, 1] using the Legendre polynomials.
   5 --
   6 module Integration.Gaussian (
   7   jacobi_matrix )
   8 where
   9
  10 import Algebra.Absolute ( abs )
  11 import qualified Algebra.Algebraic as Algebraic ( C )
  12 import qualified Algebra.Ring as Ring ( C )
  13 import Data.Vector.Fixed ( Arity )
  14 import NumericPrelude hiding ( abs )
  15
  16 import Linear.Matrix ( Mat(..), construct )
  17 import Linear.QR ( eigenvectors_symmetric )
  18
  19
  20 -- | Compute the Jacobi matrix for the Legendre polynomials over
  21 --   [-1,1]. This matrix is derived from the three-term recurrence
  22 --   relation which is expressable as a matrix equation. When the
  23 --   polynomials are ortho/normal/, its Jacobi matrix should be
  24 --   symmetric. But since the Legendre polynomials are not
  25 --   ortho/normal/, we have to make a few modifications. The procedure
  26 --   is described in Amparo, Gil, and Segura, \"Numerical Methods for
  27 --   Special Functions\".
  28 --
  29 --   Examples:
  30 --
  31 --   >>> import Linear.Matrix ( Mat2, frobenius_norm, fromList )
  32 --
  33 --   >>> let c = 1/(sqrt 3)
  34 --   >>> let expected = fromList [[0,c],[c, 0]] :: Mat2 Double
  35 --   >>> let actual = jacobi_matrix :: Mat2 Double
  36 --   >>> frobenius_norm (actual - expected) < 1e-14
  37 --   True
  38 --
  39 jacobi_matrix :: forall m a. (Arity m, Algebraic.C a, Ring.C a)
  40               => Mat m m a
  41 jacobi_matrix =
  42   construct lambda
  43   where
  44     coeff_a :: Int -> a
  45     coeff_a n = fromRational' $ (toRational (n + 1))/(toRational (2*n + 1))
  46
  47     coeff_c :: Int -> a
  48     coeff_c n = fromRational' $ (toRational n)/(toRational (2*n + 1))
  49
  50     max2 :: Ord b => b -> b -> b
  51     max2 x y = if x >= y then x else y
  52
  53     lambda i j
  54       | j == i = fromInteger 0
  55       | abs (i-j) == 1 =
  56           let k = max2 i j
  57           in  sqrt $ (coeff_a (k-1)) * (coeff_c k)
  58
  59       | otherwise = fromInteger 0