Changes between Initial Version and Version 27 of Ticket #8321


Ignore:
Timestamp:
Jul 1, 2011, 11:41:38 AM (12 years ago)
Author:
maldun
Comment:

Replying to 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.

Legend:

Unmodified
Added
Removed
Modified
  • Ticket #8321

    • Property Status changed from new to needs_work
    • Property Authors changed from to Stefan Reiterer
    • Property Cc maldun fredrik.johansson kcrisman mariah added
    • Property Keywords numerics integration added