Bessel function has problems in converting symbolic function into handle function.

I want to derive the spherical Bessel function, so I first construct the spherical Bessel function by using Bessel function.

,

Because I want to derive the spherical Bessel function, I want to first construct a symbolic expression of the spherical Bessel function, and then use the diff function to get the first derivative, and then convert it into a function handle.

However, when n is large, the symbolic expression of spherical Bessel function has problems（In the figure below, I haven't derived the spherical Bessel function yet, but there has been a case of non-convergence, so it can be seen that the derivative is also non-convergence）：

If I don't need to find the derivative of Bessel function, then I only need to use the handle to complete this task. I need to construct a function to find the first derivative of spherical Bessel function at any point. Is there any good way?

Answer 1

[Torsten]

It works for the numerical "besselj" function.

Its derivative is given here:

https://uk.mathworks.com/matlabcentral/answers/93677-how-can-i-evaluate-the-derivatives-of-a-bessel-function-at-different-points

n = 40;
x = 0:0.1:100;
y = sqrt(pi./(2*x)).*besselj(n+0.5,x);
% Derivative of y with respect to x
dy = -sqrt(pi)./(2*x).^(1.5).*besselj(n+0.5,x)+sqrt(pi./(2*x))*0.5.*(besselj(n+0.5-1,x)-besselj(n+0.5+1,x));
figure(1)
plot(x,[y;dy])
grid on

Or you work with symbolic precision:

syms x
y = sqrt(pi/(2*x))*besselj(n+0.5,x);
dy = diff(y,x);
xnum = 0.1:0.1:100;
ynum = subs(y,x,xnum);
dynum = subs(dy,x,xnum);
figure(2)
plot(xnum,[ynum;dynum])
grid on

Answer 2

[Torsten]

This code seems to work. Can you show us where you encounter problems ?

syms x
n = 12;
f = sqrt(sym(pi)/(2*x))*besselj(n+1/2,x)

f = 

df = diff(f,x)

df = 

dfnum = matlabFunction(df)

dfnum = function_handle with value:
    @(x)sqrt(2.0).*1.0./sqrt(x).*sqrt(1.0./(x.*2.0)).*(cos(x).*(1.0./x.^2.*-3.003e+3+1.0./x.^4.*1.35135e+6-1.0./x.^6.*1.9297278e+8+1.0./x.^8.*9.820936125e+9-1.0./x.^10.*1.51242416325e+11+1.0./x.^12.*3.16234143225e+11+1.0)+sin(x).*(1.0./x.^3.*6.006e+3-1.0./x.^5.*5.4054e+6+1.0./x.^7.*1.15783668e+9-1.0./x.^9.*7.8567489e+10+1.0./x.^11.*1.51242416325e+12-1.0./x.^13.*3.7948097187e+12)+(cos(x).*(1.0./x.^3.*1.925e+3-1.0./x.^5.*9.4248e+5+1.0./x.^7.*1.2087306e+8-1.0./x.^9.*4.700619e+9+1.0./x.^11.*4.0542838875e+10).*7.8e+1)./x+1.0./x.^2.*cos(x).*(1.0./x.^2.*9.625e+2-1.0./x.^4.*2.3562e+5+1.0./x.^6.*2.014551e+7-1.0./x.^8.*5.87577375e+8+1.0./x.^10.*4.0542838875e+9-1.0).*7.8e+1+(sin(x).*(1.0./x.^2.*9.625e+2-1.0./x.^4.*2.3562e+5+1.0./x.^6.*2.014551e+7-1.0./x.^8.*5.87577375e+8+1.0./x.^10.*4.0542838875e+9-1.0).*7.8e+1)./x)-(sqrt(2.0).*1.0./x.^(3.0./2.0).*(sin(x).*(1.0./x.^2.*-3.003e+3+1.0./x.^4.*1.35135e+6-1.0./x.^6.*1.9297278e+8+1.0./x.^8.*9.820936125e+9-1.0./x.^10.*1.51242416325e+11+1.0./x.^12.*3.16234143225e+11+1.0)-(cos(x).*(1.0./x.^2.*9.625e+2-1.0./x.^4.*2.3562e+5+1.0./x.^6.*2.014551e+7-1.0./x.^8.*5.87577375e+8+1.0./x.^10.*4.0542838875e+9-1.0).*7.8e+1)./x).*sqrt(1.0./(x.*2.0)))./2.0-(sqrt(2.0).*1.0./x.^(5.0./2.0).*(sin(x).*(1.0./x.^2.*-3.003e+3+1.0./x.^4.*1.35135e+6-1.0./x.^6.*1.9297278e+8+1.0./x.^8.*9.820936125e+9-1.0./x.^10.*1.51242416325e+11+1.0./x.^12.*3.16234143225e+11+1.0)-(cos(x).*(1.0./x.^2.*9.625e+2-1.0./x.^4.*2.3562e+5+1.0./x.^6.*2.014551e+7-1.0./x.^8.*5.87577375e+8+1.0./x.^10.*4.0542838875e+9-1.0).*7.8e+1)./x).*1.0./sqrt(1.0./(x.*2.0)))./4.0

[泽川]

After I plot the function, something is obviously wrong:

[Torsten]

We cannot make use of graphics. Can you include your code as plain .m-file by using the Code-button from the menu bar ?

According to the formula, your function should have a pole at x = 0, shouldn't it ? Maybe that's what the plot shows.

[泽川]

n = 40;
x = 0:0.5:100;
syms X
Formfunc   = sqrt(pi./(2*X)) .* besselj(n+0.5, X)


% Formfunc   = @(X)sqrt(pi./(2*X)) .* besselj(n+0.5, X);
% Formfunc   = besselj(n, X);

% Formfunc_d = diff(Formfunc, X, order);
Formfunc_d = Formfunc;
Formfunc_d = matlabFunction(Formfunc_d);

% 判断阶数n是否为非负整数
if ~isnumeric(n) || n < 0 || floor(n) ~= n
    error('阶数n必须为非负整数。');
end

y = zeros(size(x));   % 预分配结果数组

idx_zero = (x == 0);  % 筛选出 x=0 的索引位置
if any(idx_zero)
    y(idx_zero & (n == 0)) = 1;  % n=0 且 x=0 时，球状 Bessel 函数的值为 1
    y(idx_zero & (n ~= 0)) = 0;  % n 不等于 0 且 x=0 时，球状 Bessel 函数的值为 0
end

idx_nonzero = ~idx_zero;  % 筛选出 x!=0 的索引位置
if any(idx_nonzero)
    % 计算球状 Bessel 函数
    y(idx_nonzero) = Formfunc_d( x(idx_nonzero) );
end

plot(x,y)