[Numpy-discussion] numpy quad and maple
davidmenhur at gmail.com
Sat May 10 13:43:57 EDT 2014
On 10 May 2014 01:04, Ihab Riad <ifriad at gmail.com> wrote:
> I did the following integral with numpy quad
> quad(lambda h:
> and I get the following values (1.8368139214123403e-126,
The first value is the integral, the second is the error. You have (1.8 +-
3.3) 10^-126. That is, your integral is almost zero.
> When I do the integral with maple and mathematica I get
That is inside the bracket Scipy said. And I wouldn't even trust them to be
correct without doing some checking.
Numbers so small are numerically non trustworthy. Your function is:
E^(A - B log(h) - C log^2(h)) = E^A * h^B' * h ^ (C' log(h))
where B' and C' depend on your numbers. E^A is already 10^-507 (are you
sure this numbers are correct?), so if you strip it out of your integral,
you have a simpler expression, and more reasonable values. This is what is
killing your integral and sending it to 0.
Now, you have to normalise the values. So, changing variables: h-> sx, dh
sx^B' * (sx)^[C' (log(s) + log(x))]
Expand and simplify, find the value that makes more reasonable values
(order of one).
Last thing: giving quad a limit in infinity can make it behave crazily. If
you look at your second term, that is roughly x^(-log(x)), it decays very
quickly for large values of x (10^-5 for x=20), so you can, without
affecting your integral much, set a cut at a some reasonable h (plot it to
be sure). You can even make an estimate comparing with the analytical
result  of that particular piece.
Bottomline: don't use a lambda, use a function and expand your parameters.
That will make the expression simpler.
log10 = np.log10(h)
B = 1936.9610182100055*5
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the NumPy-Discussion