How can integrate and differentiate spherical Bessel functions in Matlab

15 views (last 30 days)
I have this integration problem and want to solve it numerically using Matlab
where κ is the Wavenumber.
I wrote the following Matlab code:
% the varaibles a, b and k must be defined before the follwong code
syms r
fun1=K*r*sqrt(pi./(2*K*r)).*besselj(n+0.5,K*r);
diffj=diff(fun1,r);
fun2=r^2*diffj;
integ_j=int(fun2, r, [a b]);
sol=double(integ_j);
when I run this code, i do not get the expected results.
Secondly, when the values of a or b go to zero, the result is NaN, due to the existence of r in the denomenator of the bessel function.
I would be grateful if someone could help me solve this problem or give me some hints how to differentiate or integrate spherical bessel functions

Answers (2)

David Goodmanson
David Goodmanson on 26 May 2021
Hello MB,
[1] Dimensionally, k can't be wavelength. Maybe you meant the wave number, 2pi / lambda.
[2] The integrand increases proportionally to r^2, which is certainly possible, but somewhat suspicious. Are you certain that an r^2 factor for the volume element has not already been inserted into ( kr j_n(kr) )^2 as the factors of r there? That happens sometimes.
You can use the identity
d/dz Jn(z) = -J(n+1,z) + (n/z)*Jn(z)
along with the identity you have already
jn(z) = sqrt((pi/2)/z) besselj(n+1/2,z)
to obtain
d/d[kr] (kr j_n(kr)) = -kr j_n+1(kr) + (n+1) j_n(kr)
The following function is the integrand, which can be plugged into the 'integral' function, integrating from a to b:
function y = fun(r)
k = 2.2; n = 3; % for example
kr = k*r;
y = ((-kr.*js(n+1,kr) + (n+1)*js(n,kr))).^2.*r.^2;
function y = js(n,x) % nested function
% spherical bessel function j, first kind
% y = js(n,x);
y = sqrt((pi/2)./x).*besselj(n+1/2,x);
% get rid of nans
if n==0
y(x==0) = 1;
elseif n >0
y(x==0) = 0;
else
y(x==0) = (-1)^(n-1)*inf;
end
end % end js
end % end fun

Hussein Thary
Hussein Thary on 17 Jun 2023
int
unction y = fun(r)
k = 2.2; n = 3; % for example
kr = k*r;
y = ((-kr.*js(n+1,kr) + (n+1)*js(n,kr))).^2.*r.^2;
function y = js(n,x) % nested function
% spherical bessel function j, first kind
% y = js(n,x);
y = sqrt((pi/2)./x).*besselj(n+1/2,x);
% get rid of nans
if n==0
y(x==0) = 1;
elseif n >0
y(x==0) = 0;
else
y(x==0) = (-1)^(n-1)*inf;
end
end % end js
end % end fun
Function definitions in a script must appear at the end of the file.
Move all statements after the "js" function definition to before the first local function definition.

Error in connector.internal.fevalMatlab

Error in connector.internal.fevalJSON

Community Treasure Hunt

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

Start Hunting!