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 different 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 defined 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 confirms 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 defining integral for this function yields (again confirmed 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 float:
> int(sqrt((x^2+1)/(x+2)), x=-1.0..5); 7.277500982
The first explanation that would occur to one here is that Maple’s integral
is wrong, but I think this is not so. One can differentiate 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 definite 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 definite integral numerically. (There are 4 elliptic integrals in
Maple’s definite 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 differ 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 first 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 first computes the algebraic expression, then evaluates the floating expression with Digits:=10. During this step, it is doing a rounding error since you are computing a small value as the difference 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 first digits, the next ones 0 are not) but you don’t have enough digits to compute an exact difference.
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 first 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