quad_qag and quad_qags fail with imagpart( bessel_i(1, %i*x)) arg
Subject: quad_qag and quad_qags fail with imagpart( bessel_i(1, %i*x)) arg
From: Dieter Kaiser
Date: Sat, 22 Oct 2011 23:45:17 +0200
Am Samstag, den 22.10.2011, 14:20 -0700 schrieb Edwin Woollett:
> For cases where integrate, numer fails, I am using
> either quad_qag or quad_qags for numerical integration
> in one dimension.
>
> I automatically split up any integrand into real and
> imaginary parts using realpart and imagpart, and
> later combine the answers.
>
> This has worked with bessel functions ok up to
> my testing of bessel_i(1,%i*x).
>
> (%i4) imagpart(bessel_i(1,%i)),numer;
> (%o4) 0.44005058574493
>
> (%i6) quad_qag(imagpart(bessel_i(1,%i*x)),x,1,3,3,limit=700);
> (%o6) quad_qag(bessel_i(1,%i*x),x,1,3,3,epsrel = 1.0E-8,epsabs = 0.0,
> limit = 700)
>
> (%i11) quad_qags(imagpart(bessel_i(1,%i*x)),x,1,3,limit=700);
> (%o11) quad_qags(bessel_i(1,%i*x),x,1,3,epsrel = 1.0E-8,epsabs = 0.0,
> limit = 700)
>
> Since imagpart(bessel_i(1,%i*x)) is a real number, why can't
> quad_qag and quad_qags cope?
The reason is, that Maxima knows symbolically that bessel_i(1,%i*x) is a
pure imaginary. Therefore the imagpart automatically simplifies:
(%i1) imagpart(bessel_i(1,%i*x));
(%o1) bessel_i(1, %i x)
But unfortunately because of rounding errors Maxima gets a small
realpart when evaluating the function numerically.
(%i3) imagpart(bessel_i(1,%i));
(%o3) bessel_i(1, %i)
(%i4) %,numer;
(%o4) .4400505857449336 %i + 2.694443716532522e-17
At some places of the code for the Bessel functions the numerical
routines takes care of this small contributions, but not for bessel_i. I
will have a look at this problem.
A workaround is to use an order, which is a floating point number.
Maxima only simplifies for an integer order.
(%i1) imagpart(bessel_i(1.0,%i*x));
(%o1) imagpart(bessel_i(1.0, %i x))
(%i2) quad_qag(imagpart(bessel_i(1.0,%i*x)),x,1,3,3,limit=700);
(%o2) [1.0252496414599, 1.138255757937468e-14, 31, 0]
Dieter Kaiser