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: Edwin Woollett
Date: Sat, 22 Oct 2011 15:14:13 -0700
On Oct. 22, Dieter Kaiser wrote:
----------------------
>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]
------------------------------------------
Thanks for the info.
Your method of using 2.0 instead of 2 for the order of
the bessel function works for me on a "one off" basis:
(%i34) quad_qags(imagpart(bessel_i(1.0,%i*x)),x,1,3,limit=700);
(%o34) [1.0252496414599,1.1382557579374678E-14,21,0]
(%i35) quad_qag(imagpart(bessel_i(1.0,%i*x)),x,1,3,3,limit=700);
(%o35) [1.0252496414599,1.1382557579374679E-14,31,0]
but for an automated version of nintegrate, it would be painful
to have to replace the integers by floats in the args of an arbitrarily
complicated expression.
A better workaround for my purposes is to use my fchop function
from the ch. 9 software:
quad_qag(fchop(imagpart(bessel_i(1.0,%i*x))),x,1,3,3,limit=700);
(%o36) [1.0252496414599,1.1382557579374679E-14,31,0]
(%i37) quad_qags(fchop(imagpart(bessel_i(1.0,%i*x))),x,1,3,limit=700);
(%o37) [1.0252496414599,1.1382557579374678E-14,21,0]
Ted