Am Dienstag, den 25.10.2011, 10:58 -0700 schrieb Edwin Woollett:
> On Oct. 24, 2011, Dieter Kaiser wrote in a bug fix message:
> -----------------------------------------
> >Bugs item #3427581, was opened at 2011-10-23 19:39
> >Message generated for change (Comment added) made by crategus
> >[snip]
> >Message:
> >This bug has been fixed in bessel.lisp revision 23.10.2011.
> >
> >(%i1) imagpart(bessel_i(1,%i*x));
> >(%o1) -%i*bessel_i(1,%i*x)
> >
> > Nevertheless, the reported example is not solved:
> >
> >(%i2) quad_qag(imagpart(bessel_i(1,%i*x)),x,1,3,3,limit=700);
> >(%o2) quad_qag(-%i*bessel_i(1,%i*x),x,1,3,3,epsrel = 1.e-8,epsabs = 0.0,
> > limit = 700)
> >
> >The reason is a small imaginary part because of rounding errors:
> >
> >(%i3) imagpart(bessel_i(1,%i));
> >(%o3) -%i*bessel_i(1,%i)
> >(%i4) rectform(%),numer;
> >(%o4) .4400505857449336-2.694443716532522e-17*%i
> >
> >For other Bessel functions Maxima knows pure real or imaginary results and
> >cuts off the rounding errors accordingly. This can be added for the
> >bessel_i function too.
> ----------------------------------------------------------
> nint.mac is now using the new bessel.lisp.
>
> Until the rounding error problem of bessel_i is fixed, we can use
> a workaround suggested by you (Dieter Kaiser), namely to use
> a float number instead of an integer for the order of the bessel
> function.
> -------------------
> (%i1) load(nint);
> (%o1) "c:/work2/nint.mac"
>
> (%i2) quad_qag(imagpart(bessel_i(1,%i*x)),x,1,3,3,limit=700);
> (%o2) quad_qag(-%i*bessel_i(1,%i*x),x,1,3,3,epsrel = 1.0E-8,epsabs = 0.0,
> limit = 700)
>
> (%i3) quad_qag(imagpart(bessel_i(1.0,%i*x)),x,1,3,3,limit=700);
> (%o3) [1.0252496414599,1.1382557579374679E-14,31,0]
>
> (%i4) quad_qag(imagpart(bessel_i(1.0,%i*x)),x,1,10^4,3,limit=700);
> (%o4) [0.77229384691138,1.1418381427522419E-12,31713,0]
>
> (%i5) quad_qag(imagpart(bessel_i(2.0,%i*x)),x,1,3,3,limit=700);
> (%o5) [8.2026594591804487E-17,9.1067813947414708E-31,31,0]
>
> (%i6) quad_qag(imagpart(bessel_i(2.0,%i*x)),x,1,10^4,3,limit=700);
> (%o6) [-2.2013154155922133E-16,5.7377810005552024E-15,43369,1]
>
> (%i7) quad_qag(imagpart(bessel_i(3.0,%i*x)),x,1,3,3,limit=700);
> (%o7) [-0.28287409015192,3.1405332795655258E-15,31,0]
>
> (%i8) quad_qag(imagpart(bessel_i(3.0,%i*x)),x,1,10^4,3,limit=700);
> (%o8) [-0.98790703708805,1.0431047896924372E-12,31713,0]
> ------------------------------
Thank you for the message. I will have a second look at the
implementation of besse_i to remove the rounding errors for floating
point evaluation in the case of pure real or imaginary results.
Dieter Kaiser