problem using fplot for small arguments of spherical bessel functions

i'm trying to plot spherical bessel functions (i. e. half integer order bessel functions) for small arguments (x << 1) and i've tried using vectors and it turns out very well; but the problem is when i try to use symbolic variable to define the function and plot it, i observe like meshing problems or diverging behavior (like in the image below). is there a way to fix it?

i've already tried modify the meshdensity but the problem remains. thanks a lot

l = 6; % order of the bessel function
syms r
rho = sqrt(pi./(2*r)) .* besselj(l + 0.5,r); % definition of spherical bessel function 
% in terms of bessel function
f1 = figure;
fplot(real(rho),[0 0.35])
ylim([-0.1,0.5])

Answer 1

[David Goodmanson]

Hi Donaldo,

As you probably know, the function works pretty well except for small values of r around the origin where you see a bunch of numerical chaff. The problem here is with sym's choice of an expression for besselj:

rho =
(2^(1/2)*(1/(2*r))^(1/2)*(sin(r)*(210/r^2 - 4725/r^4 + 10395/r^6 - 1) ...
                    - (21*cos(r)*(495/r^4 - 60/r^2 + 1))/r))/r^(1/2)

The expression contains a sin(rr) and a cos(rr) function, each multiplied by a polynomial in 1/r^2. This is exact analytically, and works for most values of r. But computationally it makes no sense at all to try this when r very is small. In that case the besselj result itself is also very small. But since (1/r^2) is very large, the expression has a polynomial with huge terms of both signs, which have to all (almost) cancel out in order to get the result. Numerically, for small enough r it won't happen.

(When r is small, there is a different expression available, a power series solution in r^2).

For sym, it's not part of the job to recognize when r is getting too small and the expression breaks down numerically. So no one is on the job. Part of what you see in the plot is digitization noise when the expression overruns the 15 digits or so of double precision arithmetic.

Unfortunately PS's code plugs the sym same expression into a function and ends up having the same issues on a plot.

For small r you could compensate for the limitations of the sym expression by using variable precision arithmetic, vpa. But the expression itself is being used in a regime it is not intended for and it's not a good approach to patch things up by throwing a ton of decimal digits at it. Besides, it will never work for r = 0.

As you have seen by using vectors, numerical besselj is on the job and knows how to work the problem. It just makes sense to go with numerical besselj in this situation.

Answer 2

[prabhat kumar sharma]

Hello Donaldo,

To plot spherical Bessel functions for small arguments using symbolic variables in MATLAB, try the following approach to avoid issues:

1. Avoid Zero: Start your range slightly above zero to prevent division by zero.
2. Use Numeric Evaluation: Convert the symbolic expression to a numeric function and plot it:

l = 6; % order of the Bessel function
syms r
rho_sym = sqrt(pi./(2*r)) .* besselj(l + 0.5, r);
rho_func = matlabFunction(rho_sym, 'Vars', r);

r_values = linspace(0.01, 0.35, 1000); % Start slightly above zero
rho_values = arrayfun(rho_func, r_values);

plot(r_values, real(rho_values));
ylim([-0.1, 0.5]);

●

This should help eliminate the meshing issues and provide a smooth plot.