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

7 visualizzazioni (ultimi 30 giorni)

Modificato: Torsten il 26 Lug 2024
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 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

### Risposta accettata

Torsten il 26 Lug 2024
Modificato: Torsten il 26 Lug 2024
It works for the numerical "besselj" function.
Its derivative is given here:
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 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

### Più risposte (1)

Torsten il 26 Lug 2024
Modificato: Torsten il 26 Lug 2024
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 CommentiMostra 1 commento meno recenteNascondi 1 commento meno recente
Torsten il 26 Lug 2024
Modificato: Torsten il 26 Lug 2024
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)

Accedi per commentare.

### Categorie

Scopri di più su Bessel functions in Help Center e File Exchange

R2024a

### Community Treasure Hunt

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

Start Hunting!

Translated by