1 {-# LANGUAGE FlexibleContexts #-}
2 {-# LANGUAGE NoImplicitPrelude #-}
3 {-# LANGUAGE ScopedTypeVariables #-}
5 -- | Gaussian quadrature over [-1, 1] using the Legendre polynomials.
7 module Integration.Gaussian (
12 import Algebra.Absolute ( abs )
13 import qualified Algebra.Algebraic as Algebraic ( C )
14 import qualified Algebra.Ring as Ring ( C )
15 import Data.Vector.Fixed ( Arity, S )
16 import NumericPrelude hiding ( abs )
18 import Linear.Matrix (
26 import Linear.QR ( eigenvectors_symmetric )
29 -- | Compute the Jacobi matrix for the Legendre polynomials over
30 -- [-1,1]. This matrix is derived from the three-term recurrence
31 -- relation which is expressable as a matrix equation. When the
32 -- polynomials are ortho/normal/, its Jacobi matrix should be
33 -- symmetric. But since the Legendre polynomials are not
34 -- ortho/normal/, we have to make a few modifications. The procedure
35 -- is described in Amparo, Gil, and Segura, \"Numerical Methods for
36 -- Special Functions\".
40 -- >>> import Linear.Matrix ( Mat2, frobenius_norm, fromList )
42 -- >>> let c = 1/(sqrt 3)
43 -- >>> let expected = fromList [[0,c],[c, 0]] :: Mat2 Double
44 -- >>> let actual = jacobi_matrix :: Mat2 Double
45 -- >>> frobenius_norm (actual - expected) < 1e-14
48 jacobi_matrix :: forall m a. (Arity m, Algebraic.C a, Ring.C a)
54 coeff_a n = fromRational' $ (toRational (n + 1))/(toRational (2*n + 1))
57 coeff_c n = fromRational' $ (toRational n)/(toRational (2*n + 1))
59 max2 :: Ord b => b -> b -> b
60 max2 x y = if x >= y then x else y
63 | j == i = fromInteger 0
66 in sqrt $ (coeff_a (k-1)) * (coeff_c k)
68 | otherwise = fromInteger 0
72 -- | Compute the nodes/weights for Gaussian quadrature over [-1,1]
73 -- using the Legendre polynomials. The test cases can all be found
78 -- >>> import Linear.Matrix ( Col2, frobenius_norm, fromList )
80 -- >>> let (ns,ws) = nodes_and_weights 1000 :: (Col2 Double, Col2 Double)
81 -- >>> let c = 1/(sqrt 3)
82 -- >>> let expected_ns = fromList [[-c],[c]] :: Col2 Double
83 -- >>> let expected_ws = fromList [[1],[1]] :: Col2 Double
84 -- >>> frobenius_norm (ns - expected_ns) < 1e-12
86 -- >>> frobenius_norm (ws - expected_ws) < 1e-12
89 nodes_and_weights :: forall m a.
90 (Arity m, Eq a, Algebraic.C a)
92 -> (Col (S m) a, Col (S m) a)
93 nodes_and_weights iterations =
98 -- A shift is needed to make sure the QR algorithm
99 -- converges. Since the roots (and thus eigenvalues) of orthogonal
100 -- polynomials are real and distinct, there does exist a shift
101 -- that will make it work. We guess lambda=1 makes it work.
103 shifted_jac = jac - identity_matrix
104 (shifted_nodes, vecs) = eigenvectors_symmetric iterations shifted_jac
107 ones = construct (\_ _ -> fromInteger 1)
110 nodes = shifted_nodes + ones -- unshift the nodes
112 -- Get the first component of each column.
113 first_components = row' vecs 0
115 -- Square it and multiply by 2; see the Golub-Welsch paper for
117 weights_row = matmap (\x -> (fromInteger 2)*x^2) first_components
119 weights = transpose $ weights_row