Re: AW: Re: Nasty bug in Integrate (version 5.0)
- To: mathgroup at smc.vnet.net
- Subject: [mg45980] Re: AW: [mg45958] Re: Nasty bug in Integrate (version 5.0)
- From: Dr Bob <drbob at bigfoot.com>
- Date: Sat, 31 Jan 2004 05:20:50 -0500 (EST)
- References: <CGEPLKMMCCOGKJGLHEKNIEBECBAA.klamser@t-online.de>
- Reply-to: drbob at bigfoot.com
- Sender: owner-wri-mathgroup at wolfram.com
Of course there's no error-free software, but did I miss something? You
talked about proving results are errorfree, but all you did is show a
couple of answers. Which is correct, and based on what? Visual inspection?
(If that's always enough, why do we need Mathematica?)
To check an integration, you need to differentiate, so let's try that in
5.0.1, shall we? First, here's the integrand, and we integrate it (I'm not
showing the result).
f[x_] = 1/x + x^c;
fiveOhOne = Integrate[f[x], {x, a, b}];
The derivative of fiveOhOne at x==a should be -f[a]. The derivative at
x==b should be f[b]. So let's try at b, for now. Is the derivative correct
at b?
FullSimplify[D[fiveOhOne, b] == f[b]]
If[(Im[a] - Im[b])*
((-Im[b])*Re[a] +
Im[a]*Re[b]) > 0 ||
Im[b]/(Im[a] - Im[b]) >=
0 || Im[a]/(-Im[a] +
Im[b]) >= 0,
(a^(1 + c) - b^(1 + c) +
Log[a] + c*Log[a] -
Log[b] - c*Log[b])/
((a - b)*(1 + c)),
Integrate[
1/(a + (-a + b)*x) +
(a + (-a + b)*x)^c,
{x, 0, 1},
Assumptions ->
!((Im[a] - Im[b])*
((-Im[b])*Re[a] +
Im[a]*Re[b]) > 0 ||
Im[b]/(Im[a] -
Im[b]) >= 0 ||
Im[a]/(-Im[a] +
Im[b]) >= 0)]] +
(-a + b)*
If[(Im[a] - Im[b])*
((-Im[b])*Re[a] +
Im[a]*Re[b]) > 0 ||
Im[b]/(Im[a] -
Im[b]) >= 0 ||
Im[a]/(-Im[a] +
Im[b]) >= 0,
(-(1/b) - c/b -
b^c*(1 + c))/
((a - b)*(1 + c)) +
(a^(1 + c) -
b^(1 + c) + Log[a] +
c*Log[a] - Log[b] -
c*Log[b])/((a - b)^2*
(1 + c)), Integrate[
-(x/(a + (-a + b)*x)^
2) + c*x*
(a + (-a + b)*x)^
(-1 + c), {x, 0, 1},
Assumptions ->
!((Im[a] - Im[b])*
((-Im[b])*Re[a] +
Im[a]*Re[b]) > 0 ||
Im[b]/(Im[a] -
Im[b]) >= 0 ||
Im[a]/(-Im[a] +
Im[b]) >= 0)]] ==
1/b + b^c
Hmm.... maybe that's correct. Maybe not. How can I be sure?
The following shows it isn't correct for ALL a,b, and c, but it leaves
open the question of where the error might be. It may be correct for all
other combinations of b and c, for all I know.
Simplify[D[fiveOhOne, b] ==
f[b] /. {b -> 0.1,
c -> 1}]
% /. a -> 0
1/(-0.1 + 1.*a) == 0
False
Setting b only (for every value I've tried) confirms the result:
Simplify[D[fiveOhOne, b] == f[b] /. {b -> 2.1}]
True
Similarly, setting a only seems to confirm the result, if I choose a>=1 or
a==1/2 (apparently).
FullSimplify[D[fiveOhOne, a] == -f[a] /. {a -> 1}]
True
But it fails for lots of other a values. So there must be b values I
should have tried but didn't. It's a hit-or-miss operation.
The 4.2 integral checks out easily:
fourTwo = -((a^(1 + c) -
b^(1 + c) + Log[a] +
c*Log[a] - Log[b] -
c*Log[b])/(1 + c));
Simplify[D[fourTwo, a] + f[a]]
Simplify[D[fourTwo, b] - f[b]]
0
0
But it can be far more complicated, even with a correct answer. Simplify
can't always identify expressions that resolve to zero.
It makes no sense to say "don't use a computer algebra program for things
you can do by hand", when the other kind can't always be checked for
accuracy. That leaves us doing EVERYTHING by hand.
It makes more sense to use the program for all problems, and use the
program to check them all for accuracy. But good luck doing that; it won't
always be easy. Checking easy problems is good practice for the harder
ones.
Bobby
On Fri, 30 Jan 2004 20:52:18 +0100, Peter Klamser <Klamser at t-online.de>
wrote:
> Hi,
>
> there is no errorfree software, a scientist should know this fact.
>
> Everyone should decide, which software he should use and prove if the
> result
> of the software is errorfree.
>
> I never trust a software, I only trust my knowledge.
>
> Integrate[1/x + x^c, {x, a, b}]
>
> returns with Mathematica 5.0.1
>
> (-a + b)*If[(Im[a] - Im[b])*((-Im[b])*Re[a] + Im[a]*Re[b]) > 0 ||
> Im[b]/(Im[a] - Im[b]) >= 0 ||
> Im[a]/(-Im[a] + Im[b]) >= 0, (a^(1 + c) - b^(1 + c) + Log[a] +
> c*Log[a] - Log[b] - c*Log[b])/
> ((a - b)*(1 + c)), Integrate[1/(a + (-a + b)*x) + (a + (-a + b)*x)^c,
> {x, 0, 1},
> Assumptions -> !((Im[a] - Im[b])*((-Im[b])*Re[a] + Im[a]*Re[b]) > 0
> ||
> Im[b]/(Im[a] - Im[b]) >= 0 || Im[a]/(-Im[a] + Im[b]) >= 0)]]
>
> Version 4.2 returns
>
> -((a^(1 + c) - b^(1 + c) + Log[a] + c*Log[a] - Log[b] -
> c*Log[b])/(1 + c))
>
> But in version 5.0 the wrong result appears.
>
> Therefore it was a bug in version 5.0
>
> Regards
>
> Peter
>
> -----Ursprüngliche Nachricht-----
> Von: Bobby R. Treat [mailto:drbob at bigfoot.com]
> Gesendet: Freitag, 30. Januar 2004 10:17
> An: mathgroup at smc.vnet.net
> Betreff: [mg45958] Re: Nasty bug in Integrate (version 5.0)
>
>
> The only reason we know the answer is wrong is BECAUSE it's trivial by
> hand. It makes no sense to trust Mathematica for problems that are too
> hard
> to
> check.
>
> Bobby
>
> Andrzej Kozlowski <akoz at mimuw.edu.pl> wrote in message
> news:<bvapk6$aat$1 at smc.vnet.net>...
>> Mathematica seems to get lost in the large number of checks that it
>> attempts to perform to find the exact conditions (in the complex
>> plane) for the convergence of this integral. Even if a and b are
>> assumed to be real there is a problem that 0 may lie in between them.
>> Only in the case when you specify that the limits are both positive or
>> both negative you get the correct answers because Mathematica can
>> quickly deal with this problem.
>> If you really want to see how much checking Mathematica attempts to do
>> just evaluate:
>>
>> Trace[Integrate[1/x + x^c, {x, a, b}], TraceInternal->True]
>>
>> and be prepared to wait and obtain a huge output.
>>
>> Anyway, this is clearly a bug. However, it is in general a good idea
>> not to ask a computer program to do what is trivial to do by hand,
>> which in this case means evaluating a definite rather than an
>> indefinite integral. Note that
>>
>>
>> Integrate[1/x + x^c, x]
>>
>> x^(c + 1)/(c + 1) + Log[x]
>>
>> and the answer is obtained instantly. If this is what you really were
>> after (you can substitute in the limits for example by using:
>>
>>
>> Subtract @@ (Integrate[1/x + x^c, x] /. x -> #1 & ) /@ {b, a}
>>
>>
>> -(a^(c + 1)/(c + 1)) - Log[a] + Log[b] + b^(c + 1)/(c + 1)
>>
>>
>> Of course this now makes sense only under certain conditions on a,b and
>> c which you now have to determine yourself.
>>
>>
>> On 28 Jan 2004, at 11:19, Math User wrote:
>>
>> > Hi,
>> >
>> > Sorry if this has been discussed before. Could anyone explain why
>> > Mathematica 5.0 is ignoring the term 1/x in the first and second
>> > answers?
>> >
>> > Thanks!
>> >
>> >
>> > In[1]:= Integrate[1/x + x^c, {x, a, b}]
>> >
>> > Out[1]= (-a^(1 + c) + b^(1 + c))/(1 + c)
>> >
>> > In[2]:= Integrate[1/x + x^c, {x, a, b}, Assumptions -> {Element[a,
>> > Reals], Element[b, Reals]}]
>> >
>> > Out[2]= (-a^(1 + c) + b^(1 + c))/(1 + c)
>> >
>> > In[3]:= Integrate[1/x + x^c, {x, a, b}, Assumptions -> {a > 0, b > 0}]
>> >
>> > Out[3]= (-a^(1 + c) + b^(1 + c))/(1 + c) + Log[b/a]
>> >
>> >
>
>
>