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_jac = jacobi_matrix - identity_matrix :: Mat (S m) (S m) a
(shifted_nodes, vecs) = eigenvectors_symmetric iterations shifted_jac
ones :: Col (S m) a
gaussian f =
gaussian' f nodes weights
where
+ coerce :: (ToRational.C b) => b -> a
coerce = fromRational' . toRational
nodes :: Col10 a