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   gaussian,
   9   gaussian',
  10   jacobi_matrix,
  11   nodes_and_weights )
  12 where
  13
  14 import Algebra.Absolute ( abs )
  15 import qualified Algebra.Absolute as Absolute ( C )
  16 import qualified Algebra.Algebraic as Algebraic ( C )
  17 import qualified Algebra.Field as Field ( C )
  18 import qualified Algebra.Ring as Ring ( C )
  19 import qualified Algebra.ToRational as ToRational ( C )
  20 import Data.Vector.Fixed ( Arity, S )
  21 import NumericPrelude hiding ( abs )
  22
  23 import Linear.Matrix (
  24   Col,
  25   Col10,
  26   Mat(..),
  27   construct,
  28   element_sum2,
  29   fromList,
  30   identity_matrix,
  31   map2,
  32   row,
  33   transpose,
  34   zipwith2 )
  35 import Linear.QR ( eigenvectors_symmetric )
  36
  37
  38 -- | Compute the Jacobi matrix for the Legendre polynomials over
  39 --   [-1,1]. This matrix is derived from the three-term recurrence
  40 --   relation which is expressable as a matrix equation. When the
  41 --   polynomials are ortho/normal/, its Jacobi matrix should be
  42 --   symmetric. But since the Legendre polynomials are not
  43 --   ortho/normal/, we have to make a few modifications. The procedure
  44 --   is described in Amparo, Gil, and Segura, \"Numerical Methods for
  45 --   Special Functions\".
  46 --
  47 --   Examples:
  48 --
  49 --   >>> import Linear.Matrix ( Mat2, frobenius_norm )
  50 --
  51 --   >>> let c = 1/(sqrt 3)
  52 --   >>> let expected = fromList [[0,c],[c, 0]] :: Mat2 Double
  53 --   >>> let actual = jacobi_matrix :: Mat2 Double
  54 --   >>> frobenius_norm (actual - expected) < 1e-14
  55 --   True
  56 --
  57 jacobi_matrix :: forall m a. (Arity m, Algebraic.C a, Ring.C a)
  58               => Mat (S m) (S m) a
  59 jacobi_matrix =
  60   construct lambda
  61   where
  62     coeff_a :: Int -> a
  63     coeff_a n = fromRational' $ (toRational (n + 1))/(toRational (2*n + 1))
  64
  65     coeff_c :: Int -> a
  66     coeff_c n = fromRational' $ (toRational n)/(toRational (2*n + 1))
  67
  68     max2 :: Ord b => b -> b -> b
  69     max2 x y = if x >= y then x else y
  70
  71     lambda i j
  72       | j == i = fromInteger 0
  73       | abs (i-j) == 1 =
  74           let k = max2 i j
  75           in  sqrt $ (coeff_a (k-1)) * (coeff_c k)
  76
  77       | otherwise = fromInteger 0
  78
  79
  80
  81 -- | Compute the nodes/weights for Gaussian quadrature over [-1,1]
  82 --   using the Legendre polynomials. The test cases can all be found
  83 --   on e.g. Wikipedia.
  84 --
  85 --   Examples:
  86 --
  87 --   >>> import Linear.Matrix ( Col2, frobenius_norm )
  88 --
  89 --   >>> let (ns,ws) = nodes_and_weights 1000 :: (Col2 Double, Col2 Double)
  90 --   >>> let c = 1/(sqrt 3)
  91 --   >>> let expected_ns = fromList [[-c],[c]] :: Col2 Double
  92 --   >>> let expected_ws = fromList [[1],[1]] :: Col2 Double
  93 --   >>> frobenius_norm (ns - expected_ns) < 1e-12
  94 --   True
  95 --   >>> frobenius_norm (ws - expected_ws) < 1e-12
  96 --   True
  97 --
  98 nodes_and_weights :: forall m a.
  99                      (Arity m, Eq a, Algebraic.C a)
 100                   => Int
 101                   -> (Col (S m) a, Col (S m) a)
 102 nodes_and_weights iterations =
 103   (nodes, weights)
 104   where
 105     -- A shift is needed to make sure the QR algorithm
 106     -- converges. Since the roots (and thus eigenvalues) of orthogonal
 107     -- polynomials are real and distinct, there does exist a shift
 108     -- that will make it work. We guess lambda=1 makes it work.
 109
 110     shifted_jac = jacobi_matrix - identity_matrix :: Mat (S m) (S m) a
 111     (shifted_nodes, vecs) = eigenvectors_symmetric iterations shifted_jac
 112
 113     ones :: Col (S m) a
 114     ones = construct (\_ _ -> fromInteger 1)
 115
 116     nodes :: Col (S m) a
 117     nodes = shifted_nodes + ones -- unshift the nodes
 118
 119     -- Get the first component of each column.
 120     first_components = row vecs 0
 121
 122     -- Square it and multiply by 2; see the Golub-Welsch paper for
 123     -- this magic.
 124     weights_row = map2 (\x -> (fromInteger 2)*x^2) first_components
 125
 126     weights = transpose $ weights_row
 127
 128
 129 -- | Integrate a function over [-1,1] using Gaussian quadrature. This
 130 --   is the \"simple\" version of the 'gaussian'' function; here the ten
 131 --   nodes and weights are fixed.
 132 --
 133 --   Examples:
 134 --
 135 --   >>> let f x = x * (sin x)
 136 --   >>> let actual = gaussian f :: Double
 137 --   >>> let expected = 0.6023373578795136
 138 --   >>> abs (actual - expected) < 1e-12
 139 --   True
 140 --
 141 gaussian :: forall a.
 142             (Absolute.C a, Algebraic.C a, ToRational.C a, Field.C a)
 143          => (a -> a) -- ^ The function to integrate.
 144          -> a
 145 gaussian f =
 146   gaussian' f nodes weights
 147   where
 148     coerce :: (ToRational.C b) => b -> a
 149     coerce = fromRational' . toRational
 150
 151     nodes :: Col10 a
 152     nodes = fromList $ map (map coerce) node_list
 153
 154     node_list :: [[Double]]
 155     node_list = [[-0.9739065285171726],
 156                  [-0.8650633666889886],
 157                  [-0.6794095682990235],
 158                  [-0.43339539412924943],
 159                  [-0.14887433898163027],
 160                  [0.14887433898163094],
 161                  [0.4333953941292481],
 162                  [0.6794095682990247],
 163                  [0.8650633666889846],
 164                  [0.9739065285171717]]
 165
 166     weights :: Col10 a
 167     weights = fromList $ map (map coerce) weight_list
 168
 169     weight_list :: [[Double]]
 170     weight_list =  [[0.06667134430868775],
 171                     [0.14945134915057925],
 172                     [0.21908636251598282],
 173                     [0.26926671930999607],
 174                     [0.295524224714752],
 175                     [0.295524224714754],
 176                     [0.26926671930999635],
 177                     [0.219086362515982],
 178                     [0.14945134915058048],
 179                     [0.06667134430868814]]
 180
 181
 182
 183 -- | Integrate a function over [-1,1] using Gaussian quadrature. The
 184 --   nodes and the weights must be given (due to the design of the
 185 --   library, the number of nodes/weights must be requested at the
 186 --   type level so this function couldn't just compute e.g. @n@ of
 187 --   them for you).
 188 --
 189 --   Examples:
 190 --
 191 --   >>> import Linear.Matrix ( Col5 )
 192 --
 193 --   >>> let (ns, ws) = nodes_and_weights 100 :: (Col5 Double, Col5 Double)
 194 --   >>> let f x = x * (sin x)
 195 --   >>> let actual = gaussian' f ns ws
 196 --   >>> let expected = 0.6023373578795136
 197 --   >>> abs (actual - expected) < 1e-8
 198 --   True
 199 --
 200 gaussian' :: forall m a.
 201              (Arity m, ToRational.C a, Ring.C a)
 202           => (a -> a)    -- ^ The function @f@ to integrate.
 203           -> Col (S m) a -- ^ Column matrix of nodes
 204           -> Col (S m) a -- ^ Column matrix of weights
 205           -> a
 206 gaussian' f nodes weights =
 207   -- The one norm is just the sum of the entries, which is what we
 208   -- want.
 209   element_sum2 weighted_values
 210   where
 211     function_values = map2 f nodes
 212     weighted_values = zipwith2 (*) weights function_values