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   fromList,
  29   identity_matrix,
  30   map2,
  31   row,
  32   transpose,
  33   zipwith2 )
  34 import Linear.QR ( eigenvectors_symmetric )
  35 import Normed ( Normed(..) )
  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     jac = jacobi_matrix
 106
 107     -- A shift is needed to make sure the QR algorithm
 108     -- converges. Since the roots (and thus eigenvalues) of orthogonal
 109     -- polynomials are real and distinct, there does exist a shift
 110     -- that will make it work. We guess lambda=1 makes it work.
 111
 112     shifted_jac = jac - identity_matrix
 113     (shifted_nodes, vecs) = eigenvectors_symmetric iterations shifted_jac
 114
 115     ones :: Col (S m) a
 116     ones = construct (\_ _ -> fromInteger 1)
 117
 118     nodes :: Col (S m) a
 119     nodes = shifted_nodes + ones -- unshift the nodes
 120
 121     -- Get the first component of each column.
 122     first_components = row vecs 0
 123
 124     -- Square it and multiply by 2; see the Golub-Welsch paper for
 125     -- this magic.
 126     weights_row = map2 (\x -> (fromInteger 2)*x^2) first_components
 127
 128     weights = transpose $ weights_row
 129
 130
 131 -- | Integrate a function over [-1,1] using Gaussian quadrature. This
 132 --   is the \"simple\" version of the 'gaussian'' function; here the ten
 133 --   nodes and weights are fixed.
 134 --
 135 --   Examples:
 136 --
 137 --   >>> let f x = x * (sin x)
 138 --   >>> let actual = gaussian f :: Double
 139 --   >>> let expected = 0.6023373578795136
 140 --   >>> abs (actual - expected) < 1e-12
 141 --   True
 142 --
 143 gaussian :: forall a.
 144             (Absolute.C a, Algebraic.C a, ToRational.C a, Field.C a)
 145          => (a -> a) -- ^ The function to integrate.
 146          -> a
 147 gaussian f =
 148   gaussian' f nodes weights
 149   where
 150     coerce = fromRational' . toRational
 151
 152     nodes :: Col10 a
 153     nodes = fromList $ map (map coerce) node_list
 154
 155     node_list :: [[Double]]
 156     node_list = [[-0.9739065285171726],
 157                  [-0.8650633666889886],
 158                  [-0.6794095682990235],
 159                  [-0.43339539412924943],
 160                  [-0.14887433898163027],
 161                  [0.14887433898163094],
 162                  [0.4333953941292481],
 163                  [0.6794095682990247],
 164                  [0.8650633666889846],
 165                  [0.9739065285171717]]
 166
 167     weights :: Col10 a
 168     weights = fromList $ map (map coerce) weight_list
 169
 170     weight_list :: [[Double]]
 171     weight_list =  [[0.06667134430868775],
 172                     [0.14945134915057925],
 173                     [0.21908636251598282],
 174                     [0.26926671930999607],
 175                     [0.295524224714752],
 176                     [0.295524224714754],
 177                     [0.26926671930999635],
 178                     [0.219086362515982],
 179                     [0.14945134915058048],
 180                     [0.06667134430868814]]
 181
 182
 183
 184 -- | Integrate a function over [-1,1] using Gaussian quadrature. The
 185 --   nodes and the weights must be given (due to the design of the
 186 --   library, the number of nodes/weights must be requested at the
 187 --   type level so this function couldn't just compute e.g. @n@ of
 188 --   them for you).
 189 --
 190 --   The class constraints on @a@ could be loosened significantly if
 191 --   we implemented column sum outside of the 1-norm.
 192 --
 193 --   Examples:
 194 --
 195 --   >>> import Linear.Matrix ( Col5 )
 196 --
 197 --   >>> let (ns, ws) = nodes_and_weights 100 :: (Col5 Double, Col5 Double)
 198 --   >>> let f x = x * (sin x)
 199 --   >>> let actual = gaussian' f ns ws
 200 --   >>> let expected = 0.6023373578795136
 201 --   >>> abs (actual - expected) < 1e-8
 202 --   True
 203 --
 204 gaussian' :: forall m a.
 205              (Arity m, Absolute.C a, Algebraic.C a, ToRational.C a, Ring.C a)
 206           => (a -> a)    -- ^ The function @f@ to integrate.
 207           -> Col (S m) a -- ^ Column matrix of nodes
 208           -> Col (S m) a -- ^ Column matrix of weights
 209           -> a
 210 gaussian' f nodes weights =
 211   -- The one norm is just the sum of the entries, which is what we
 212   -- want.
 213   norm_p (1::Int) weighted_values
 214   where
 215     function_values = map2 f nodes
 216     weighted_values = zipwith2 (*) weights function_values