defint or gamma_incomplete bug?



Numerically evaluating this expression results in subtractive cancellation; try:

(%i9) e : integrate(integrate(cos(x*y*log(2)),y,1,10),x,1,10)$

(%i10) makelist(block([fpprec : 15], bfloat(block([fpprec : 10 * k], expand(bfloat(e))))), k,1,10);

(%o10) [-2.238107472896b14,6.40685300762984b-16*%i+6.58088693446743b6,-9.59754643393143b-1,
        -9.5975928150589b-1,-9.59759281505811b-1,-9.59759281505811b-1,-9.59759281505811b-1,
        -9.59759281505811b-1,-9.59759281505811b-1,-9.59759281505811b-1]

The nfloat function (defined in hypergeometric) knows nothing about gamma_incomplete :(

--Barton

-----maxima-bounces at math.utexas.edu wrote: -----
To: "maxima mailing list" <maxima at math.utexas.edu>
From: "Edwin Woollett" 
Sent by: maxima-bounces at math.utexas.edu
Date: 12/20/2011 02:58PM
Subject: defint or gamma_incomplete bug?

integrate gives a wildly wrong answer for two
2d integrals:
--------------------------------------------------
Maxima 5.25.1 http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (a.k.a. GCL)

(%i1) quad_qag(quad_qag(cos(x*y*log(2)),y,1,10,3)[1],x,1,10,3);

(%o1) [-0.95975928150581,2.9857740100970423E-14,217,0]

(%i2) quad_qag(quad_qag(sin(x*y*log(2)),y,1,10,3)[1],x,1,10,3);

(%o2) [0.061314144038473,2.928537406327685E-14,217,0]

(%i3) integrate(integrate(cos(x*y*log(2)),y,1,10),x,1,10);

(%o3) -(%i*gamma_incomplete(0,100*%i*log(2))
?? ? ? -2*%i*gamma_incomplete(0,10*%i*log(2))+%i*gamma_incomplete(0,%i*log(2))
?? ? ? -%i*gamma_incomplete(0,-%i*log(2))
?? ? ? +2*%i*gamma_incomplete(0,-10*%i*log(2))
?? ? ? -%i*gamma_incomplete(0,-100*%i*log(2)))
?? ? ? /(2*log(2))

(%i4) expand(float(rectform(%)));

(%o4) 3.3117308623575085E+10

(%i5) integrate(integrate(sin(x*y*log(2)),y,1,10),x,1,10);

(%o5) (gamma_incomplete(0,100*%i*log(2))-2*gamma_incomplete(0,10*%i*log(2))
?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?+gamma_incomplete(0,%i*log(2))
?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?+gamma_incomplete(0,-%i*log(2))
?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?-2*gamma_incomplete(0,-10*%i*log(2))
?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?+gamma_incomplete(0,-100*%i*log(2)))
?? ? ? /(2*log(2))

(%i6) expand(float(rectform(%)));

(%o6) 3.298208192195924E+10
----------------------------------------------
Ted Woollett

_______________________________________________
Maxima mailing list
Maxima at math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima