1 {-# LANGUAGE NoImplicitPrelude #-}
2 {-# LANGUAGE ScopedTypeVariables #-}
4 -- | Gaussian quadrature over [-1, 1] using the Legendre polynomials.
6 module Integration.Gaussian (
10 import Algebra.Absolute ( abs )
11 import qualified Algebra.Algebraic as Algebraic ( C )
12 import qualified Algebra.Ring as Ring ( C )
13 import Data.Vector.Fixed ( Arity )
14 import NumericPrelude hiding ( abs )
16 import Linear.Matrix ( Mat(..), construct )
17 import Linear.QR ( eigenvectors_symmetric )
20 -- | Compute the Jacobi matrix for the Legendre polynomials over
21 -- [-1,1]. This matrix is derived from the three-term recurrence
22 -- relation which is expressable as a matrix equation. When the
23 -- polynomials are ortho/normal/, its Jacobi matrix should be
24 -- symmetric. But since the Legendre polynomials are not
25 -- ortho/normal/, we have to make a few modifications. The procedure
26 -- is described in Amparo, Gil, and Segura, \"Numerical Methods for
27 -- Special Functions\".
31 -- >>> import Linear.Matrix ( Mat2, frobenius_norm, fromList )
33 -- >>> let c = 1/(sqrt 3)
34 -- >>> let expected = fromList [[0,c],[c, 0]] :: Mat2 Double
35 -- >>> let actual = jacobi_matrix :: Mat2 Double
36 -- >>> frobenius_norm (actual - expected) < 1e-14
39 jacobi_matrix :: forall m a. (Arity m, Algebraic.C a, Ring.C a)
45 coeff_a n = fromRational' $ (toRational (n + 1))/(toRational (2*n + 1))
48 coeff_c n = fromRational' $ (toRational n)/(toRational (2*n + 1))
50 max2 :: Ord b => b -> b -> b
51 max2 x y = if x >= y then x else y
54 | j == i = fromInteger 0
57 in sqrt $ (coeff_a (k-1)) * (coeff_c k)
59 | otherwise = fromInteger 0