X-Git-Url: http://gitweb.michael.orlitzky.com/?p=numerical-analysis.git;a=blobdiff_plain;f=src%2FIntegration%2FGaussian.hs;h=6eb46f6683a9c58d0105da1f1519d732c53bdd64;hp=3169b2091b67faad01d1c34bae375ee4e7c9e1d5;hb=46f18f6dd356db75b77bd95a15260c7c1a817bff;hpb=cf94f2c074a13ca9fd011fae5c0a668afefd3a18 diff --git a/src/Integration/Gaussian.hs b/src/Integration/Gaussian.hs index 3169b20..6eb46f6 100644 --- a/src/Integration/Gaussian.hs +++ b/src/Integration/Gaussian.hs @@ -1,20 +1,38 @@ +{-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE NoImplicitPrelude #-} {-# LANGUAGE ScopedTypeVariables #-} -- | Gaussian quadrature over [-1, 1] using the Legendre polynomials. -- module Integration.Gaussian ( - jacobi_matrix ) + 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 Data.Vector.Fixed ( Arity ) +import qualified Algebra.ToRational as ToRational ( C ) +import Data.Vector.Fixed ( Arity, S ) import NumericPrelude hiding ( abs ) -import Linear.Matrix ( Mat(..), construct ) +import Linear.Matrix ( + Col, + Col10, + Mat(..), + construct, + fromList, + identity_matrix, + map2, + row, + transpose, + zipwith2 ) import Linear.QR ( eigenvectors_symmetric ) +import Normed ( Normed(..) ) -- | Compute the Jacobi matrix for the Legendre polynomials over @@ -28,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 @@ -37,7 +55,7 @@ import Linear.QR ( eigenvectors_symmetric ) -- True -- jacobi_matrix :: forall m a. (Arity m, Algebraic.C a, Ring.C a) - => Mat m m a + => Mat (S m) (S m) a jacobi_matrix = construct lambda where @@ -57,3 +75,142 @@ jacobi_matrix = in sqrt $ (coeff_a (k-1)) * (coeff_c k) | otherwise = fromInteger 0 + + + +-- | Compute the nodes/weights for Gaussian quadrature over [-1,1] +-- using the Legendre polynomials. The test cases can all be found +-- on e.g. Wikipedia. +-- +-- Examples: +-- +-- >>> import Linear.Matrix ( Col2, frobenius_norm ) +-- +-- >>> let (ns,ws) = nodes_and_weights 1000 :: (Col2 Double, Col2 Double) +-- >>> let c = 1/(sqrt 3) +-- >>> let expected_ns = fromList [[-c],[c]] :: Col2 Double +-- >>> let expected_ws = fromList [[1],[1]] :: Col2 Double +-- >>> frobenius_norm (ns - expected_ns) < 1e-12 +-- True +-- >>> frobenius_norm (ws - expected_ws) < 1e-12 +-- True +-- +nodes_and_weights :: forall m a. + (Arity m, Eq a, Algebraic.C a) + => Int + -> (Col (S m) a, Col (S m) a) +nodes_and_weights iterations = + (nodes, weights) + where + jac = jacobi_matrix + + -- A shift is needed to make sure the QR algorithm + -- converges. Since the roots (and thus eigenvalues) of orthogonal + -- polynomials are real and distinct, there does exist a shift + -- that will make it work. We guess lambda=1 makes it work. + + shifted_jac = jac - identity_matrix + (shifted_nodes, vecs) = eigenvectors_symmetric iterations shifted_jac + + ones :: Col (S m) a + ones = construct (\_ _ -> fromInteger 1) + + nodes :: Col (S m) a + nodes = shifted_nodes + ones -- unshift the nodes + + -- Get the first component of each column. + first_components = row vecs 0 + + -- Square it and multiply by 2; see the Golub-Welsch paper for + -- this magic. + weights_row = map2 (\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 = map2 f nodes + weighted_values = zipwith2 (*) weights function_values