7.83 bug in integration 3, Maple 5 to Maple 6 (7.12.98)

7.83.1 Douglas B. Meade

NOTE: This is, essentially, a bug report. I welcome help proposing a complete patch but expect that Maple’s Technical Support will take responsibility for seeing that similar mistakes have not been made elsewhere in the integration system and for seeing that my suggestions do not cause additional problems.

I have discovered a bug deep within Maple’s integration procedures. This problem was discovered during a conversation with Jon Johnson at last month’s ICTCM.

Jon’s question was: why does Maple (R5) fail to evaluate the following iterated integral

> f := sqrt( (cos(x)-cos(y))^2 + (sin(x)-sin(y))^2 ); 
> A := Int( Int( f, y=0..2*Pi ), x=0..2*Pi ); 
> value ( A );
 

The error message that is produced by the last command is:

Error, (in int/ellalg/trxstandard/4) int/ellalg/trxstandard/4 uses a 
7th argument, L, which is missing
 

A little exploration (with tracelast, showstat, and debug) led me to discover that `int/ellalg/trxstandard` calls `int/ellalg/trxstandard/4` with only 6 arguments.

Arguments 7 and 8 (L and U) are required.

> showstat( `int/ellalg/trxstandard` );
 

The ninth argument is used at the end of the procedure

> showstat( `int/ellalg/trxstandard/4`, 1..20 );
 

My conjecture is that the default values for L and U are -infinity and infinity, respectively. I do not know what \(z\) represents so cannot determine a suitable default value. (I tried 0, but this was not successful.)

Here is the code that I used to

> `int/ellalg/trxstandard/4orig` := op( `int/ellalg/trxstandard/4` ); 
 
> `int/ellalg/trxstandard/4` := proc(lc, a, b, c, d, x, L, U, z) 
>     local LL, UU, zz; 
>     if nargs<7 then LL:=-infinity fi; 
>     if nargs<8 then UU:= infinity fi; 
> #    if nargs<9 then zz:= 0        fi; 
>     `int/ellalg/trxstandard/4`(lc, a, b, c, d, x, LL, UU, z); 
> end;
 

With this "fix" the evaluation of the integrals progresses a little further but halts with a message similar to the original message.

If someone knows what value to assign to \(z\), I’ll be glad to continue my testing. Otherwise, I hope Maple’s Technical Support will be able to come up with a comprehensive solution to this problem.

With Maple 6 I got no error message, but also no result. It is corrected with Maple 7 (U. Klein)

7.83.2 Joe Riel (15.12.98)

Note that the call to `int/ellalg/trxstandard/4` in `int/ellalg/trxstandard` is

  `int/ellalg/trxstandard/4`(lc, op(srL), x, L, U, z)
 

This suggests that the problem is in op(srL), it should expand to 4 elements. Looking at the assignment to srL, it appears that the culprit is the procedure `int/ellalg/trxstandard/roots`.

An obvious error in that procedure is the assignment to r after the if then else statement. It is

 r := traperror(sort(f, `int/ellalg/trxstandard/compare`));
 

This is clearly wrong because f is an unassigned local variable. Sorting it makes no sense. Either a previous assignment to f is missing or this line should probably be

 r := traperror(sort(rts, `int/ellalg/trxstandard/compare`));
 

Making this change returned an answer for the inner integral, 0, which is obviously wrong (and the same wrong answer that R4 returns).

So I cannot give you a solution, but I’m confident that the problem lies in `int/ellalg/trxstandard/roots`.

If you restrict \(x\) to the open range \(0\) to \(\pi \), i.e.

assume(x, RealRange(Open(0),Open(Pi)):
 

then the `.../compare` routine used in `.../roots` is able to sort the solutions. However, because the expression has a double root at sin(x)/(1+cos(x)), only three values are in the list. The `.../trxstandard` procedure generates an error because it attempts to access a [nonexistent] fourth element in the list returned by `.../roots`. This is with the substitution of rts for f in `.../roots`.

When `.../roots` [and possibly `.../trxstandard] is eventually correctly patched, I suspect that if you are ever to evaluate the double integral you will have to split the outer integral into two pieces, 0..Pi and Pi..2*Pi, and integrate these separately, using assumptions on \(x\).

7.83.3 Graham P. McCauley (18.1.99)

On 07 Dec 1998 Douglas Meade submitted a report on what he thought was a bug in the R5 integration procedures. I have tried to follow this up, and suspect that the bug may really be in the routine `solve`.

Douglas found that

> f:=sqrt((cos(x)-cos(y))^2+(sin(x)-sin(y))^2); 
> A:=Int(Int(f,x=0..2*Pi),y=0..2*Pi); 
> value(A);
 

gave rise to the message

Error, (in int/ellalg/trxstandard/4) int/ellalg/trxstandard/4 uses a 7th 
argument, L, which is missing
 

He used `tracelast` to locate the offending call

#(`int/ellalg/trxstandard`,26): `int/ellalg/trxstandard/4`(lc,op(srL),x,L,U,z)

and tried to add extra arguments at the end of the list.

It seems to me that `op(srL)` should provide four values (the roots of a quartic polynomial in ‘x‘ whose coefficicents are trigonometric expressions in ‘y‘). Use of the debugger shows that ‘srL‘ in the present case is an unassigned variable ‘f‘.

This variable is in fact local to `int/ellalg/trxstandard/roots`. In the present case it remains unassigned in a nest of conditional statements. These appear to be an attempt to deal with the fact that ‘solve‘ might fail to report some roots of a quartic ‘sr‘ in ‘x‘ whose coefficents are trigonometric expressions in ‘y‘. Adding a final clause

else f := [op(rts), simplify(eval(tcoeff(sr,x)/lc/convert(rts,`*`)))]
 

to the nest ensures that if only one root is missing (and none of them are zero!) it is added to the list ‘rts‘. [I have checked that it is indeed a slightly disguised repeat of ‘rts[3]‘ ]. The new list is assigned to ‘f‘ which is sorted in some way and RETURNed to become ‘srL‘ above. Now MapleV soldiers on happily to produce the value zero.

Two mysteries remain:

(1) why does ‘solve‘ miss the repeated root?

(2) why does the call of `int/ellalg/trxstandard` have a final argument ‘z‘ that is given the value ‘droot‘ that seems to suggest a double root was already suspected? [I haven’t found where this argument is used. It is the one that was troubling Douglas Meade in his proposed "fix".] A final word on the final result. The value zero is presumably not what the originator of the integral expected. We can see how it comes about if we tidy up the integrand and divide the square region of integration diagonally by writing

> g:=combine(simplify(f,trig),trig); 
 
          g := (2-2*cos(x-y))^(1/2) 
 
> A1:=Int(Int(g,x=0..y),y=0..2*Pi): 
> A2:=Int(Int(g,x=y..2*Pi),y=0..2*Pi):
 

Both A1 and A2 give the value 4*Pi*4^(1/2). [This appears to come out without invoking `int/ellalg/trxstandard`.] One naturally reads these values as 8*Pi since ‘sqrt‘ of a real positive quantity is taken to be positive. However when we ask MapleV to integrate over the whole square at once, it effectively insists on an analytic interpretation of the integrand that makes the two triangles cancel!