CODEWITHSUNDEEP
matlab
matlablearner
matlablearner

Reputation: 21

Really simple MATLAB code with illogical results

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.

Upvotes: 2

Views: 193

Answers (2)

oalah
oalah

Reputation: 89

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.

Upvotes: 1

Eitan T
Eitan T

Reputation: 32920

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.

Upvotes: 2

Related Questions