[Scipy-tickets] [SciPy] #1466: More accurate roots/weights for Gauss-Hermite quadrature

SciPy Trac scipy-tickets@scipy....
Tue Jun 28 01:02:32 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:                 |  
---------------------------+------------------------------------------------

Comment(by charris):

 You are at about at the machine precision limit and the roots aren't
 symmetrical, but a single pass of Newton cleans things up nicely:


 {{{
 In [51]: from numpy.polynomial import Hermite as H

 In [52]: p = H([0]*40 + [1])

 In [53]: r = p.roots()

 In [54]: r1 = r - p(r)/(p.deriv()(r))

 In [55]: r1
 Out[55]:
 array([-8.0987611392508505 , -7.41158253148546908, -6.84023730524935569,
        -6.32825535122008187, -5.85409505603039992, -5.40665424797012761,
        -4.97926097854525551, -4.56750207284439469, -4.16825706683250008,
        -3.77920675343522339, -3.3985582658596285 , -3.02487988390128448,
        -2.65699599844289569, -2.29391714187508322, -1.93479147228229587,
        -1.578869894931614  , -1.22548010904628901, -0.8740066123570881 ,
        -0.52387471383227713, -0.1745372145975824 ,  0.1745372145975824 ,
         0.52387471383227713,  0.8740066123570881 ,  1.22548010904628901,
         1.578869894931614  ,  1.93479147228229564,  2.29391714187508322,
         2.65699599844289569,  3.02487988390128448,  3.3985582658596285 ,
         3.77920675343522339,  4.16825706683250008,  4.56750207284439469,
         4.97926097854525551,  5.40665424797012761,  5.85409505603039992,
         6.32825535122008187,  6.84023730524935569,  7.41158253148546908,
         8.0987611392508505 ])


 }}}



 {{{
 Recursive roots:
 [-8.0987611392508505  -7.41158253148546908 -6.84023730524935569
  -6.32825535122008187 -5.85409505603039992 -5.40665424797012761
  -4.97926097854525551 -4.56750207284439469 -4.16825706683250008
  -3.77920675343522339 -3.3985582658596285  -3.02487988390128448
  -2.65699599844289569 -2.29391714187508367 -1.93479147228229564
  -1.578869894931614   -1.22548010904628901 -0.87400661235708799
  -0.52387471383227724 -0.1745372145975824   0.1745372145975824
   0.52387471383227724  0.87400661235708799  1.22548010904628901
   1.578869894931614    1.93479147228229564  2.29391714187508367
   2.65699599844289569  3.02487988390128448  3.3985582658596285
   3.77920675343522339  4.16825706683250008  4.56750207284439469
   4.97926097854525551  5.40665424797012761  5.85409505603039992
   6.32825535122008187  6.84023730524935569  7.41158253148546908
   8.0987611392508505 ]
 }}}

-- 
Ticket URL: <http://projects.scipy.org/scipy/ticket/1466#comment:1>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.


More information about the Scipy-tickets mailing list