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

http://www.cs.earlham.edu/~timm/ellipticIntBug.mws

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.

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.

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`