Message Boards Message Boards

[✓] Avoid getting easy integral wrong?

GROUPS:

Seems impossible and so either we have user error (most likely) or Mathematica has an integration bug. Here's the integral:

 Integrate[(1 - Exp[(-k)*x^a])/(1 - Exp[-k]), {x, 0, 1}, Assumptions -> {k > 0, a > 0}]

And here's the answer I get.

 (E^k*(a*k^a^(-1) - Gamma[a^(-1)] + Gamma[a^(-1), k]))/(a*(-1 + E^k)*k^a^(-1))

Let's evaluate this expression at a=0.1 and k=0.1 We get zero! But this can't be. If we plot the function, we see that it is non-negative on the interval from 0 to 1 and assumes positive values in lots of places. And if I do numerical integration on the original expression, I get 0.912835, which looks about right.

What's going on? Is there a fix? Am I missing something?

By the way, there are other integrals similar to this, which also yield results that do not seem possible. Here's the integral on which I discovered what sure looks like a problem. I believe this integral should never be greater than 1/2. And, yet, when I evaluate it for various values of k and a I get values well over 1/2.

Integrate[(1 - Exp[(-k)*x^a])/(1 - Exp[-k]) - x, {x, 0, 1}, Assumptions -> {k > 0, a > 0}]
POSTED BY: Seth Chandler
Answer
4 months ago

It's floating point cancellation, we can use arbitrary precision or exact values:

In[7]:= (E^k k^(-1/a) (a k^(1/a) - Gamma[1/a] + Gamma[1/a, k]))/(
 a (-1 + E^k)) /. {a -> 1/10, k -> 1/10}

Out[7]= (100000000000 (-(36287999999999999/100000000000) + 
   401044422751291/(1000000000 E^(1/10))) E^(1/10))/(-1 + E^(1/10))

In[8]:= N[%, 20]

Out[8]= 0.91283470967626027646

In[9]:= N[%7]

Out[9]= 0.
POSTED BY: Ilian Gachevski
Answer
4 months ago

Group Abstract Group Abstract