How can integrate and differentiate spherical Bessel functions in Matlab

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

Answer 1

[David Goodmanson]

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 

Answer 2

[Hussein Thary]

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