7.146 buggy integration in Maple 6 and Maple 7 (8.2.01)

7.146.1 Rouben Rostamian
7.146.2 Helmut Kahovec (10.2.01)
7.146.3 Preben Alsholm (12.2.01)
7.146.4 Robert Israel (12.2.01)

7.146.1 Rouben Rostamian

This either shows an integration bug in Maple 6 or it shows that I am abusing it:

> expr := exp(-(x-1)^2) / sqrt(x); 
                                                2 
                                    exp(-(x - 1) ) 
                            expr := -------------- 
                                          1/2 
                                         x 
 
> evalf(int(expr,x=0..+infinity)); 
                                 .4118623312 I 
 
> evalf(Int(expr,x=0..+infinity)); 
                                  1.973732150
 

It is corrected with Maple 8. (U. Klein)

7.146.2 Helmut Kahovec (10.2.01)

It’s a bug. A couple of days ago the following integration bug was discussed in the German newsgroup de.sci.mathematik:

> restart; 
> expr:=1/(1+(4*x-2)^2);
 

Maple6 returns the wrong antiderivative:

> int(expr,x); 
 
                         1/2 arctan(4 x - 2) 
 
> simplify(expr-diff(%,x)); 
 
                                   1 
                          - ---------------- 
                                    2 
                            5 + 16 x  - 16 x
 

The correct antiderivative is:

> 1/4*arctan(4*x-2); 
 
> simplify(expr-diff(%,x)); 
 
                                  0
 

Thus we can produce an analogous example to yours:

> evalf(int(expr,x=0..1)); 
 
                             1.107148718 
 
> evalf(Int(expr,x=0..1)); 
 
                             .5535743589
 

This is corrected with Maple 7. (U. Klein)

7.146.3 Preben Alsholm (12.2.01)

This is clearly a bug. It also disturbs me that doing your computations several times (always starting with a restart) I occasionally get the answer .4118623312 WITHOUT the I.

7.146.4 Robert Israel (12.2.01)

Certainly it’s a bug. In this case it occurs in `int/definite/meijerg`. It seems to be basically a branch-cut problem.

> int(exp(-(x-1)^2) / sqrt(x), x=0..infinity); 
 
1/2 I exp(-1) (-3 sqrt(Pi) sqrt(2) exp(1/2) BesselK(3/4, 1/2) 
 
         + sqrt(Pi) sqrt(2) exp(1/2) BesselK(7/4, 1/2))/sqrt(Pi)
 

Or more generally:

> int(exp(-(x-t)^2) / sqrt(x), x=0..infinity); 
 
                             / 
                 2           | 
       1/2 exp(-t ) sqrt(-t) | 
                             | 
                             \ 
 
                                     2                    2 
           sqrt(Pi) sqrt(2) exp(1/2 t ) BesselK(3/4, 1/2 t ) 
        -3 ------------------------------------------------- 
                                   2 
                                  t 
 
                                                            \ 
                                     2                    2 | 
         + sqrt(Pi) sqrt(2) exp(1/2 t ) BesselK(7/4, 1/2 t )|/sqrt(Pi) 
                                                            | 
                                                            /
 

For \(t < 0\) this seems to be OK. Now BesselK(3/4, z) and BesselK(7/4, z) both have branch points at \(z=0\), and the usual branch cuts are on the negative real axis, which would correspond to the imaginary axis for BesselK(3/4, 1/2 t^2) and BesselK(7/4, 1/2 t^2). And Maple 6’s answer above does have branch cuts on the imaginary axis. Of course the integral itself is an entire function of t.

Maple V Release 5.1 had:

> S:= int(exp(-(x-t)^2) / sqrt(x), x=0..infinity); 
 
                          /  3/2  1/2                          2 
               2      1/2 |Pi    2    hypergeom([1/4], [1/2], t ) 
S := 1/2 exp(-t ) (-t)    |-------------------------------------- 
                          |                         1/2 
                          \          GAMMA(3/4) (-t) 
 
                                                             \ 
           1/2                1/2                          2 |   /   1/2 
     - 2 Pi    GAMMA(3/4) (-t)    hypergeom([3/4], [3/2], t )|  /  Pi 
                                                             | / 
                                                             / 
> S1:= convert(S, StandardFunctions); 
 
                          /  3/2   2 1/4          2                     2 
               2      1/2 |Pi    (t )    exp(1/2 t ) BesselI(-1/4, 1/2 t ) 
S1:= 1/2 exp(-t ) (-t)    |----------------------------------------------- 
                          |                        1/2 
                          \                    (-t) 
 
         3/2     1/2          2                    2 \ 
       Pi    (-t)    exp(1/2 t ) BesselI(1/4, 1/2 t )|   /   1/2 
     - ----------------------------------------------|  /  Pi 
                            2 1/4                    | / 
                          (t )                       / 
 
> evalf(subs(t=1,S1)); 
 
                                  1.973732150
 

I believe that this version is correct for all t.