This forum question led to a discussion of a bitwise magazine review that compared Mathematica 5.2 and Maple 10. In that review the author struggled to get the following numeric integral to compute accurately and quickly in Maple.

evalf(Int(BesselJ(0, 50001*x)*x*exp(I*(355*x^2*1/2)), x = .35 .. 1));

Below, I reproduce an attempt at computing an accurate result quickly in Maple. I'm copying it here because that thread got quite long and messy.

I'll use three things. A fast hardware Bessel J0. The NAG routine d01akc specialized for oscillatory integrands, but with the max_num_subints optional parameter exposed for use. And a comment by Doug Meade about splitting the integrand into real and complex parts.

Some readers commented on an earlier version of this blog entry which contained a source fragment containing a similar define_external call involving a function in the GNU Scientific Library (which is covered by the GNU General Public License). If such use truly implies that that code fragment is also covered by the license then it should be mentioned that the GSL Library comes with no warranty and may not be used in proprietary software. However this link may indicate that such a code fragment may not be allowed. I've exchanged the code fragment that would set up a call to the GSL double precision Bessel J0 function with one which would call the equivalent function from a standard Unix math runtime library such as libm.so.

> kernelopts(printbytes=false):
>
> # The OS's runtime math library libm.so may have Bessel J0.
> BesselJ0 := define_external('j0',
> 'r'::float[8],'RETURN'::float[8] ,
> LIB="/usr/lib64/libm.so"):
>
# a quick sanity check.
> evalf[trunc(evalhf(Digits))](BesselJ(0,50000.0));
                             -0.00256784217783324
 
> BesselJ0(50000.0);
                            -0.00256784217783323976

>
#BesselJ0(50001*x)*x*exp(I*(355*x^2*1/2));
> fR := proc(x)
> global oparam;
> option hfloat;
>   BesselJ0(oparam*x)*x*cos(355/2*x^2);
> end proc:
>
> fI := proc(x)
> global oparam;
> option hfloat;
>   BesselJ0(oparam*x)*x*sin(355/2*x^2);
> end proc:
>
> NAGd01akcM := ExternalCalling:-DefineExternal(
>     'hw_d01akc',
>      ExternalCalling:-ExternalLibraryName("int",'HWFloat')):
>                                                                                 
> d01akc:=proc(f,Left,Right,epsabs,epsrel,max_num_subint)
> local res;
>   res:=NAGd01akcM(eval(f),Left,Right,epsabs,epsrel,max_num_subint,
>              'abserr','num_subint','fun_count',0);
>   res, abserr, num_subint, fun_count;
> end proc:
>
> Digits:=trunc(evalhf(Digits)):
>
> oparam:=50001:
> gc():(st,ba,bu):=time(),kernelopts(bytesalloc),kernelopts(bytesused):
> (solR,abserr,num_subint,fun_count):=
>     d01akc(fR,.35,1.,Float(5,-12),Float(5,-12),5000):
> (solI,abserr,num_subint,fun_count):=
>     d01akc(fI,.35,1.,Float(5,-12),Float(5,-12),5000):
> time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;
                           3.096, 36169248, 68621344
 
> solR+I*solI;
                                   -7                          -7
            0.318175350197895498 10   - 0.779830112408473511 10   I
 
>
# With the value of oparam for which bitwise article shows a Mma result
# of -2.00775e-7 + I*4.27539e-8 .
> oparam:=10000:
> (solR,abserr,num_subint,fun_count):=
> d01akc(fR,.35,1.,Float(5,-12),Float(5,-12),5000):
> (solI,abserr,num_subint,fun_count):=
> d01akc(fI,.35,1.,Float(5,-12),Float(5,-12),5000):
> solR+I*solI;
                                   -6                          -6
           -0.200775234034556560 10   + 0.427538846202058781 10   I
 
>
> oparam:=300000:
> gc(): st:=time():
> (solR,abserr,num_subint,fun_count):=
> d01akc(fR,.35,1.,Float(5,-12),Float(5,-12),5000):
> time()-st, solR, abserr;
                                        -8                         -11
         34.060, 0.219400784174614172 10  , 0.472640230574086231 10

Perhaps the really high precision needed by Maple to compute the integral at high values of the parameter is only due to evalf(Int()) trying to use its _Dexp or _Gquad methods. Perhaps the results above are accurate, when using hardware precision J0 alongside a quadrature routine specialized for oscillatory integrands. Next question: how to verify them all.

I was pretty happy with 3sec to compute a result for oparam=50001. Checking the fun_count values shows that there were about 125000 calls to BesselJ0 in those 3sec. I think that it would be even less if Maple's evalhf allowed call_external. If that were true then, since procedure fR has option hfloat, the results from BesselJ0 might never need to get converted from hardware to software and back to hardware. The evidence that they do get converted is that quite a bit of garbage collection takes place and bytesused goes up.

acer