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.
1 个评论
Walter Roberson
2019-1-14
You should switch to symbolic calculation. The intermediate values in calculations such as you are doing can easily exceed what double precision can hold.
采纳的回答
David Goodmanson
2019-1-14
编辑：David Goodmanson
2019-1-14
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
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Bessel functions 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
