12 views (last 30 days)

Show older comments

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

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

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

Start Hunting!