src/Integration/Gaussian.hs

   1 {-# LANGUAGE FlexibleContexts #-}
   2 {-# LANGUAGE NoImplicitPrelude #-}
   3 {-# LANGUAGE ScopedTypeVariables #-}
   4
   5 -- | Gaussian quadrature over [-1, 1] using the Legendre polynomials.
   6 --
   7 module Integration.Gaussian (
   8   jacobi_matrix,
   9   nodes_and_weights )
  10 where
  11
  12 import Algebra.Absolute ( abs )
  13 import qualified Algebra.Algebraic as Algebraic ( C )
  14 import qualified Algebra.Ring as Ring ( C )
  15 import Data.Vector.Fixed ( Arity, S )
  16 import NumericPrelude hiding ( abs )
  17
  18 import Linear.Matrix (
  19   Col,
  20   Mat(..),
  21   construct,
  22   identity_matrix,
  23   matmap,
  24   row',
  25   transpose )
  26 import Linear.QR ( eigenvectors_symmetric )
  27
  28
  29 -- | Compute the Jacobi matrix for the Legendre polynomials over
  30 --   [-1,1]. This matrix is derived from the three-term recurrence
  31 --   relation which is expressable as a matrix equation. When the
  32 --   polynomials are ortho/normal/, its Jacobi matrix should be
  33 --   symmetric. But since the Legendre polynomials are not
  34 --   ortho/normal/, we have to make a few modifications. The procedure
  35 --   is described in Amparo, Gil, and Segura, \"Numerical Methods for
  36 --   Special Functions\".
  37 --
  38 --   Examples:
  39 --
  40 --   >>> import Linear.Matrix ( Mat2, frobenius_norm, fromList )
  41 --
  42 --   >>> let c = 1/(sqrt 3)
  43 --   >>> let expected = fromList [[0,c],[c, 0]] :: Mat2 Double
  44 --   >>> let actual = jacobi_matrix :: Mat2 Double
  45 --   >>> frobenius_norm (actual - expected) < 1e-14
  46 --   True
  47 --
  48 jacobi_matrix :: forall m a. (Arity m, Algebraic.C a, Ring.C a)
  49               => Mat (S m) (S m) a
  50 jacobi_matrix =
  51   construct lambda
  52   where
  53     coeff_a :: Int -> a
  54     coeff_a n = fromRational' $ (toRational (n + 1))/(toRational (2*n + 1))
  55
  56     coeff_c :: Int -> a
  57     coeff_c n = fromRational' $ (toRational n)/(toRational (2*n + 1))
  58
  59     max2 :: Ord b => b -> b -> b
  60     max2 x y = if x >= y then x else y
  61
  62     lambda i j
  63       | j == i = fromInteger 0
  64       | abs (i-j) == 1 =
  65           let k = max2 i j
  66           in  sqrt $ (coeff_a (k-1)) * (coeff_c k)
  67
  68       | otherwise = fromInteger 0
  69
  70
  71
  72 -- | Compute the nodes/weights for Gaussian quadrature over [-1,1]
  73 --   using the Legendre polynomials. The test cases can all be found
  74 --   on e.g. Wikipedia.
  75 --
  76 --   Examples:
  77 --
  78 --   >>> import Linear.Matrix ( Col2, frobenius_norm, fromList )
  79 --
  80 --   >>> let (ns,ws) = nodes_and_weights 1000 :: (Col2 Double, Col2 Double)
  81 --   >>> let c = 1/(sqrt 3)
  82 --   >>> let expected_ns = fromList [[-c],[c]] :: Col2 Double
  83 --   >>> let expected_ws = fromList [[1],[1]] :: Col2 Double
  84 --   >>> frobenius_norm (ns - expected_ns) < 1e-12
  85 --   True
  86 --   >>> frobenius_norm (ws - expected_ws) < 1e-12
  87 --   True
  88 --
  89 nodes_and_weights :: forall m a.
  90                      (Arity m, Eq a, Algebraic.C a)
  91                   => Int
  92                   -> (Col (S m) a, Col (S m) a)
  93 nodes_and_weights iterations =
  94   (nodes, weights)
  95   where
  96     jac = jacobi_matrix
  97
  98     -- A shift is needed to make sure the QR algorithm
  99     -- converges. Since the roots (and thus eigenvalues) of orthogonal
 100     -- polynomials are real and distinct, there does exist a shift
 101     -- that will make it work. We guess lambda=1 makes it work.
 102
 103     shifted_jac = jac - identity_matrix
 104     (shifted_nodes, vecs) = eigenvectors_symmetric iterations shifted_jac
 105
 106     ones :: Col (S m) a
 107     ones = construct (\_ _ -> fromInteger 1)
 108
 109     nodes :: Col (S m) a
 110     nodes = shifted_nodes + ones -- unshift the nodes
 111
 112     -- Get the first component of each column.
 113     first_components = row' vecs 0
 114
 115     -- Square it and multiply by 2; see the Golub-Welsch paper for
 116     -- this magic.
 117     weights_row = matmap (\x -> (fromInteger 2)*x^2) first_components
 118
 119     weights = transpose $ weights_row