[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