Plotting error for Bessel Function for "large" input values.

Sample of code in question:

%% Inputs

InputAngle = 4;
Wavelength = 0.65;

z = 10100;

FiberCoreRadius = 100;

%% Window
r = 0:10:10000;

%% Dependents
k = 2*pi/Wavelength;
AngleRad= InputAngle/180*pi;
w0 = FiberCoreRadius*0.67; % Assuming fiber is long enough to generate

Zr = k*w0^2/2;
wz = w0*sqrt(1+z^2/Zr^2);
Rz = (Zr^2+z^2)/(z);
alpha = k*sin(AngleRad);
l = alpha*FiberCoreRadius;

%% BG Code
Gaussian = exp(-(1/wz^2-1i*k/(2*Rz)).*(r.^2+alpha^2*z^2/k^2));
Bessel = besselj(l,(alpha.*r)./(1+1i*z/Zr));


FullBeamProfile = abs(Gaussian.*Bessel).^2;


%% Plotting

hold on;
plot(r,FullBeamProfile/max(FullBeamProfile),'b');

For larger Z values, the plot is returning zero for all values and not even attempting to plot larger values of r which is preventing me from modeling what I need to for my research. For smaller Z values there is no issue, but I need larger values to confirm what I'm working on.

The posted example gives the gaussian curve but does not go the full 10000 'r' distance I requested which becomes problematic when I need to see beyond that (for larger z values). It looks as though my Bessel function is going beyond the limits of a 64bit float and returning INF, and my gaussian function is returning below the limit, thus returning 0. I'm not 100% sure if this is causing this issue though.

Does anyone have any idea on how to fix this kind of error or have any suggested work arounds? I would love to have the full Gaussian curve generated, but if it helps, having the distance to peak value would work as well.

[Walter Roberson]

You should switch to symbolic calculation. The intermediate values in calculations such as you are doing can easily exceed what double precision can hold.

Answer 1

[David Goodmanson]

Hi Greg,

Time to look at logs. The bessel function has a scaled version that keeps it from blowing up, and you can look at the log of the whole works as is done below. The plot shows the (complex) log of (Gaussian*Bessel)^2 and is free of overflow or underflow problems. Exponentiating that quantity and taking abs will get rid of the imaginary part at which point only the real part matters. The real part has a maximum at about r = 700, but it drops down to around -3e4 at r = 10000. exp(-3e4) is so ridiculously small that I don't think this a double precision problem, really. It's more an issue of what is even physically meaningful.

%% Inputs
InputAngle = 4;
Wavelength = 0.65;
z = 10100;
FiberCoreRadius = 100;
%% Window
r = 0:10:10000;
%% Dependents
k = 2*pi/Wavelength;
AngleRad= InputAngle/180*pi;
w0 = FiberCoreRadius*0.67; % Assuming fiber is long enough to generate
Zr = k*w0^2/2;
wz = w0*sqrt(1+z^2/Zr^2);
Rz = (Zr^2+z^2)/(z);
alpha = k*sin(AngleRad);
l = alpha*FiberCoreRadius;
% -------------same up to this point ----------------
%% BG Code
G = exp(-(1/wz^2-1i*k/(2*Rz)).*(r.^2+alpha^2*z^2/k^2));
logG = (-(1/wz^2-1i*k/(2*Rz)).*(r.^2+alpha^2*z^2/k^2));
% Bessel = besselj(l,(alpha.*r)./(1+1i*z/Zr));
argB = (alpha.*r)./(1+1i*z/Zr);
logsfac = abs(imag(argB));     % scaled bessel is bessel divided by exp(logsfac) 
logBscaled = log(besselj(l,(alpha.*r)./(1+1i*z/Zr),1));       % note third input
%FullBeamProfile = abs(Gaussian.*Bessel).^2;
logFBP = 2*(logG +logsfac +logBscaled);
figure(1)
plot(r,real(logFBP),r,imag(logFBP))
return
   

[Greg Lovell]

I figured the logs wouldn't work since the Bessel function itself was reaching infinity rather quickly. Thanks for the tips and code.