definite integral and indefinite integral different (~Gaussian)

Hello, I get strange behaviour on the following function. It persists when changing constants, but is not present for very simple functions, so any hint for a limit would be appreciated.

f1 = exp(-1/2*(y-a)^2/N)+exp(-1/2*(y+a)^2/N); f1

f2 = f1*y^2; f2

f4 = f2(a=3,N=0.7); f4
>(e^(-0.714285714285714*(y - 3)^2) + e^(-0.714285714285714*(y + 3)^2))*y^2

Z = integral(f4,y,-100,100); Z.n(digits=5)
>9.7000

F = integrate(f4,y); H = F(y=100) - F(y=-100); H.n(digits=5)
>-2.4622e-2917


As you see, I get right results when using the definite integral, while calculation the indefinite integral and manually evaluating it gains wrong results.

By the way, when integrating from -inf to inf, I should get N+a^2. Is there any way to see this?

edit retag close merge delete

What version of Sage are you using?

( 2012-08-12 08:49:49 -0500 )edit

Thank you for the answers, both of you!

( 2012-08-12 10:39:25 -0500 )edit

Sort by ยป oldest newest most voted

One can obtain an exact value of the integral (with factor a^2+N):

var('y a N')
assume(N>0)
I1=integrate(exp(-1/2*(y-a)^2/N)*y^2,(y,-oo,oo))
I2=integrate(exp(-1/2*(y+a)^2/N)*y^2,(y,-oo,oo))
I=(I1+I2).factor()
print I
print I(a=3.0,N=0.7).n()
#2*(a^2 + N)*sqrt(pi)*sqrt(2)*sqrt(N)
#40.6855961680184

more

I'm getting 40.686 for Z here.

Your integrals depend on the error function, and I know there are some issues with erf() in Sage. You can see some of the anomalies in the post here.

This could be the issue.

more

numerical_integral(f4,-oo,oo) gives: (40.685596168018435, 5.067965214014384e-06)

( 2012-08-12 10:37:58 -0500 )edit

Oh: you are correct, I also get 40.686. The 9.7 is when I have a different step f3 between f2 and f4.

( 2012-08-12 12:07:36 -0500 )edit