#### 7.58 Bug in Elliptic Integral, Maple V to Maple 8 (25.7.01)

##### 7.58.1 Tim McLarnan

I’m not a big expert on special functions, but it appears to me that Maple has been incorrectly calculating the values of incomplete elliptic integrals at least since Maple V. The problem I’ll describe here exists on my Intel box under both Maple 6 and Maple 7. The values computed by Maple V release 5 on a Mac are slightly diﬀerent from Maple 7’s answers, but also wrong.

The errors are not small - in one example, Maple computes a value of 1.80+0.62i when the correct answer is 0.27+0.97i.

Here’s an example with the EllipticE function. Maple reports the value

> evalf(EllipticE(1/5*sqrt(70+35*I),2/5*sqrt(5)-1/5*I*sqrt(5)));

1.801868551 + .6195716925 I



On the other hand, the EllipticE function is deﬁned by an integral which Maple can evaluate explicitly,

> z := evalf(1/5*sqrt(70+35*I)):
k := evalf(2/5*sqrt(5)-1/5*I*sqrt(5)):
evalf(Int(sqrt(1-k^2*t^2)/sqrt(1-t^2),t=0..z));

.2742693842 + .9709289578 I



(This second calculation takes Maple close to 10 minutes on my computer; MuPAD conﬁrms the result in a couple of seconds.)

I think the results of these calculations should be equal, and they are not.

Similarly, with the EllipticF function, Maple gives the value

> evalf(EllipticF(1/5*sqrt(70+35*I),2/5*sqrt(5)-1/5*I*sqrt(5)));

-.008504510023 + 1.374884128 I



Working out numerically the deﬁning integral for this function yields (again conﬁrmed by MuPAD)

> z := evalf(1/5*sqrt(70+35*I)):
k := evalf(2/5*sqrt(5)-1/5*I*sqrt(5)):
evalf(Int(1/sqrt(1-t^2)/sqrt(1-k^2*t^2),t=0..z));

2.193595226 + .5178375879 I



Again, I think these two calculations should have given the same result, and they do not.

Now, it would be easy to dismiss these observations as

(1) Of no consequence, because when was the last time anybody used an incomplete elliptic integral, and

(2) Probably a misunderstanding on the part of someone who doesn’t know that much about special functions. (Guilty.)

Let me argue that instead, this is a problem that could bite any user, and that my values for these functions are right (and Maple’s wrong).

Consider the function sqrt((x^2+1)/(x+2)). This is a positive real function for $$x > -2$$. A plot suggests that int(sqrt((x^2+1)/(x+2)), x=-1..5) should be about 8.

Maple computes this integral using incomplete elliptic integrals, including the two discussed above.

> int(sqrt((x^2+1)/(x+2)), x=-1..5);

1/2
2 182                  1/2            1/2             1/2             1/2
-------- - 2/6825 I 182    (70 + 35 I)    (-45 + 35 I)    (-45 - 35 I)
3

1/2     1/2
(70 + 35 I)     2 5             1/2              1/2
EllipticF(--------------, ------ - 1/5 I 5   ) + 4/6825 182
5           5

1/2             1/2             1/2
(70 + 35 I)    (-45 + 35 I)    (-45 - 35 I)

1/2     1/2                    1/2
(70 + 35 I)     2 5             1/2    2 2              1/2
EllipticE(--------------, ------ - 1/5 I 5   ) - ------ + 2/75 I 2
5           5                      3

1/2           1/2           1/2
(10 + 5 I)    (15 + 5 I)    (15 - 5 I)

1/2     1/2
(10 + 5 I)     2 5             1/2          1/2           1/2
EllipticF(-------------, ------ - 1/5 I 5   ) - 4/75 2    (10 + 5 I)
5          5

1/2     1/2
1/2           1/2           (10 + 5 I)     2 5             1/2
(15 + 5 I)    (15 - 5 I)    EllipticE(-------------, ------ - 1/5 I 5   )
5          5



Maple then asserts that the numerical value of the integral is

> evalf(%);

14.18222460 + 4.701625434 I



an obviously absurd answer. A correct answer can be obtained by replacing a limit of integration by a ﬂoat:

> int(sqrt((x^2+1)/(x+2)), x=-1.0..5);

7.277500982



The ﬁrst explanation that would occur to one here is that Maple’s integral is wrong, but I think this is not so. One can diﬀerentiate Maple’s value for the integral int(sqrt((t^2+1)/(t+2)), t=-1..x) by hand, and one seems to get sqrt((x^2+1)/(x+2)) back, just as one ought. The real clincher, though, is that if one takes Maple’s expression for the deﬁnite integral and plugs in the values I computed above for the EllipticE and EllipticF functions, the result is 7.277500984, which agrees exactly with the result of computing the deﬁnite integral numerically. (There are 4 elliptic integrals in Maple’s deﬁnite integral; Maple seems to compute 2 of them right, and 2 of them wrong.)

A Maple worksheet showing a bit more detail on all this is at

I’ve checked the calculations here with Maple 6 under Windows and Linux, Maple 7 under Linux, and Maple V release 5 under MacOS. The exact wrong answers diﬀer slightly between Maple V and Maple 7, but all these versions show the same problem.

I hope I’m just doing something wrong here, but if not, it’s extremely alarming to me that one seems not to be able to rely on Maple correctly to compute the values of fundamental functions, and that bugs of this sort can persist for so many years.

##### 7.58.2 Frederic PASCAL (14.9.01)

for computing the numerical value of the integral

Int(sqrt((x^2+1)/(x+2)), x=-1..5)

you ﬁrst compute an algebraic expression of the integral then you ask for a numerical evaluation of the later expression.

You are increasing all the chances to get a wrong result and it is the case. Indeed it is a problem of rounding errors.

Let’s take a more simple integral :

Int(x^n/(x+a),x=0..1) with n=10 and a=10

The integral should be less than 1.

The algebraic expression is the following :

>a:=10: n:=10: int(x^n/(x+a),x=0..1);

10000000000*ln(11)-10000000000*ln(10)-300227066381/315



Now when Maple (MapleV5) computes

> a:=10; n:=10; evalf(int(x^n/(x+a),x=0..1));

2.0



Maple ﬁrst computes the algebraic expression, then evaluates the ﬂoating expression with Digits:=10. During this step, it is doing a rounding error since you are computing a small value as the diﬀerence of large numbers.

>evalf(10000000000*ln(11));evalf(10000000000*ln(10));evalf(300227066381/315);
11
.2397895273 10

11
.2302585093 10

9
.9531017980 10



The computation of ln is correct (well the 10 ﬁrst digits, the next ones 0 are not) but you don’t have enough digits to compute an exact diﬀerence.

Here the exact value is (now Maple doesn’t compute an algebraic expression)

> evalf(Int(x^n/(x+a),x=0..1));

.008327965519



Remarks with Digits:=25 you get an almost correct value

> Digits:=25:a:=10:n:=10:evalf(int(x^n/(x+a),x=0..1));

.0083279655188951



In your case, the problem is of the same order (i also don’t know much about special functions) :

> r:=int(sqrt((x^2+1)/(x+2)), x=-1..5):
> Digits:=10: evalf(r);
> Digits:=25: evalf(r);
evalf(Int(sqrt((x^2+1)/(x+2)), x=-1..5));
-8
7.277500973 + .3 10   I

-22
7.277500981950899644491168 - .19 10    I

7.277500981950899644491212



I could not reproduce the ﬁrst two results (U. Klein)

Using numerics with Maple must be done with precaution.

##### 7.58.3 James R. FitzSimons (24.9.01)

Maple has a bug.

Maple is giving the wrong answer using elliptic integrals.

> evalf(Int(sqrt((x^2+1)/(x+2)), x=-1..5));

7.277500982

> evalf(int(sqrt((x^2+1)/(x+2)), x=-1..5));

14.18222461 + 4.701625434 I



This is the numeric result using DERIVE.

7.2775009819508996444