Drawing Bessel function is not possible in that it give me an error

I need to draw First kind of Bessel Function. I have provided the following code but it does not work.

%Bessel polynomial
%Bessel polynomial of the first kind
x=0:0.01:5
for m=0:5
    i=0:1*10^6
    j=sum((((-1).^i)/(factorial(i).*factorial(gamma(i+m+1)))).*(x/2).^(2.*i+m))
    hold on
    plot(x,j)
end
legend('J0','J1','J2','J3','J4')
    

It gives me the following error:

错误使用  .^ 
对于此运算，数组的大小不兼容。

出错 differential_equation (第 27 行)
    j=sum((((-1).^i)/(factorial(i).*factorial(gamma(i+m+1)))).*(x/2).^(2.*i+m))

    

Please help

[Torsten]

You want to compute the factorial of 10^6 and factorial(gamma(1e6)) ? Did you think about how large these numbers will be ?

Why don't you use the implemented MATLAB functions, e.g. besselj, bessely, besselh, besseli and besselk ?

[Zeyuan]

I will look into the functions you mention. But I still would like to know what is wrong with the function I proposed. I followed strictly to the first kind of bessel function. I followed this equation:

Answer 1

[Torsten]

Maybe this is what you want:

x=(0:0.01:5).';
i=0:20;
hold on
for m=0:5
    j=sum((-1).^i./(factorial(i).*gamma(i+m+1)).*(x/2).^(2*i+m),2);
    plot(x,j)
end
legend('J0','J1','J2','J3','J4','J5')
grid on

●

Answer 2

[Walter Roberson]

x=0:0.01:5

That is a row vector of length 501.

i=0:1*10^6

That is a row vector of length 1000001

j=sum((((-1).^i)/(factorial(i).*factorial(gamma(i+m+1)))).*(x/2).^(2.*i+m))

That (-1).^i is a row vector of length 1000001

The factorial(i) is a row vector of length 1000001; the factorial(gamma(i+m+1)) is the same length. Using .* between two row vectors the same size is valid in MATLAB, so the factorial(i).*factorial(gamma(i+m+1)) subexpression is valid at the MATLAB level.

Now, factorial(1000001) is far beyond what can be calculated in double precision; but it is possible if we were to switch over to symbolic calculations.

gamma(i+m+1) would be the same as factorial(i+m) because i and m are both non-negative integers. With the largest i and largest m, the gamma() term would be as much as factorial(1000006)

You then factorial() the gamma() result. The fact that you call factorial() here tells us that you expect the gamma() result to be a non-negative integer... which makes it even more strange that you called gamma() instead of factorial(). But that is "odd" rather than "wrong".

Now, for the largest i and largest m value, you are trying to calculate factorial(factorial(1000006)). That is going to take a rather long time if it is possible at all, considering that factorial(1000006) is roughly 8 * 10^5565708 . But let us imagine that it is indeed possible and continue on with the analysis.

You now have the (-1).^i row vector of length 1000001 and that is / against a row vector of length 1000001 . Now, technically that is valid MATLAB, giving a scalar least squared solution to X * (the first vector) = (the second vector) -- but is that really what is wanted??

Unfortunately, if you have already switched to symbolic work in order to be able to represent factorial(1000000), you run into the problem that the symbolic / operator cannot handle the case of (rowvector / rowvector), returning a warning. For the case of rowvector / twodimensionalmatrix the symbolic / operator returns infinity as many times as there are rows in the two dimensional matrix if there are fewer rows than there are elements in the leading row vector. For the case where the two dimensional matrix is square with as many rows and columns as is in the leading row vector, then the symbolic / operator can produce meaningful results.

So... if you really do mean the / matrix-right-divide operator between those two large vectors, then you will need to impliment it differently in order to use symbolic numbers.

After all of that comes the (x/2).^(2.*i+m) term. The x is a row vector of length 501. The (2.*i+m) is a row vector of length 1000001 . When you have non-scalar operator to the left side of .^ then the right side needs to be either a scalar or else a matrix that is a "compatible" size . [1 501] is not compatible with [1 1000001] .

I really cannot tell what your intended calculation is. Maybe you intend to use x as a column vector and to use the ./ operator instead of the / operator; in which case the result j would be 501 x 1000001. It is possible to plot a vector of length 501 against an array 501 x 1000001 -- slow, but possible. Kinda meaningless in terms of density of information considering limited screen size, but possible.

[Torsten]

It's intended to implement the infinite series for J_alpha(x) under

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

But the OP should carefully read your explanations - they show how much can go wrong in a single line of code.

[Zeyuan]

x=(0:0.01:15).';
i=0:50;
hold on
set(gca,'YLim',[-0.5,1.2]); %y坐标显示限制(视窗)
set(gca,'XLim',[0,15]); %x坐标显示限制(视窗)
for m=0:2
    j=sum((((-1).^i)./(factorial(i).*factorial(i+m))).*(x/2).^(2*i+m),2);
    plot(x,j)
end
legend('J0','J1','J2')

Thank you for your comment, that is very helpful. I updated my code to the above. It works well now.

[Walter Roberson]

syms x
syms m integer
assume(m >= 0)
syms i integer
assume(i >= 0)
syms K
symsum((((-1).^i)./(factorial(i).*factorial(i+m))).*(x/2).^(2*i+m), i, 0, inf)

ans = 

[Zeyuan]

https://www.mathworks.com/matlabcentral/answers/2175778-bessel-funtion-drawing-problem-part-two?s_tid=srchtitle

Can you please help me take a look at this?