Series solution to modified Bessel function second kind zeroth order

Im trying to show that the series solution to the bessel coefficient of 2. kind. and zero order gives the same result as the Matlab function besselk(0,x) - But I cant

The series solution looks like following:

K0(x)=-(gam+log(x/2))*I0(x)+sum(from_n=0 to infinity{x^(2*(n))/(2^(2*(n))*factorial(n)^2)}(1/(n+1))

It should approach zero at any x, but I cant make it do so. The error should not be within the I0(x), since i get this to fit perfectly with besseli(0,x).

My code look like following

gam=0.5772156649
x=5;
x2=16;

MATLAB_BESSELi= besseli(0,x);
MATLAB_BESSELi2= besseli(0,x2);
MATLAB_BESSELk= besselk(0,x);

for n=0:10
    I(n+1)=x^(2*n)/(2^(2*n)*factorial(n)^2) % Equation 2.3.3 in MYERS
    sumI= sum(I);

      I2(n+1)=x2^(2*n)/(2^(2*n)*factorial(n)^2) % Equation 2.3.3 in MYERS
      sumI2= sum(I2);

%     figure(1)
%     subplot(1,2,1)
%     plot(n,sumI,'xb','markersize',10); hold on
%     plot(n,MATLAB_BESSELi,'or'); hold on
%     legend('Analytical expression','Matlab values','location','se')
end

for n=0:10;

      LED1=-(gam+log(x/2))*sumI
      LED2(n+1)=x^(2*(n))/(2^(2*(n))*factorial(n)^2)*(1/(n+1))
      sumLED2=sum(LED2)

      K=LED1+sumLED2

          figure(2)
          plot(n, K,'xb'); hold on
          plot(n, besselk(0,x),'xr')

end

Does anyone know how to solve this?

Please excuse my bad programming behaviors when looking at my script.

Thanks

[Star Strider]

I can’t match what you wrote with my reference (Abramowitz and Stegun, Dover 1972) P. 360, 9.1.13. At the very least, you seem to be missing a factor of (2/pi). I can’t tell if that’s the only problem.

[Simon Sand]

If I look in that reference, my equations should match equation 9.6.12 and 9.6.13. As far as i understand, they should be applicable for Bessel first and second Kind of zeroth order.

Thanks

Answer 1

[Star Strider]

I decided to code that on my own, in part so I could understand the process. (I haven’t had to code my own series evaluations since I got MATLAB a couple decades ago.)

Here’s my solution (Abramowitz and Stegun, 9.6.12 and 9.6.13):

gam=0.5772156649;

x = 5;

I0 = 1;
for k1 = 1:10
    I0 = I0 + (((x.^2)./4).^k1)/(factorial(k1)^2);
end
I0m = besseli(0,x);
I0cmp = [I0m  I0]

K0 = -(log(x/2)+gam)*I0;
frx = 0;
for k1 = 1:10
    frx = frx + 1/k1;
    K0 = K0 + (frx*(((x.^2)./4).^k1)/(factorial(k1)^2));
end
K0m = besselk(0,x);
K0cmp = [K0m  K0]

NOTE: My notation corresponds to that in Abramowitz and Stegun, except for using x where they use z.

I tend to use parentheses more than others might, because compilers never manage to guess what I’m thinking.