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

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

Answer 1

[Walter Roberson]

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.

[Brenden]

Unfortunately I do not have that tool box (however I will have to get it!). So are you saying that if I increase the precision of the constant then I will prevent the instability? Can I do this by any other means?

Also why is it that my precision can not be better than 6.01349*10^(-8)?

[Walter Roberson]

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]

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.