integrate bessel_k(2,%i*y) surprise



Am Donnerstag, den 27.10.2011, 21:18 +0200 schrieb Dieter Kaiser:
> Am Donnerstag, den 27.10.2011, 10:13 -0700 schrieb Edwin Woollett:
> > On Oct. 26, 2011, Raymond Toy wrote:
> > -------------------------------
> > >> (%o26) 1.736293756153568-2.1367387806763253E+42*%i
> > >> ----------------------------------------------------------------
> > >> Where does the 10^42 come from??
> > >
> > >It comes from the bessel_y(1,100*%i) term that integrate produces.   I
> > >checked the value of this and I think it's correct.  I assume that
> > >integrate produces the incorrect result.  I didn't check that.
> > ----------------------------------------------
> > I agree that bessel_y(1,100*%i) produces 10^42.
> > I also assume that integrate produces the incorrect result.
> > 
> > Wolfram alpha gives a numerical value
> >   NIntegrate[BesselK[2,y*I],{y,1,100}]  --->
> > 
> >     -1.42033 + 1.62942*I
> > 
> > and gives the indefinite integral in terms of the hypergeomentric
> > function: the Meijer G-function.
> 
> I have already corrected the integrals for bessel_i and bessel_y. I will
> check the integrals for bessel_k, too. I suppose an error in the
> implementation of the formula.

I have found an error in the integral of bessel_k. After correcting the
formula I get:

(%i2) integrate(bessel_k(2, %i*y),y,1,100);

(%o2) (2*%i-50*%pi*struve_l(0,100*%i))*bessel_k(1,100*%i)
 +((%pi*struve_l(0,%i)-4*%i)*bessel_k(1,%i)
  +%pi*struve_l(-1,%i)*bessel_k(0,%i))
  /2-50*%pi*struve_l(-1,100*%i)*bessel_k(0,100*%i)

(%i3) rectform(%),numer;
(%o3) 1.629424764063958*%i-1.420329286158196

Wolfram alpha gets the same numerical result.
I will commit the correction soon.

Dieter Kaiser