MATLAB Answers

How can integrate and differentiate spherical Bessel functions in Matlab

조회 수: 10(최근 30일)
M Bakry
M Bakry 2021년 5월 23일
답변: David Goodmanson 2021년 5월 26일
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

답변(1개)

David Goodmanson
David Goodmanson 2021년 5월 26일
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