[Scipy-tickets] [SciPy] #1466: More accurate roots/weights for Gauss-Hermite quadrature
SciPy Trac
scipy-tickets@scipy....
Mon Jun 27 23:28:03 CDT 2011
#1466: More accurate roots/weights for Gauss-Hermite quadrature
---------------------------+------------------------------------------------
Reporter: Bogdan | Owner: pv
Type: enhancement | Status: new
Priority: normal | Milestone: Unscheduled
Component: scipy.special | Version: devel
Keywords: |
---------------------------+------------------------------------------------
System: OSX 10.6.7, Python 2.6.1 (the one provided with system), Numpy
2.0.0.dev-3071eab and SciPy 0.10.0dev (installed via SciPy Superpack)
I am using scipy.special.orthogonal.hermite and .h_roots to perform
decomposition into basis of harmonic oscillator eigenfunctions (which is
essentially G-H quadrature). I noticed that decomposition is not behaving
well when the number of modes is big (>40 in my case). After some testing
I narrowed it down to inaccurate polynomial roots provided by h_roots().
As far as I understand from the code, scipy implementation uses Golub-
Welsch algorithm, which is general and fast because it employs linear
algebra. But in our imperfect world with finite precision numbers it
starts to lose stability with large matrix sizes. So I took recursive root
finding algorithm from Numerical Recipes (3rd ed., section 4.6.1) and it
showed much better accuracy (see attached script with comparison).
I understand that this algorithm is less general (because it uses good
guesses for first roots of Hermite polynomials) and definitely slower than
G-W, but there are applications when its accuracy is necessary. I think it
should be present in SciPy, at least as an alternative to the current
method.
--
Ticket URL: <http://projects.scipy.org/scipy/ticket/1466>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.
More information about the Scipy-tickets
mailing list