Numerical round off / instability help? How to get more precision?

4 次查看（过去 30 天）
Brenden
Brenden 2013-3-13
I am trying to calculate the complex magnetic attenuation coefficient using the code below (for 3 different materials: Al, Cu & Steel):
clear all
close all
num_pt = 10000; % number of data points
fmax = 30000; % max frequency in Hz
f = (1:fmax/num_pt:fmax); % frequencies Hz
w = 2*pi.*f; % angular frequency
u = 4*pi*10^(-7); % permeiability of the metal
s = 1/(6.01349*10^(-8)); % Al condictivity of the metal
% s = 1/(2.37230*10^(-8)); % Cu condictivity of the metal
% s = 1/(1.55363*10^(-7)); % Steel condictivity of the metal
d0 = sqrt(2./(w.*u*s)); % skin depth
k0 = 1./d0;
R1 = 0.0181; % inner radius of pipe
R2 = 0.0193; % outer radius of pipe
v1 = k0.*R1*(1i+1);
v2 = k0.*R2*(1i+1);
% bessel functions of first and second kind for calculation of a the
% complex attenuation coefficient.
J0_v1 = besselj(0,v1);
J0_v2 = besselj(0,v2);
J1_v1 = besselj(1,v1);
J1_v2 = besselj(1,v2);
Y0_v1 = bessely(0,v1);
Y0_v2 = bessely(0,v2);
Y1_v1 = bessely(1,v1);
Y1_v2 = bessely(1,v2);
% break complex attenuation into parts for calculation
A = (J0_v1.*Y1_v1 - J1_v1.*Y0_v1);
B = (J0_v2.*Y1_v1 - J1_v1.*Y0_v2);
C = (k0.*R1./2);
D = (J0_v2.*Y0_v1 - J0_v1.*Y0_v2);
% complex attenuation coefficient
a = A./(B+C.*D.*(1+1i));
atacc = (a.*conj(a));
% high frequency approximation to a
% p = 2*sqrt(R2/R1).*exp(-k0.*(R2-R1))./sqrt(1+k0.*R1+k0.^2.*R1^2/2);
plot(f,atacc);
axis([0,fmax,0,1])
My issue is that I am getting numerical instability for high frequency. Running the above code you will see erratic spikes (starting around 15kHz) which should not be there. It is clear this is a precision issue as one can convert all values to 'single' and the erratic behaviour moves to lower frequencies. Is there anything that I can do to increase the precision to allow calculation at higher frequencies without error.
Thank you very much for your time and comments. BN

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2013-3-13
If you have the symbolic toolbox you can do the calculations in it.
Remember: your precision can be no better than the precision of your 6.01349*10^(-8) constant.
  3 个评论
显示 1更早的评论
Walter Roberson
Walter Roberson 2013-3-13
Remember standard numeric analysis: the precision of your answer can never be better than the precision of your inputs. We can presume that your "4" is infinitely precise and your "pi" is as precise as numerically possible, so looking at your other inputs, your lowest precision values are the R1 and R2, only 3 digits each.
For example your line
R1 = 0.0181;
means that R1 is between [0.01805, 0.01815). You use those low-precision radii to calculate v1 and v2, so v1 and v2 can be not better than 3 digits either. Does it then make sense to calculate your bessel functions to more than 3 digits? Well, possibly it does, but output of those bessel calculations should be no better than what you would get at the limits of what those 3 digits of accuracy of v1 and v2 could represent.
Doing a small test: by the time you get to v1(end) your precision for J0_v1 is down to 1 digit, times 10^9. Checking with a high-precision package, I see that MATLAB's output for besselj(0) near v1(end) is 14 digits of precision -- so if you could actually supply inputs more precise then 3 digits, you could get quite enough precision from your besselj calls.
GIGO. Garbage In, Garbage Out. Your input values do not deserve even as much accuracy as you think you are seeing.
Brenden
Brenden 2013-3-13
编辑：Brenden 2013-3-13
Yes I see what you are saying. If then I take measurements with higher precision I should postpone the effects of this instability.
Thanks again.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by