Hi,
I want to have a simple script to obtain the Gaussian quadrature
points and weights.
I started with points, that are just roots of Legendre polynomials [1, 2]:
In [2]: p = legendre_poly(64, x)
In [4]: r1 = nsolve(p, 0.0243)
This gives the root close to "0.0243". However, as can be seen below, it's not
as precise as it could be:
In [5]: r2 = Float("0.024350292663424433")
In [8]: p.subs(x, r1)
Out[8]: -1.05002809488412e-16
In [9]: p.subs(x, r2)
Out[9]: -7.85866634120666e-18
The r2 is clearly much more accurate. Does anyone know how to increase
the accuracy?
I tried things like:
In [10]: import sympy.mpmath
In [15]: sympy.mpmath.mp.dps
Out[15]: 15
In [16]: sympy.mpmath.mp.dps = 200
In [17]: r1 = nsolve(p, 0.0243)
In [20]: r1
Out[20]:
0.0243502926634244496317653613537851156236027440792911104610944223230041344547
657467006373545311269982730883817820285820001244657144755607813544705987728320
14115658100725011616718802395041644889394938962
In [21]: "%e" % p.subs(x, r1)
Out[21]: -1.097728e-16
or
In [22]: r1 = nsolve(p, 0.0243, tol=1e-100)
In [23]: r1
Out[23]:
0.0243502926634244496317653613537851156236027440792911104610944223230041344547
657467006373545311269982730883817820285820001244657144755607813544705987728320
14115658100725011616718802395041644889394938962
In [24]: "%e" % p.subs(x, r1)
Out[24]: -1.097728e-16
But no luck.
Even with using 200 digits, I still get a root which is only 1e-16 accurate.
Ondrej
[1] http://en.wikipedia.org/wiki/Gaussian_quadrature
[2]
http://processingjs.nihongoresources.com/bezierinfo/legendre-gauss-values.php
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to sympy@googlegroups.com.
To unsubscribe from this group, send email to
sympy+unsubscr...@googlegroups.com.
For more options, visit this group at
http://groups.google.com/group/sympy?hl=en.
