MATLAB Answers

How can integrate and differentiate spherical Bessel functions in Matlab

8 views (last 30 days)
M Bakry
M Bakry on 23 May 2021
Answered: David Goodmanson on 26 May 2021
I have this integration problem and want to solve it numerically using Matlab
where κ is the wavelength.
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 goes 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 (1)

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

Community Treasure Hunt

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

Start Hunting!

Translated by