Which is the right way of initialization of variables in this code?

I have been solving a set of 5 equations with 5 unknowns by lsqnonlin in a 'for' loop as given below:

for i=35:50
        f1= @(w) [0.001*(0.9-0.3*Pa/w(2))*Av*Pa*sqrt(2*Ga/(R*T*(Ga-1))*(Pin/Pa)^((Ga-1)/Ga)*((w(2)/Pa)^((Ga-1)/Ga)-1))-w(1);
             (Mi-(w(1)+mv(i-1)))*R*T/(M*w(3))-w(2);
             (Mi-(w(1)+mv(i-1)))*R*T/(M*Pa)*(1+0.1*((w(2)-Pa)*0.000145)^.623)-w(3);
             pi*w(4)*w(5)*L/4-w(3);
             pi/4*(w(5)+w(4))*(1+3*((w(5)-w(4))/(w(5)+w(4)))^2/(10+sqrt(4-3*((w(5)-w(4))/(w(5)+w(4)))^2)))-pi*r/2];
        [x,residual,resnorm,exit,output]= lsqnonlin(f1,[0.00015,P,V,a,b],[0,0,0,0,r],[Mi,P,V,a,Inf],optimset('Display','off','MaxFunEvals', 400000,'MaxIter',40000))
        mv(i)=x(1)+mv(i-1)
        P=x(2)
        V=x(3)
        a=x(4) 
        b=x(5)
    end

Apart from unknown variables which are found, the code requires values of other variables, and initial values and bounds supplied for the function. In initializing these values, I have found that I get different results based on how these are initialized. Below are the two ways I have been confronted with:

Method 1:
    d=.035;Av=2*pi*d^2/4;
    L=.5;r=.38;a=r;b=r;V=pi*a*b*L/4;
    Mi=0.05;R=8.314;T=550;M=0.0289;P=Mi*R*T/(M*V);Pin=P;
    mv(34)=0;
    Pa=101325;Ga=1.4;

Method 2:
    d=.035;Av=0.00192;
    L=.5;r=.38;a=0.38;b=0.38;V=0.0567;
    Mi=0.05;R=8.314;T=550;M=0.0289;P=139528;Pin=139528;
    mv(34)=0;
    Pa=101325;Ga=1.4;

The only difference in both methods is that:

In method 1, certain variables are found by some formulation.

In method 2, all variables are supplied as numeric values.

But the values of the initialized variables are the same irrespective of the methods. But upon solving with the values supplied by these two methods, the results are different. I have tried debugging to find the reason, but in vain. I know that the reason might be pretty simple, but I'm unable to find it. Someone please give which method is correct in this context. I herewith also attach the code for these two methods.

Thanks in advance.

[Adam]

How different are your "different results"?

[Thanigaivel Raja T]

Method 1 result at i=50(end):

mv =

Columns 1 through 8

         0         0         0         0         0         0         0         0

Columns 9 through 16

         0         0         0         0         0         0         0         0

Columns 17 through 24

         0         0         0         0         0         0         0         0

Columns 25 through 32

         0         0         0         0         0         0         0         0

Columns 33 through 40

         0         0    0.0001    0.0003    0.0005    0.0006    0.0008    0.0009

Columns 41 through 48

    0.0011    0.0012    0.0014    0.0015    0.0017    0.0018    0.0020    0.0021

Columns 49 through 50

    0.0023    0.0024

P =

   1.3863e+05

V =

    0.0543

a =

    0.2635

b =

    0.4828

Method 2:

mv =

Columns 1 through 8

         0         0         0         0         0         0         0         0

Columns 9 through 16

         0         0         0         0         0         0         0         0

Columns 17 through 24

         0         0         0         0         0         0         0         0

Columns 25 through 32

         0         0         0         0         0         0         0         0

Columns 33 through 40

         0         0    0.0001    0.0003    0.0005    0.0006    0.0008    0.0009

Columns 41 through 48

    0.0011    0.0012    0.0014    0.0015    0.0017    0.0018    0.0020    0.0021

Columns 49 through 50

    0.0023    0.0024

P =

   1.3604e+05

V =

    0.0553

a =

    0.3303

b =

    0.4266

Similarly, the value of P,V,a,b are different for every value of i.

Answer 1

[John D'Errico]

Congratulations!!!!!! You are the 10 millionth person to learn about floating point numbers! First prize for having learned that is a free trip to Newark, New jersey. By free, I mean that we pay nothing, you pay for the travel, hotel, food, etc.

Ok, seriously, look at what you did.

format long g
V=pi*a*b*L/4
V =
        0.0567057473972958

So in fact, V is not 0.0567, but a number with more significant digits. There is a difference. Those extra digits are very often significant, as you have found here.

Which is correct? Obviously you need to not cut off all of your numbers at 3 significant digits. Is 1/3 really 0.333? Of course not. Nor does pi=3.14 (unless you live in certain backwater places in the US, or for US politicians where it may be only 3.) Using such poor approximations will almost always get you in trouble.

Enjoy your trip to the vacation hot spot of Newark. :) (I'll admit I grew up in New Jersey, not far from Newark.)

[Thanigaivel Raja T]

Thank you for your suggestion. There is some considerable difference in values with double precision.

Based on supplying the numeric values upto double precision, the initialized variables are:

Method 2:

    d=.035;Av=0.00192422550032375;
    L=.5;r=.38;a=0.38;b=0.38;V=0.0567057473972958;
    Mi=0.05;R=8.314;T=550;M=0.0289;P=139514.000570562;Pin=139514.000570562;
    mv(34)=0;
    Pa=101325;Ga=1.4;

But by solving again by both methods, the results are still largely different (especially value of variables a and b) as given below:

At i=50,

By method 1,

mv(50)=

0.00240792391703642

P =

          138627.299203453

V =

        0.0543201242795223

a =

         0.263535129795466

b =

         0.482824540830555

By method 2,

mv(50)=

0.00239369626785816

P =

          138005.633489696

V =

        0.0545811290217641

a =

         0.379503213649327

b =

         0.380000000002385

So I want to know does MATLAB does calculations at even higher accuracy and hence the difference in results or there is something else that needs to be fixed in the code?

I herewith also attach Method 2 in updated format.

[Thanigaivel Raja T]

Upon increasing the step size and reducing the function tolerance from default value, i get approximately similar results.

[x,residual,resnorm,exit,output]= lsqnonlin(f1,[0.00015,P,V,a,b],[0,0,0,0,r],[Mi,P,V,a,Inf],optimset('FinDiffRelStep',1e-2,'Tolfun',1e-12,'Display','off','MaxFunEvals', 400000,'MaxIter',40000))

Method 1 result:

mv(50)=

0.00239636426157751

P =

           136128.86080041

V =

        0.0553305236043039

a =

          0.33016417835987

b =

         0.426751388297412

Method 2:

mv(50)=

0.00239620412514647

P =

            136297.4580337

V =

         0.055262266751178

a =

         0.328904597110712

b =

         0.427857220821132

Am I right now? Or something else need to be done?

[David Goodmanson]

Hardly seems possible, assuming that you have checked that all your initial variables are exactly equal between the two methods, and that lsqnonlin is deterministic.

[Thanigaivel Raja T]

Yes, looks like the problem is with approximation. Even if values are supplied at double precision, the results aren't the same. So I want to know, upto what precision does MATLAB hold? And why was I getting better closer results by increasing the step size and reducing function tolerance, even though the input values of double precision and actual input values aren't the same?

[John D'Errico]

You have two different starting points. Since you have not shown me CAREFULLY that they are in fact the same down to double precision, then you probably have not been careful to ensure that. All you have given us is the result. They are different down to about the 3rd digit. That leads me to believe that in fact, your parameters were also different down in the 3rd decimal place.

Double precision is 52 binary bits of precision, so ROUGHLY 16 digits.

As for your results due to a change in the step size and tolerances, how can I guess, since we have been given insufficient information from you as to what was really done here.

I'm sorry, but this happens over and over again. I see people saying they did something and got a strange result. In fact, what really happened is they did something completely different from what they claimed. People do sloppy things when they are starting to learn MATLAB. They have different values in a variable than they think they do. And since they don't really understand floating point arithmetic, or how numbers are stored, they look at 4 digits of output and think the number is exactly what they see.

So, IF you need help here, you need to supply the TRUE values of your variables, in a .mat file, attached to a comment. Then you need to show the explicit code that you are comparing. Show the results that you got. I'll take your code,with your variables, and run the same tests, to verify that you really did run the tests you claim to have done. Only then can I explain to you what you did and what is happening. Or you can do the same things I would do. I might put debug stops in your objective function, to look carefully at the iterations. I would turn the display on, look at the iterations, to see if the optimizer is converging happily, and how tolerances are impacting the result. Again, you can do all of the above.

Or you don't have to do any of the above, in which case, you would simply need to accept that yes, things will happen when you are not sufficiently careful about tolerances, about using numbers to their full precision.