Bessel Function Calculation in Matlab

34 次查看（过去 30 天）
Mashrur Zawad
Mashrur Zawad 2023-2-22
评论： Mashrur Zawad 2023-2-23
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

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2023-2-22
Ran in:
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
  3 个评论
显示 1更早的评论
Mashrur Zawad
Mashrur Zawad 2023-2-22
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
Walter Roberson 2023-2-22
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.

请先登录，再进行评论。

更多回答（3 个）

David Hill
David Hill 2023-2-22
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
John D'Errico
John D'Errico 2023-2-22
编辑：John D'Errico 2023-2-23
Ran in:
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
the cyclist
the cyclist 2023-2-22
Ran in:
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
  2 个评论
Mashrur Zawad
Mashrur Zawad 2023-2-22
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Bessel functions 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by