From: SourceForge.net <no...@so...> - 2011-05-24 16:22:08
|
Bugs item #3292033, was opened at 2011-04-23 11:24 Message generated for change (Comment added) made by dgildea You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=104933&aid=3292033&group_id=4933 Please note that this message will contain a full copy of the comment thread, including the initial issue submission, for this request, not just the latest update. Category: Lisp Core - Integration Group: None Status: Open Resolution: None Priority: 5 Private: No Submitted By: https://www.google.com/accounts () Assigned to: Nobody/Anonymous (nobody) Summary: error in integrating exp(-x)*sinh(sqrt(x)) Initial Comment: Maxima 5.24.0 http://maxima.sourceforge.net using Lisp SBCL 1.0.24 (%i3) integrate(exp(-x)*sinh(sqrt(x)),x,0,inf); (%o3) 0 (%i4) quad_qagi(exp(-x)*sinh(sqrt(x)),x,0,inf); (%o4) [1.137937897234377, 5.171862937913829e-11, 345, 0] ---------------------------------------------------------------------- >Comment By: Dan Gildea (dgildea) Date: 2011-05-24 12:22 Message: A simpler test case: (%i16) integrate(exp(sqrt(x)-x), x, 0, inf); (%o16) %e^(1/4)*gamma_incomplete(1,1/4)-%e^(1/4)*gamma_incomplete(1/2,1/4)/2 +2*%e^(1/4)*sqrt(%pi) (%i17) float(%); (%o17) 5.006110228172447 (%i18) quad_qagi(exp(sqrt(x)-x), x, 0, inf); (%o18) [2.730234433703704,9.207745677031198e-11,345,0] The problem seems to be caused by the fact that the indefinite integral: (%i15) integrate(exp(sqrt(x)-x), x); (%o15) %e^(1/4)*%i *(%i*gamma_incomplete(1,(1-2*sqrt(x))^2/4)*(1-2*sqrt(x))^2 /(2*sqrt(x)-1)^2 -%i*gamma_incomplete(1/2,(1-2*sqrt(x))^2/4)*(1-2*sqrt(x)) /(2*abs(2*sqrt(x)-1))) has a discontinuity at x=1/4. This integral is computed by case m2-exp-type-5 in integrate-exp-special. Mathematica gives the indefinite integral as: -exp(sqrt(x)-x)+exp(1/4)*sqrt(%pi)*erf((-1+2*sqrt(x))/2)/2; which results in the correct value for the definite integral in this case. ---------------------------------------------------------------------- Comment By: Aleksas (aleksasd) Date: 2011-05-18 14:20 Message: solving with Maxima 5.24.0: Dwight Tables of integrals 862.04 (%i1) declare(integrate,linear)$ assume(y>0)$ (%i3) S:'integrate(exp(-x)*sinh(sqrt(x)),x,0,inf)$ (%i4) exponentialize(%); (%o4) integrate((%e^sqrt(x)-%e^(-sqrt(x)))*%e^(-x),x,0,inf)/2 (%i5) expand(%); (%o5) integrate(%e^(sqrt(x)-x),x,0,inf)/2-integrate(%e^(-x-sqrt(x)),x,0,inf)/2 (%i6) S1:changevar(%, y=sqrt(x), y, x); (%o6) integrate(y*%e^(y-y^2),y,0,inf)-integrate(y*%e^(-y^2-y),y,0,inf) (%i10) changevar(part(S1,1),2*y-1=z,z,y)+ changevar(part(S1,2),2*y+1=z,z,y); (%o10) integrate((z+1)*%e^(-(z^2-1)/4),z,-1,inf)/4-integrate((z-1)*%e^(-(z^2-1)/4),z,1,inf)/4 (%i8) ev(%, nouns),ratsimp; (%o8) (%e^(1/4)*sqrt(%pi))/2 (%i9) float(%), numer; (%o9) 1.137937897234373 Where is the error? wrong: (%i1) R:integrate(y*%e^(y-y^2),y,0,inf); (%o1) (%e^(1/4)*gamma_incomplete(1,1/4))/2-(%e^(1/4)*gamma_incomplete(1/2,1/4))/4 (%i2) float(R), numer; (%o2) 0.22717931961748 (%i3) quad_qagi(y*%e^(y-y^2), y, 0, inf); (%o3) [1.365117216851849,7.6857545173630158*10^-10,165,0] correct: (%i4) R1:integrate(y*%e^(-y-y^2),y,0,inf); (%o4) (%e^(1/4)*gamma_incomplete(1,1/4))/2-(%e^(1/4)*gamma_incomplete(1/2,1/4))/4 (%i5) float(R1), numer; (%o5) 0.22717931961748 (%i6) quad_qagi(y*%e^(-y-y^2), y, 0, inf); (%o6) [0.22717931961748,1.4118104693665376*10^-9,135,0] Aleksas D ---------------------------------------------------------------------- You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=104933&aid=3292033&group_id=4933 |