Numerical round off / instability help? How to get more precision?
5 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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
0 Commenti
Risposta accettata
Walter Roberson
il 13 Mar 2013
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 Commenti
Walter Roberson
il 13 Mar 2013
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.
Più risposte (0)
Vedere anche
Prodotti
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!