#8321: numerical integration with arbitrary precision
-------------------------+--------------------------------------------------
Reporter: burcin | Owner: burcin
Type: defect | Status: needs_work
Priority: major | Milestone: sage-4.7.1
Component: symbolics | Keywords: numerics,integration
Work_issues: | Upstream: N/A
Reviewer: | Author: Stefan Reiterer
Merged: | Dependencies:
-------------------------+--------------------------------------------------
Old description:
> From the sage-devel:
>
> {{{
> On Feb 20, 2010, at 12:40 PM, John H Palmieri wrote:
> ...
> > I was curious about this, so I tried specifying the number of digits:
> >
> > sage: h = integral(sin(x)/x^2, (x, 1, pi/2)); h
> > integrate(sin(x)/x^2, x, 1, 1/2*pi)
> > sage: h.n()
> > 0.33944794097891573
> > sage: h.n(digits=14)
> > 0.33944794097891573
> > sage: h.n(digits=600)
> > 0.33944794097891573
> > sage: h.n(digits=600) == h.n(digits=14)
> > True
> > sage: h.n(prec=50) == h.n(prec=1000)
> > True
> >
> > Is there an inherit limit in Sage on the accuracy of numerical
> > integrals?
> }}}
>
> The `_evalf_` function defined on line 179 of
> `sage/symbolic/integration/integral.py` calls the gsl
> `numerical_integral()` function and ignores the precision.
>
> We should raise a `NotImplementedError` for high precision, or find a way
> to do arbitrary precision numerical integration.
New description:
From the sage-devel:
{{{
On Feb 20, 2010, at 12:40 PM, John H Palmieri wrote:
...
> I was curious about this, so I tried specifying the number of digits:
>
> sage: h = integral(sin(x)/x^2, (x, 1, pi/2)); h
> integrate(sin(x)/x^2, x, 1, 1/2*pi)
> sage: h.n()
> 0.33944794097891573
> sage: h.n(digits=14)
> 0.33944794097891573
> sage: h.n(digits=600)
> 0.33944794097891573
> sage: h.n(digits=600) == h.n(digits=14)
> True
> sage: h.n(prec=50) == h.n(prec=1000)
> True
>
> Is there an inherit limit in Sage on the accuracy of numerical
> integrals?
}}}
The `_evalf_` function defined on line 179 of
`sage/symbolic/integration/integral.py` calls the gsl
`numerical_integral()` function and ignores the precision.
We should raise a `NotImplementedError` for high precision, or find a way
to do arbitrary precision numerical integration.
--
Comment(by maldun):
Replying to [comment:26 burcin]:
> I hadn't actually done the computation in Sage. I just wanted to note
that web site here. :)
>
> With the patch attached to this ticket:
>
> {{{
> sage: sigma = 1e-4
> sage: ff = exp(-x**2/(2*sigma)) / sqrt(2*pi*sigma)
> sage: from sage.symbolic.integration.integral import definite_integral
> sage: definite_integral(ff, x, -20, 10, hold=True)
> integrate(70.7106781186548*e^(-5000.00000000000*x^2)/sqrt(pi), x, -20,
10)
> sage: definite_integral(ff, x, -20, 10, hold=True).n()
> 5.38249970415053e-939
> }}}
>
> Without the patch:
>
> {{{
> sage: definite_integral(ff, x, -20, 10, hold=True).n()
> 2.1074458151264474e-45
> }}}
>
> We get a better result if we allow maxima to evaluate the integral
symbolically:
>
> {{{
> sage: definite_integral(ff, x, -20, 10).n()
> 1.00000000000000
> sage: definite_integral(ff, x, -20, 10)
> 0.353553390593*(sqrt(2)*sqrt(pi)*erf(500*sqrt(2)) +
sqrt(2)*sqrt(pi)*erf(1000*sqrt(2)))/sqrt(pi)
> }}}
>
This is indeed a great example!
This is again a grid problem and not a precision problem. If you make
trapezoidal or gauss rule
on finer grids you get also better results.
{{{
sage: from numpy import *
sage: from scipy.integrate import trapz
sage: sigma = 1e-4
sage: def ff(x): return exp(-x**2/(2*sigma))/sqrt(2*pi*sigma)
....:
sage: ffv = vectorize(ff)
sage: x = arange(-20,10,1)
sage: y = ffv(x)
sage: trapz(y,x)
39.894228040143268
sage: x = arange(-20,10,0.5)
sage: y = ffv(x)
sage: trapz(y,x)
19.947114020071634
sage: x = arange(-20,10,0.05)
sage: y = ffv(x)
sage: trapz(y,x)
1.9947262692023391
sage: x = arange(-20,10,0.005)
sage: y = ffv(x)
sage: trapz(y,x)
1.0
from scipy.integrate import fixed_quad
sage: fixed_quad(ff,-20,10,n=int(10))
(0.0, None)
sage: fixed_quad(ff,-20,10,n=int(100))
(2.6290056634068843e-58, None)
sage: fixed_quad(ff,-20,10,n=int(1000))
(0.8616135058547989, None)
}}}
The reason for this is simply that the function is
approximately 1 in a small region arround zero and nearly zero elsewhere.
Thats also the reason why maxima works so well here: The main part is
included in the analytical solution.
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/8321#comment:27>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en.
