Bessel Function Calculation in Matlab

I was trying to get the result for Bessel function of first kind by using the bessel function main equation and Matlab builtin bessel function.But I am getting good result for builtin function but while using the actual equationa at certain point my series starts divergeing instead of converging.I am giving the codes I haved used

n=1
JnE=0;
for z=11:60
    for k=0:100
        JnE = JnE+(((-1)^k))/(factorial(k)*factorial(k+n+1))*(z/2)^(2*k+n);
    end
    fprintf('JnEz=%f\n',JnE)
end

Can anyone tell me where I am doing it wrong?

I am also giving my code by using besselj function

format short
n=1
for z=11:60
    JnE=besselj(n,z)
end

Answer 1

[Walter Roberson]

Are you sure that you want to be totaling from z = 11 to z = 60? You only set JnE to zero before the z loop.

You are running into problems with factorial(k) losing precision.

But your big problem was that you had factorial(k+n+1) when you needed gamma(k+n+1) -- when you rewrote as factorial you forgot to subtract the 1.

syms k

n=1

n = 1

JnE=0;
for z=11:60
    part = symsum((((-1)^k))/(factorial(k)*gamma(k+n+1))*(z/2)^(2*k+n), k, 0, 100);
    JnE = JnE + double(part);
    bess = besselj(n,z);
    fprintf('(%d) part = %f, cumulative = %f, besselj = %f\n', z, part, JnE, bess);
end

(11) part = -0.176785, cumulative = -0.176785, besselj = -0.176785
(12) part = -0.223447, cumulative = -0.400232, besselj = -0.223447
(13) part = -0.070318, cumulative = -0.470550, besselj = -0.070318
(14) part = 0.133375, cumulative = -0.337175, besselj = 0.133375
(15) part = 0.205104, cumulative = -0.132071, besselj = 0.205104
(16) part = 0.090397, cumulative = -0.041674, besselj = 0.090397
(17) part = -0.097668, cumulative = -0.139343, besselj = -0.097668
(18) part = -0.187995, cumulative = -0.327337, besselj = -0.187995
(19) part = -0.105701, cumulative = -0.433039, besselj = -0.105701
(20) part = 0.066833, cumulative = -0.366206, besselj = 0.066833
(21) part = 0.171120, cumulative = -0.195085, besselj = 0.171120
(22) part = 0.117178, cumulative = -0.077908, besselj = 0.117178
(23) part = -0.039519, cumulative = -0.117427, besselj = -0.039519
(24) part = -0.154038, cumulative = -0.271465, besselj = -0.154038
(25) part = -0.125350, cumulative = -0.396815, besselj = -0.125350
(26) part = 0.015046, cumulative = -0.381770, besselj = 0.015046
(27) part = 0.136585, cumulative = -0.245185, besselj = 0.136585
(28) part = 0.130551, cumulative = -0.114633, besselj = 0.130551
(29) part = 0.006934, cumulative = -0.107699, besselj = 0.006934
(30) part = -0.118751, cumulative = -0.226450, besselj = -0.118751
(31) part = -0.133024, cumulative = -0.359475, besselj = -0.133024
(32) part = -0.026589, cumulative = -0.386064, besselj = -0.026589
(33) part = 0.100620, cumulative = -0.285444, besselj = 0.100620
(34) part = 0.132971, cumulative = -0.152473, besselj = 0.132971
(35) part = 0.043991, cumulative = -0.108482, besselj = 0.043991
(36) part = -0.082330, cumulative = -0.190812, besselj = -0.082330
(37) part = -0.130580, cumulative = -0.321392, besselj = -0.130580
(38) part = -0.059162, cumulative = -0.380554, besselj = -0.059162
(39) part = 0.064056, cumulative = -0.316497, besselj = 0.064056
(40) part = 0.126038, cumulative = -0.190459, besselj = 0.126038
(41) part = 0.072101, cumulative = -0.118358, besselj = 0.072101
(42) part = -0.045994, cumulative = -0.164352, besselj = -0.045994
(43) part = -0.119540, cumulative = -0.283892, besselj = -0.119540
(44) part = -0.082803, cumulative = -0.366695, besselj = -0.082803
(45) part = 0.028349, cumulative = -0.338347, besselj = 0.028349
(46) part = 0.111291, cumulative = -0.227056, besselj = 0.111291
(47) part = 0.091269, cumulative = -0.135787, besselj = 0.091269
(48) part = -0.011329, cumulative = -0.147116, besselj = -0.011329
(49) part = -0.101506, cumulative = -0.248622, besselj = -0.101506
(50) part = -0.097512, cumulative = -0.346134, besselj = -0.097512
(51) part = -0.004862, cumulative = -0.350996, besselj = -0.004862
(52) part = 0.090414, cumulative = -0.260582, besselj = 0.090414
(53) part = 0.101566, cumulative = -0.159017, besselj = 0.101566
(54) part = 0.020030, cumulative = -0.138986, besselj = 0.020030
(55) part = -0.078250, cumulative = -0.217237, besselj = -0.078250
(56) part = -0.103485, cumulative = -0.320721, besselj = -0.103485
(57) part = -0.033996, cumulative = -0.354717, besselj = -0.033996
(58) part = 0.065260, cumulative = -0.289458, besselj = 0.065260
(59) part = 0.103347, cumulative = -0.186111, besselj = 0.103347
(60) part = 0.046598, cumulative = -0.139513, besselj = 0.046598

[Mashrur Zawad]

I got it now.Thanks a lot

[Mashrur Zawad]

Another Thing if I use z:11:100 for example.After certain point I am getting mismatch in the result when it comes to 78.I am showing you the result from your code.Could you please tell is there anything I am missing

(77) part = 0.748634, cumulative = 0.463405, besselj = 0.066561

(78) part = 9.420296, cumulative = 9.883702, besselj = 0.087528

(79) part = 123.533271, cumulative = 133.416973, besselj = 0.028279

(80) part = 1581.984719, cumulative = 1715.401692, besselj = -0.056057

(81) part = 19631.756398, cumulative = 21347.158090, besselj = -0.088134

(82) part = 236185.355933, cumulative = 257532.514023, besselj = -0.039294

(83) part = 2756889.839186, cumulative = 3014422.353209, besselj = 0.044857

(84) part = 31244561.628974, cumulative = 34258983.982184, besselj = 0.087011

(85) part = 344049715.403874, cumulative = 378308699.386057, besselj = 0.049151

(86) part = 3683427773.000552, cumulative = 4061736472.386610, besselj = -0.033185

(87) part = 38366272816.132011, cumulative = 42428009288.518623, besselj = -0.084239

(88) part = 389032778462.679688, cumulative = 431460787751.198303, besselj = -0.057712

(89) part = 3842601154446.739746, cumulative = 4274061942197.937988, besselj = 0.021271

(90) part = 36993168098345.445312, cumulative = 41267230040543.382812, besselj = 0.079926

(91) part = 347312956427861.750000, cumulative = 388580186468405.125000, besselj = 0.064863

(92) part = 3181718541405431.500000, cumulative = 3570298727873836.500000, besselj = -0.009338

(93) part = 28456023042380144.000000, cumulative = 32026321770253980.000000, besselj = -0.074199

(94) part = 248588646214490144.000000, cumulative = 280614967984744128.000000, besselj = -0.070519

(95) part = 2122261570236178688.000000, cumulative = 2402876538220922880.000000, besselj = -0.002393

and goes on

[Walter Roberson]

Good question, and I do not know. When I examine the terms, they all appear to be reasonable to me, and since each is a symbolic rational, they sum to a symbolic rational, and I cannot see any reason why that rational would not be the 9.4-ish that is being reported. So I do not know what is going wrong.

Answer 2

[David Hill]

You need to increase your k as z gets larger.

n=sym(1);
F(1)=n;
for k=1:202
   F(k+1)=F(k)*k;
end
J=zeros(1,60);
for z=1:100
    JnE=sym(0);
    for k=0:200%this is suppose to be 0 to infinity! 
        JnE = JnE+(((-1)^k))/(F(k+1)*F(k+2))*(sym(num2str(z))/2)^(2*k+1);
    end
    J(z)=double(JnE);
    JJ(z)=besselj(1,z);
    J(z)-JJ(z)
end

[Mashrur Zawad]

I have increased the value of k and now I am getting good answer.Thank you

Answer 3

[John D'Errico]

Several common problems in this, regardless of whether you got the factiorial/gamma problem solved.

You want to use a series solution for a problem with large values of the argument. When the terms are alternating in sign, how well does that work? LOOK AT THE TERMS!

First, do you have the correct terms in the series?

https://en.wikipedia.org/wiki/Bessel_function

Next, NEVER use factorial on problems like this. In fact, a good idea is to NEVER use gamma either. Instead, form the terms by taking the log of each term. That will use the gammaln function. At least this way you can avoid overflows and underflows. Remember that a double precision number has a limited dynamic range.

Your task was to compute the series for besselj(1,z), where z is moderately large.

k = (0:100)';
z = 10:10:60; % I'll just write the problem using a few values for z.
n = 1;

% compute the natural logs of all terms, from 0:100
besseltermsln = -gammaln(k+1) - gammaln(k + n + 1) + (2*k+n)*log(z/2);

Now, remember, these terms will be exponentiated, and they alternate in sign.

plot(k,besseltermsln,'-')
legend('10','20','30','40','50','60')
grid on

Looking at that plot, you probably had a chance to compute the bessel function for z=10. But z=60? Even z=30 is probably questionable, but you might just make it.

Maybe what I need to do is to exponentiate the terms next.

semilogy(k,exp(besseltermsln),'-')
legend('10','20','30','40','50','60')
grid on

The problem will be massive subtractive cancellation. Adding and subtracting numbers, each of which is on the order of 1e16 or larger, means you loose all information content in the result. Now lets compute the sum of that series.

S = (-1).^k;
seriessum = cumsum(S.*exp(besseltermsln),1);

Now, for example, look at the series for z==10. The value should be

format long g
besselj(n,10)

ans = 
        0.0434727461688614

plot(seriessum(:,1))

so those alternating sums have terms as large as oughly 350. But the final series result should be tiny. That tells me you will never get a sum that is accurate to better than roughly 1e-12 or so.

seriessum(end,1)

ans = 
        0.0434727461676611

seriessum(end,1) - besselj(n,10)

ans = 
     -1.20026211192226e-12

Now, do the same computations for z=20. First, again I'll plot the series, to see how well it is converging.

plot(seriessum(:,2))

seriessum(end,2)

ans = 
        0.0668331175582756

besselj(n,20)

ans = 
        0.0668331241758499

seriessum(end,2) - besselj(n,20)

ans = 
     -6.61757428022103e-09

Here, we see that we will be losing roughly 7 of the lower order digits in the result.

And as I predicted, by the time you get to z=30, you will have only the top couple of the highest order digits correct.

[besselj(n,30),seriessum(end,3),seriessum(end,3) - besselj(n,30)]

ans = 1×3

-0.118751062616623 -0.118634922892277 0.000116139724346015

And finally, you should see that besselj(1,40) cannot be computed in double precision using that series as you want to do. The result is complete garbage at that point, and certainly beyond. Of course, this is probably the reason you were asked to do this homework assignment. You would need to use higher precision arithmetic to gain any hope of convergence.

[besselj(n,40),seriessum(end,4),seriessum(end,4) - besselj(n,40)]

ans = 1×3

0.126038318037585 21.3782131175142 21.2521747994766

[Mashrur Zawad]

Thanks for the insight.Now I have got some idea on this.

Answer 4

[the cyclist]

There are a handful of mistakes in your code:

• In one place, you used factorial instead of the gamma function. (You need to subtract 1 from the input argument, to use factorial).
• You should have set the starting value JnE to zero for each value of z (i.e. instead the loop over z).
• You have ignored the numerical precision that MATLAB can handle here. When the factorials and other factors get big, they are no longer stored precisely enough to be numerically stable.

In the code below, I have fixed the first two problems, and maxxed out your value of z to 30. You can see the two values are close, but also you are just starting to see the round-off error sneaking in.

n=1;
for z=11:30
    JnE = 0;
    for k=0:100
        JnE = JnE+(((-1)^k))/(factorial(k)*factorial(k+n))*(z/2)^(2*k+n);
    end
end
fprintf('JnEz=%f\n',JnE)

JnEz=-0.118764

n=1;
for z=11:30
    JnE=besselj(n,z);
end
fprintf('JnEz=%f\n',JnE)

JnEz=-0.118751

[Mashrur Zawad]

I got it now.Thanks a lot

[Mashrur Zawad]

I have tried with the code but when I increase the value of z for say z:11:100 after 78 I start getting big answers for the bessel series whereas the bessel function is giving good result.It is because of the accuracy of the matlab.I am giving some of the results of series from your code for your convenience

JnEz=-0.010977

JnEz=0.262996

JnEz=-2.437566

JnEz=-19.514872

JnEz=-9.709984

JnEz=-332.654043

JnEz=-389.703318

JnEz=-878.351073

JnEz=1446.201295

JnEz=5018.039351

JnEz=5637.749553

JnEz=46782.763903