# How can integrate and differentiate spherical Bessel functions in Matlab

조회 수: 10(최근 30일)
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 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