Bessel Function Calculation in Matlab
17 次查看（过去 30 天）
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
JnE = JnE+(((-1)^k))/(factorial(k)*factorial(k+n+1))*(z/2)^(2*k+n);
Can anyone tell me where I am doing it wrong?
I am also giving my code by using besselj function
Walter Roberson 2023-2-22
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.
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);
David Hill 2023-2-22
You need to increase your k as z gets larger.
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);
John D'Errico 2023-2-22
编辑：John D'Errico 2023-2-23
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?
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.
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.
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
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) - besselj(n,10)
Now, do the same computations for z=20. First, again I'll plot the series, to see how well it is converging.
seriessum(end,2) - besselj(n,20)
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)]
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)]
the cyclist 2023-2-22
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.
JnE = 0;
JnE = JnE+(((-1)^k))/(factorial(k)*factorial(k+n))*(z/2)^(2*k+n);