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

Question

泽川 2024-7-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2140396-bessel-function-has-problems-in-converting-symbolic-function-into-handle-function

编辑： Torsten 2024-7-26

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2024-7-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2140396-bessel-function-has-problems-in-converting-symbolic-function-into-handle-function#answer_1491046

编辑：Torsten 2024-7-26

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Torsten 2024-7-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2140396-bessel-function-has-problems-in-converting-symbolic-function-into-handle-function#answer_1490976

编辑：Torsten 2024-7-26

在 MATLAB Online 中打开

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

Torsten 2024-7-26

编辑：Torsten 2024-7-26

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.

泽川 2024-7-26

在 MATLAB Online 中打开

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)