rtoy / maxima

A Clone of Maxima's repo
Other
0 stars 0 forks source link

integrate((1-cos(x))^(3/2), x, 0, 2*%pi) wrong by factor of 2, correct in earlier version #2703

Open rtoy opened 3 months ago

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-06 18:11:35 Created by tomasriker on 2024-04-08 07:26:21 Original: https://sourceforge.net/p/maxima/bugs/4282


Current Maxima gives an answer that's half the correct answer. An earlier version of Maxima (5 41.0) gives the correct answer, which is 2^(9/2)/3.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-06 18:11:36 Created by willisbl on 2024-04-08 10:15:01 Original: https://sourceforge.net/p/maxima/bugs/4282/#0900


I think this definite integral gets evaluated by first finding an antiderivative and second applying the fundamental theorem of calculus. The integrand is continuous on [0, 2 %pi], so the antiderivative is differentiable on [0, 2 %pi]. But a graph of the antiderivative suggests that the antiderivative isn't differentiable at %pi. In such cases, the integration code is supposed to notice this and integrate over the required intervals automatically--something about that process is broken, I think.

A bit of evidence: integrating from [0, %pi] and from [%pi,2%pi] gives the correct value

(%i5)   xxx : integrate((1-cos(x))^(3/2), x,0,%pi) +  integrate((1-cos(x))^(3/2), x,%pi,2*%pi);

(%o5)   2^(9/2)/3
(%i6)   yyy  : integrate((1-cos(x))^(3/2), x,0,2*%pi);
(%o6)   2^(7/2)/3

Bugs in the limit code might have a role in this bug, but I think it's more than that. There have been some changes to the simplification code for atan2 and maybe for limits of atan2 expressions as well.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-06 18:11:40 Created by tomasriker on 2024-04-08 11:19:27 Original: https://sourceforge.net/p/maxima/bugs/4282/#b4b6


@willisbl I'm not sure about Maxima using the Fundamental Theorem of Calculus here. The antiderivative has the same value at 0 and 2*%pi. Shouldn't Maxima then give an answer of zero?

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-06 18:11:43 Created by rtoy on 2024-04-08 15:23:53 Original: https://sourceforge.net/p/maxima/bugs/4282/#67ab


Traced through some of the code. intsc1 is called which eventually calls intsc0. This eventually calls antideriv. This returns:

-(((sqrt(2)*cos(3*x)-9*sqrt(2)*cos(2*x)-9*sqrt(2)*cos(x)+sqrt(2))
 *sin((3*atan(tan(x))+3*%pi)/2)
 +(-(sqrt(2)*sin(3*x))+9*sqrt(2)*sin(2*x)+9*sqrt(2)*sin(x))
  *cos((3*atan(tan(x))+3*%pi)/2))
 /6)

Using this result, sin-cos-intsubs is called. This eventually calls intsubs.

This looks for discontinuities in the result and finds none so it does a substitution using same-sheet-subs.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-06 18:11:47 Created by tomasriker on 2024-04-08 16:45:08 Original: https://sourceforge.net/p/maxima/bugs/4282/#808b


Interesting. Plotting that antiderivative between 0 and 2*%pi, it seems like it's just wrong. Also, it is different from the antiderivative that we get by simply doing integrate((1-cos(x))^(3/2), x). Why is that?

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-06 18:11:50 Created by tomasriker on 2024-04-08 17:34:28 Original: https://sourceforge.net/p/maxima/bugs/4282/#531d


I identified the commit that broke this integral, made by @robert_dodier in 2022: [fe1d8cf427aa69331147cc94e8ac993225333434] Surprisingly, it is not directly related to integration, but to auto-loading.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-06 18:11:54 Created by rtoy on 2024-04-08 22:45:42 Original: https://sourceforge.net/p/maxima/bugs/4282/#808b/4c76


I think the main difference is that antideriv sets generate-atan2 to nil. Using integrate has generate-atan2 set to t, the default.

I also have a hard time getting maxima to figure out if the derivative is actually equal to the integrand. But plotting the integrand and the derivative shows that they appear to be the same for 0 to %pi, but the derivative is the negative of the integrand from %pi to 2*%pi.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-06 18:11:57 Created by robert_dodier on 2024-04-09 16:55:11 Original: https://sourceforge.net/p/maxima/bugs/4282/#531d/95c5


Thanks for investigating. I see that with commit fe1d8c~1, the result is (2^(9/2))/3 the first time, but (2^(7/2))/3 for the second and following times.

Commit fe1d8c says something about MEVAL* cleaning stuff up via CLEARSIGN and VARLIST. I see that VARLIST has a bunch of stuff on it after calling integrate in this case -- I don't know if that's specific to this case or if that usually happens.

I tried :lisp (meval* $%i1) where %i1 is the integral in question, but I get the second, incorrect result. Not sure what to make of that.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-06 18:12:01 Created by robert_dodier on 2024-04-17 04:19:51 Original: https://sourceforge.net/p/maxima/bugs/4282/#0892


I looked at this some more and I think that there are some different things going on. (1) integrate((1-cos(x))^(3/2), x, 0, 2*%pi) calls SININT to compute an antiderivative, which you can see via display2d: false and trace(?sinint). But that antiderivative isn't valid over the whole range (0, 2*%pi). (2) Under some conditions that invalid antiderivative gets simplified to another invalid derivative, this time containing atan(tan(x)) instead of atan2(sin(x), cos(x)). That other pseudo-antiderivative has somewhat different properties. In particular I think limit finds a different limit at 2*%pi. (3) The conditions for simplifying the antiderivative aren't clear. The stuff about cleaning up CLEARSIGN and VARLIST might come into play.

I think it's just luck that Maxima ever got the correct result -- there is too much that's wrong with this picture. I think the place to start is to try to get SININT to return a valid antiderivative, or, failing that, to detect that it can't, and indicate a failure (I don't remember if that means to return NIL or to throw something or what). Independently, we can investigate why limit is returning different results for nominally equivalent expressions, and also why the trig expressions get simplified in different ways -- returning different results for the same problem is a bug. We might consider opening separate bug reports for these two or three aspects.

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-06 18:12:04 Created by rtoy on 2024-04-20 14:24:40 Original: https://sourceforge.net/p/maxima/bugs/4282/#0892/ce87


FWIW, the antiderivative comes from the Risch integrator. I traced rischint and see that it's called and returns the antiderivative that we see. The antiderivative appears to be symmetrical about %pi, just like the integrand. However integrate((1-cos(x))^(3/2),x,0,%pi) returns -2^(7/2)/3, which doesn't even have the correct sign since the integrand is positive over that interval.