Sage integral failing on the Watt Fission Spectrum
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.
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;
- Webscraping Roblox
- Remove the mandatory field label 'This field is required.' and fix the bug with 'clean_email'
- How to plot a Probability Density Function in Python?
- How large is a fresh install of Python?
- Appending new elements into an empty list
- Simple way to measure cell execution time in ipython notebook
- PyAudio working, but spits out error messages each time
- Reportlab show page number and page count IF there is more than one page in a document
- How to read SharePoint Online (Office365) Excel files into Python specifically pandas with Work or School Account?
- How to set a column which suffix name is based on a value in another column
- Debugging Python C++ extension from Visual Studio Code on Linux
- How can I get all users on Google admin_sdk?
- csv.Error: iterator should return strings, not bytes
- How to check if an object has an attribute?
- How to use selenium with proxy auth in headless mode?
- Is there a way to exit a pytest test and continue to the next one?
- Returning the lowest index for the first non whitespace character in a string in Python
- Formatting exceptions as Python does
- Prime factorization using list comprehension in Python
- Why does the power spectrum E(k) of my velocity field follow 𝑘 ^(−(n−1)) instead of 𝑘^(−n)?
- How to merge dataframes over multiple columns and split rows?
- How to create a Sympy IndexedBase using a custom subclass of Symbol?
- Removing dynamically an element from a list
- Returning boolean if set is empty
- Can variables be decorated?
- Fast(est) exponentiation of numpy 3D matrix
- Removing an element from a list based on a condition
- Printing elements of dictionary line by line
- Matplotlib does not display the hatch of a patch in a legend
- Python win32com - Class not registered error