Really simple MATLAB code with illogical results

Question

Below is my MATLAB code to generate imaginary part of wavespeed (c), using function oscalcpcf when Reynolds number is 249 and I want to run for wavenumber (alpha) between 0.1 to 2 in steps of 2.

c1holder = [];
c2holder = [];
reyholder = [];
alphaholder = [];
wantedrey = [];
wantedalpha = [];
for Rey = 249
    for alpha = 0.1:0.1:2
        c1=oscalcpcf(Rey,alpha,100);
        c2=oscalcpcf(Rey,alpha,200);
        c1holder = [c1holder c1];
        c2holder = [c2holder c2];    
        reyholder = [reyholder Rey];
        alphaholder = [alphaholder alpha];      
    end
end
vectors = [c1holder' c2holder' reyholder' alphaholder'];

The above code is not difficult at all in my opinion but my laptop went cranky for some pairs of Reynolds number and alpha. Just naming one of them, Reynolds number = 249 and alpha = 0.3.

When I run the above code, I get c1 = 6.06002472332094E-08 and c2 = 0.0000010870344982811.

Now here's the problem. If I run from 2 to 0.1 with steps -0.1, I get c1 = -0.337584041016646 and c2 = 0.0000364854401656638.

AND if I were to check manually using oscalcpcf, ie oscalcpcf(249,0.3,100) and oscalcpcf(249,0.3,200), I get c1 = -0.337583911335139 and c2 = -0.337577395716528.

I have really no idea what's going on here can someone please help!

EDIT

alpha: 2.000000000000000000000000000000
alpha: 1.899999999999999900000000000000
alpha: 1.800000000000000000000000000000
alpha: 1.700000000000000000000000000000
alpha: 1.600000000000000100000000000000
alpha: 1.500000000000000000000000000000
alpha: 1.399999999999999900000000000000
alpha: 1.299999999999999800000000000000
alpha: 1.200000000000000000000000000000
alpha: 1.100000000000000100000000000000
alpha: 1.000000000000000000000000000000
alpha: 0.899999999999999910000000000000
alpha: 0.799999999999999820000000000000
alpha: 0.699999999999999960000000000000
alpha: 0.599999999999999870000000000000
alpha: 0.500000000000000000000000000000
alpha: 0.399999999999999910000000000000
alpha: 0.299999999999999820000000000000
alpha: 0.199999999999999960000000000000
alpha: 0.099999999999999867000000000000

and for 0.1 to 2

alpha: 0.100000000000000010000000000000
alpha: 0.200000000000000010000000000000
alpha: 0.300000000000000040000000000000
alpha: 0.400000000000000020000000000000
alpha: 0.500000000000000000000000000000
alpha: 0.599999999999999980000000000000
alpha: 0.700000000000000070000000000000
alpha: 0.800000000000000040000000000000
alpha: 0.900000000000000020000000000000
alpha: 1.000000000000000000000000000000
alpha: 1.100000000000000100000000000000
alpha: 1.200000000000000200000000000000
alpha: 1.300000000000000300000000000000
alpha: 1.400000000000000100000000000000
alpha: 1.500000000000000200000000000000
alpha: 1.600000000000000100000000000000
alpha: 1.700000000000000200000000000000
alpha: 1.800000000000000300000000000000
alpha: 1.900000000000000100000000000000
alpha: 2.000000000000000000000000000000

OMG, why doesn't my computer give precise steps of 0.1 when I told it to do so. The function oscalcpcf is very sensitive to small changes in the alpha and when I check these values that my script uses, it matches if I do it manually by oscalcpcf. Could you suggest a way for my computer to give precise steps of 0.1? Thank you.

Without knowing what `oscalcpcf` is doing we can only guess. Depending on how sensitive the function is to its input your problem could be precision. Try adding the line `disp(['alpha: ' sprintf('%.30f',alpha)]);`to your inner for-loop to see what the actual value of `alpha` is. — Deve, Sep 01 '12 at 08:26
Actually, perfectly logical results, once you understand how computers do arithmetic. — , Sep 01 '12 at 13:55

score 2 · Answer 1 · edited May 23 '17 at 10:34

2

I believe you have a floating point error because of your colon-generated 0.1:0.1:2 vector.

You get inaccurate values for alpha, because computers can't represent all numbers exactly given a fixed storage size such as double precision. This is especially true with the colon operator, which makes the inaccuracy propagate through the vector elements.

Now, I'm not sure whether this improves your results or not, but based on how the COLON operator works, I would suggest that you try to run your loop as following (similar to this answer of mine):

for alpha = ((1:20) / 10)

Also, if oscalcpcf() is a function which is allowed to be tampered with, I suggest you look into it and improve its robustness/sensitivity to small changes in the input. A 10^-14% inaccuracy should by no means have a significant impact on your results.

edited May 23 '17 at 10:34

Community

1
1

answered Sep 01 '12 at 09:10

Eitan T

32,660
14
72
109

hey thanks for the suggestion but inaccuracies still exist... I wonder if there's any way to create precise vectors of step size 0.1 at all.. – matlablearner Sep 01 '12 at 09:27
Unfortunately not. I suggest you work on `oscalcpcf` to be less sensitive. – Eitan T Sep 01 '12 at 09:49
3

NO. Read the responses. You cannot represent the number 0.1 EXACTLY in floating point arithmetic. – Sep 01 '12 at 10:56

score 1 · Answer 2 · edited Sep 01 '12 at 10:58

try this:

alpha = ((1:20) *1e-1). you should get :

0.100000000000000005551115123126 
0.200000000000000011102230246252 
0.300000000000000044408920985006 
0.400000000000000022204460492503 
0.500000000000000000000000000000 
0.600000000000000088817841970013 
0.700000000000000066613381477509 
0.800000000000000044408920985006 
0.900000000000000022204460492503 
1.000000000000000000000000000000 
1.100000000000000088817841970013 
1.200000000000000177635683940025 
1.300000000000000044408920985006 
1.400000000000000133226762955019 
1.500000000000000000000000000000 
1.600000000000000088817841970013 
1.700000000000000177635683940025 
1.800000000000000044408920985006 
1.900000000000000133226762955019 
2.000000000000000000000000000000

The results would be accurate to the double precision standard. Also, sym and vpa functions in matlab symbolic toolbox may help if u have it.

Really simple MATLAB code with illogical results

2 Answers2