Opened 13 years ago

Closed 13 years ago

# [with patch, positive review] 'integrate' produces incorrect answer

Reported by: Owned by: emchristiansen burcin major sage-4.2 calculus integrate, integral, incorrect sage-4.2.alpha0 Karl-Dieter Crisman Mike Hansen

### Description

Affects: x64 Ubuntu 9.04 (Jaunty) Sage 4.0, compiled from source, and updated (sage -upgrade) as of this posting

By simply casting a number using 'n()', I can cause integrate to return a vastly different result. See below:

var = sage.calculus.calculus.var

def NormalPDF(x,mu,sigma):
return 1/sqrt(2*pi*sigma^2)*exp(-1/(2*sigma**2)*(x-mu)^2)

x = var('x')
mu = var('m')
sigma = var('s')
assume(sigma>0)
child1 = NormalPDF(x,0,1)
child2 = NormalPDF(x,n(0),n(1))
parent = NormalPDF(x,mu,sigma)

# this expansion helps Sage to do the integral
integrand1 = ((parent-child1)^2).expand()
integrand2 = ((parent-child2)^2).expand()

int1 = integrate(integrand1,x,-infinity,infinity)
int2 = integrate(integrand2,x,-infinity,infinity)

print "THIS EXPRESSION:"
print int1
print "\nSHOULD BE VERY SIMILIAR TO THIS EXPRESSION:"
print int2

print "\nMAKING THIS NUMBER:"
print int1.subs({mu:0,sigma:1}).n()
print "\nVERY SIMILAR TO THIS NUMBER:"
print int2.subs({mu:0,sigma:1}).n()


The above produces the output:

THIS EXPRESSION:
1/2*((s + 1)*sqrt(s^2 + 1)*e^(1/2*m^2/(s^2 + 1)) - 2*sqrt(2)*s)*e^(-1/2*m^2/(s^2 + 1))/(sqrt(s^2 + 1)*sqrt(pi)*s)

SHOULD BE VERY SIMILIAR TO THIS EXPRESSION:
1/2*(sqrt(0.5*s^2 + 0.5)*e^(1/2*m^2/(s^2 + 1)) - sqrt(2)*s)*e^(-1/2*m^2/(s^2 + 1))/(sqrt(0.5*s^2 + 0.5)*sqrt(pi)*s)

MAKING THIS NUMBER:
0.000000000000000

VERY SIMILAR TO THIS NUMBER:
-0.116847488627555


Mathematica claims the correct integral is:

\frac{1+\sigma }{2 \sqrt{\pi } \sigma }-\frac{\sqrt{2} e^{-\frac{\mu ^2}{2+2 \sigma ^2}}}{\sqrt{\pi +\pi  \sigma ^2}}


### comment:1 Changed 13 years ago by kcrisman

• Summary changed from 'integrate' produces incorrect answer to [with patch, needs review] 'integrate' produces incorrect answer

This is now fixed, presumably in the Maxima upgrade. Note that the integral in fact computes without expand(), and in that case there is no 'experimental error'! Attached patch verifies this.

### Changed 13 years ago by kcrisman

Based on 4.1.2.alpha2

### comment:2 Changed 13 years ago by kcrisman

Despite what it says, actually based on 4.1.2.alpha4.

### comment:3 Changed 13 years ago by mhansen

• Authors set to Karl-Dieter Crisman
• Merged in set to sage-4.2.alpha0
• Reviewers set to Mike Hansen
• Status changed from needs_review to positive_review
• Summary changed from [with patch, needs review] 'integrate' produces incorrect answer to [with patch, positive review] 'integrate' produces incorrect answer

Looks good to me.

### comment:4 Changed 13 years ago by mhansen

• Milestone set to sage-4.2
• Resolution set to fixed
• Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.