+
+
+
+-- | 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, fromList )
+--
+-- >>> 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 = matmap (\x -> (fromInteger 2)*x^2) first_components
+
+ weights = transpose $ weights_row