Sage integral failing on the Watt Fission Spectrum

Question

I am trying to integrate the Watt Fission Spectrum using Sage (the commands are below), but when I did I noticed that the result from the integral does not make sense. The Watt Fission Spectrum is a pdf and when integrated from 0 to 10 should be approximately one (confirmed by Wolfram Alpha). With Sage, when the lower bound is close to 0 the integral is negative, and then jumps positive when the lower bound is closer to 1.

Below is a plot of the integral of the Watt Fission Spectrum with the bounds x to 10 with x from 0 to 1. The jump is very clear.

How can I correct this, and be confident in the future that integrals are evaluated correctly?

var('o')
assume(x>10)
plot(integral(0.453*exp(-1.036*o)*sinh(sqrt(2.29*o)), (o, x, 10)), (x,0,1))

Per Robert Dodier's suggestion I have tried integration as follows:

map(var, 'abcde')
assume(e-d>0)
foo(a,b,c,d,e) = integral(a*exp(b*o)*sinh(sqrt(c*o)), (o, d, e))
float(foo(*map(QQ, [0.453,-1.036,2.29,0,10])))

This still returns a value of -0.0010615612261540579.

Answer 1

I don't know why Sage isn't relaying the questions from Maxima, but without those, you can't tell (in Sage) that the result is different depending on the parameters.

Here's a derivation in Maxima. I've rephrased the integral as being from 0 to x instead of x to 10. Since the integrand is a pdf, we should see the integral go from 0 to 1 as x increases.

Original formula.

(%i2) foo : 0.453*exp(-1.036*o)*sinh(sqrt(2.29*o)) $

Express that with symbols instead of numbers.

(%i3) foo1 : a * exp(-b*o) * sinh(sqrt(c*o)) $

Forestall some questions from integrate. It's not absolutely necessary to do this; integrate will just ask additional questions if these assumptions are not given.

(%i4) assume (equal (a, ratsimp (0.453)), equal (b, ratsimp (1.036)), equal (c, ratsimp (2.29))) $

rat: replaced 0.453 by 453/1000 = 0.453

rat: replaced 1.036 by 259/250 = 1.036

rat: replaced 2.29 by 229/100 = 2.29
(%i5) assume (x > 0) $

Compute the integral. integrate asks about x in relation to the parameters a, b, and c. I'll give each answer in turn and see what I get. I1, I2, and I3 are the results for all answers. I'll answer p (positive) first.

(%i6) I1 : integrate (foo1, o, 0, x);
Is 4*b^2*x-c positive, negative or zero?

p;
(%o6) a*(-gamma_incomplete(1/2,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b))
        *sqrt(c)*%e^(c/(4*b))*sqrt(x)
        /(2*sqrt(b)*abs(2*b*sqrt(x)-sqrt(c)))
        +gamma_incomplete(1/2,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b))
         *c*%e^(c/(4*b))
         /(4*b^(3/2)*abs(2*b*sqrt(x)-sqrt(c)))
        -gamma_incomplete(1/2,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b))
         *sqrt(c)*%e^(c/(4*b))
         /(4*b^(3/2))+sqrt(%pi)*sqrt(c)*%e^(c/(4*b))/(2*b^(3/2))
        +gamma_incomplete(1,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b))
         *%e^(c/(4*b))
         /(2*b)
        -gamma_incomplete(1,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b))
         *%e^(c/(4*b))
         /(2*b))

Now I'll answer n (negative).

(%i7) I2 : integrate (foo1, o, 0, x);
Is 4*b^2*x-c positive, negative or zero?

n;
(%o7) a*(-gamma_incomplete(1/2,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b))
        *sqrt(c)*%e^(c/(4*b))*sqrt(x)
        /(2*sqrt(b)*abs(2*b*sqrt(x)-sqrt(c)))
        +gamma_incomplete(1/2,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b))
         *c*%e^(c/(4*b))
         /(4*b^(3/2)*abs(2*b*sqrt(x)-sqrt(c)))
        -gamma_incomplete(1/2,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b))
         *sqrt(c)*%e^(c/(4*b))
         /(4*b^(3/2))
        +gamma_incomplete(1,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b))
         *%e^(c/(4*b))
         /(2*b)
        -gamma_incomplete(1,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b))
         *%e^(c/(4*b))
         /(2*b))

And now z (zero).

(%i8) I3 : integrate (foo1, o, 0, x);
Is 4*b^2*x-c positive, negative or zero?

z;
(%o8) a*(-gamma_incomplete(1/2,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b))
        *sqrt(c)*%e^(c/(4*b))*sqrt(x)
        /(2*sqrt(b)*abs(2*b*sqrt(x)-sqrt(c)))
        +gamma_incomplete(1/2,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b))
         *c*%e^(c/(4*b))
         /(4*b^(3/2)*abs(2*b*sqrt(x)-sqrt(c)))
        -gamma_incomplete(1/2,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b))
         *sqrt(c)*%e^(c/(4*b))
         /(4*b^(3/2))
        +gamma_incomplete(1,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b))
         *%e^(c/(4*b))
         /(2*b)
        -gamma_incomplete(1,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b))
         *%e^(c/(4*b))
         /(2*b))

I can't tell just by looking if, perhaps, some of these results are the same. Let's check.

(%i9) is (I1 = I2);
(%o9) false
(%i10) is (I1 = I3);
(%o10) false
(%i11) is (I2 = I3);
(%o11) true

OK, so there are two distinct results: I1 (answered p) and I2 (answered n and z). I'll put these two into a single conditional expression, in which the test is the question that integrate asked. Note that c/(4*b^2) is about 0.533, which appears to be the point at which the plot shown by OP has a jump.

(%i12) I : if x < c/(4*b^2) then ''I2 else ''I1;
(%o12) if x < c/(4*b^2)
           then a*(-gamma_incomplete(1/2,
                                     -(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b))
                  *sqrt(c)*%e^(c/(4*b))*sqrt(x)
                  /(2*sqrt(b)*abs(2*b*sqrt(x)-sqrt(c)))
                  +gamma_incomplete(1/2,
                                    -(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b))
                   *c*%e^(c/(4*b))
                   /(4*b^(3/2)*abs(2*b*sqrt(x)-sqrt(c)))
                  -gamma_incomplete(1/2,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b))
                   *sqrt(c)*%e^(c/(4*b))
                   /(4*b^(3/2))
                  +gamma_incomplete(1,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b))
                   *%e^(c/(4*b))
                   /(2*b)
                  -gamma_incomplete(1,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b))
                   *%e^(c/(4*b))
                   /(2*b))
           else a*(-gamma_incomplete(1/2,
                                     -(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b))
                  *sqrt(c)*%e^(c/(4*b))*sqrt(x)
                  /(2*sqrt(b)*abs(2*b*sqrt(x)-sqrt(c)))
                  +gamma_incomplete(1/2,
                                    -(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b))
                   *c*%e^(c/(4*b))
                   /(4*b^(3/2)*abs(2*b*sqrt(x)-sqrt(c)))
                  -gamma_incomplete(1/2,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b))
                   *sqrt(c)*%e^(c/(4*b))
                   /(4*b^(3/2))+sqrt(%pi)*sqrt(c)*%e^(c/(4*b))/(2*b^(3/2))
                  +gamma_incomplete(1,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b))
                   *%e^(c/(4*b))
                   /(2*b)
                  -gamma_incomplete(1,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b))
                   *%e^(c/(4*b))
                   /(2*b))

Now I'll plot that. plot2d is evaluated with a, b, and c temporarily bound to the specific values. Plot looks as expected.

(%i13) plot2d (I, [x, 0, 10]), a=0.453, b=1.036, c=2.29;