From: Michael Orlitzky <michael@orlitzky.com>
Date: Wed, 5 Feb 2014 20:30:47 +0000 (-0500)
Subject: Implement Gaussian quadrature and test it.
X-Git-Url: http://gitweb.michael.orlitzky.com/?p=numerical-analysis.git;a=commitdiff_plain;h=6f237bfcd5890513cd390000c1ac8cdbdbaffb8c

Implement Gaussian quadrature and test it.
---

diff --git a/src/Integration/Gaussian.hs b/src/Integration/Gaussian.hs
index 5d683f6..8fcca4f 100644
--- a/src/Integration/Gaussian.hs
+++ b/src/Integration/Gaussian.hs
@@ -5,25 +5,34 @@
 -- | Gaussian quadrature over [-1, 1] using the Legendre polynomials.
 --
 module Integration.Gaussian (
+  gaussian,
+  gaussian',
   jacobi_matrix,
   nodes_and_weights )
 where
 
 import Algebra.Absolute ( abs )
+import qualified Algebra.Absolute as Absolute ( C )
 import qualified Algebra.Algebraic as Algebraic ( C )
+import qualified Algebra.Field as Field ( C )
 import qualified Algebra.Ring as Ring ( C )
+import qualified Algebra.ToRational as ToRational ( C )
 import Data.Vector.Fixed ( Arity, S )
 import NumericPrelude hiding ( abs )
 
 import Linear.Matrix (
   Col,
+  Col10,
   Mat(..),
+  colzipwith,
   construct,
+  fromList,
   identity_matrix,
   matmap,
   row',
   transpose )
 import Linear.QR ( eigenvectors_symmetric )
+import Normed ( Normed(..) )
 
 
 -- | Compute the Jacobi matrix for the Legendre polynomials over
@@ -37,7 +46,7 @@ import Linear.QR ( eigenvectors_symmetric )
 --
 --   Examples:
 --
---   >>> import Linear.Matrix ( Mat2, frobenius_norm, fromList )
+--   >>> import Linear.Matrix ( Mat2, frobenius_norm )
 --
 --   >>> let c = 1/(sqrt 3)
 --   >>> let expected = fromList [[0,c],[c, 0]] :: Mat2 Double
@@ -75,7 +84,7 @@ jacobi_matrix =
 --
 --   Examples:
 --
---   >>> import Linear.Matrix ( Col2, frobenius_norm, fromList )
+--   >>> import Linear.Matrix ( Col2, frobenius_norm )
 --
 --   >>> let (ns,ws) = nodes_and_weights 1000 :: (Col2 Double, Col2 Double)
 --   >>> let c = 1/(sqrt 3)
@@ -117,3 +126,91 @@ nodes_and_weights iterations =
     weights_row = matmap (\x -> (fromInteger 2)*x^2) first_components
 
     weights = transpose $ weights_row
+
+
+-- | Integrate a function over [-1,1] using Gaussian quadrature. This
+--   is the \"simple\" version of the 'gaussian'' function; here the ten
+--   nodes and weights are fixed.
+--
+--   Examples:
+--
+--   >>> let f x = x * (sin x)
+--   >>> let actual = gaussian f :: Double
+--   >>> let expected = 0.6023373578795136
+--   >>> abs (actual - expected) < 1e-12
+--   True
+--
+gaussian :: forall a.
+            (Absolute.C a, Algebraic.C a, ToRational.C a, Field.C a)
+         => (a -> a) -- ^ The function to integrate.
+         -> a
+gaussian f =
+  gaussian' f nodes weights
+  where
+    coerce = fromRational' . toRational
+
+    nodes :: Col10 a
+    nodes = fromList $ map (map coerce) node_list
+
+    node_list :: [[Double]]
+    node_list = [[-0.9739065285171726],
+                 [-0.8650633666889886],
+                 [-0.6794095682990235],
+                 [-0.43339539412924943],
+                 [-0.14887433898163027],
+                 [0.14887433898163094],
+                 [0.4333953941292481],
+                 [0.6794095682990247],
+                 [0.8650633666889846],
+                 [0.9739065285171717]]
+
+    weights :: Col10 a
+    weights = fromList $ map (map coerce) weight_list
+
+    weight_list :: [[Double]]
+    weight_list =  [[0.06667134430868775],
+                    [0.14945134915057925],
+                    [0.21908636251598282],
+                    [0.26926671930999607],
+                    [0.295524224714752],
+                    [0.295524224714754],
+                    [0.26926671930999635],
+                    [0.219086362515982],
+                    [0.14945134915058048],
+                    [0.06667134430868814]]
+
+
+
+-- | Integrate a function over [-1,1] using Gaussian quadrature. The
+--   nodes and the weights must be given (due to the design of the
+--   library, the number of nodes/weights must be requested at the
+--   type level so this function couldn't just compute e.g. @n@ of
+--   them for you).
+--
+--   The class constraints on @a@ could be loosened significantly if
+--   we implemented column sum outside of the 1-norm.
+--
+--   Examples:
+--
+--   >>> import Linear.Matrix ( Col5 )
+--
+--   >>> let (ns, ws) = nodes_and_weights 100 :: (Col5 Double, Col5 Double)
+--   >>> let f x = x * (sin x)
+--   >>> let actual = gaussian' f ns ws
+--   >>> let expected = 0.6023373578795136
+--   >>> abs (actual - expected) < 1e-8
+--   True
+--
+gaussian' :: forall m a.
+             (Arity m, Absolute.C a, Algebraic.C a, ToRational.C a, Ring.C a)
+          => (a -> a)    -- ^ The function @f@ to integrate.
+          -> Col (S m) a -- ^ Column matrix of nodes
+          -> Col (S m) a -- ^ Column matrix of weights
+          -> a
+gaussian' f nodes weights =
+  -- The one norm is just the sum of the entries, which is what we
+  -- want.
+  norm_p (1::Int) weighted_values
+  where
+    function_values = matmap f nodes
+    weighted_values = colzipwith (*) weights function_values