From: Michael Orlitzky Date: Tue, 4 Feb 2014 22:41:19 +0000 (-0500) Subject: Add the Integration.Gaussian module where we begin to implement Gaussian quadrature. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=cf94f2c074a13ca9fd011fae5c0a668afefd3a18;p=numerical-analysis.git Add the Integration.Gaussian module where we begin to implement Gaussian quadrature. --- diff --git a/.ghci b/.ghci index 942d73d..e6f7087 100644 --- a/.ghci +++ b/.ghci @@ -4,6 +4,7 @@ -- Load everything. :{ :load src/BigFloat.hs + src/Integration/Gaussian.hs src/Integration/Simpson.hs src/Integration/Trapezoid.hs src/Linear/Iteration.hs @@ -19,6 +20,7 @@ :} import BigFloat +import Integration.Gaussian import Integration.Simpson import Integration.Trapezoid import Linear.Iteration diff --git a/numerical-analysis.cabal b/numerical-analysis.cabal index fe47ead..92b716e 100644 --- a/numerical-analysis.cabal +++ b/numerical-analysis.cabal @@ -18,6 +18,7 @@ data-files: makefile library exposed-modules: + Integration.Gaussian, Integration.Simpson, Integration.Trapezoid, Linear.Iteration, diff --git a/src/Integration/Gaussian.hs b/src/Integration/Gaussian.hs new file mode 100644 index 0000000..3169b20 --- /dev/null +++ b/src/Integration/Gaussian.hs @@ -0,0 +1,59 @@ +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE ScopedTypeVariables #-} + +-- | Gaussian quadrature over [-1, 1] using the Legendre polynomials. +-- +module Integration.Gaussian ( + jacobi_matrix ) +where + +import Algebra.Absolute ( abs ) +import qualified Algebra.Algebraic as Algebraic ( C ) +import qualified Algebra.Ring as Ring ( C ) +import Data.Vector.Fixed ( Arity ) +import NumericPrelude hiding ( abs ) + +import Linear.Matrix ( Mat(..), construct ) +import Linear.QR ( eigenvectors_symmetric ) + + +-- | Compute the Jacobi matrix for the Legendre polynomials over +-- [-1,1]. This matrix is derived from the three-term recurrence +-- relation which is expressable as a matrix equation. When the +-- polynomials are ortho/normal/, its Jacobi matrix should be +-- symmetric. But since the Legendre polynomials are not +-- ortho/normal/, we have to make a few modifications. The procedure +-- is described in Amparo, Gil, and Segura, \"Numerical Methods for +-- Special Functions\". +-- +-- Examples: +-- +-- >>> import Linear.Matrix ( Mat2, frobenius_norm, fromList ) +-- +-- >>> let c = 1/(sqrt 3) +-- >>> let expected = fromList [[0,c],[c, 0]] :: Mat2 Double +-- >>> let actual = jacobi_matrix :: Mat2 Double +-- >>> frobenius_norm (actual - expected) < 1e-14 +-- True +-- +jacobi_matrix :: forall m a. (Arity m, Algebraic.C a, Ring.C a) + => Mat m m a +jacobi_matrix = + construct lambda + where + coeff_a :: Int -> a + coeff_a n = fromRational' $ (toRational (n + 1))/(toRational (2*n + 1)) + + coeff_c :: Int -> a + coeff_c n = fromRational' $ (toRational n)/(toRational (2*n + 1)) + + max2 :: Ord b => b -> b -> b + max2 x y = if x >= y then x else y + + lambda i j + | j == i = fromInteger 0 + | abs (i-j) == 1 = + let k = max2 i j + in sqrt $ (coeff_a (k-1)) * (coeff_c k) + + | otherwise = fromInteger 0