Matlab for loop, different values in Python

Question

I have code in matlab and Im trying "rewrite" it into Python. The problem is with different values of for loop. Unfortunately I can't find a solution. Problem is with value in every array: nu theta mu M from i element (32 in matlab and 31 in python).

For both codes n=30

32 element in matlab = 31 element in python

In Python DTOR or RTOD are degrees to radians and vice versa.

Here is the whole python code

-------------------------------------Matlab:-------------------------------------


%here is problematic for loop
for t = 1:2
    i=n+2;
    h = i;
    k=0;
    z=2;
        for j= 1:n-1
            while (i<=h+n-k-1)
                if (i==h)
                    if (i==node_number-1)
                        y(i)=0;
                        theta(i)=0;
                        M(i) = Me;
                        nu(i) =PrandtlMeyer(M(i),gamma);
                        mu(i)= asin(1/M(i))*180/pi;
                        mI = tand((mu(i-n+k)-theta(i-n+k)+mu(i))/2);
                        x(i)= (y(i-n+k)/mI)+x(i-n+k);
                        if (t==2)
                            plot([x(i-n+k) x(i)],[y(i-n+k) y(i)]);
                            hold on;
                        end
                    else 
%here is 32 element
                        y(i)=0;
                        theta(i)=0;
                        nu(i) = theta(i-n+k)+nu(i-n+k)+(1/(sqrt((M(i-n+k)^2)-1)-cotd(theta(i-n+k))));
                        M(i) = InversePrandtlMeyer(nu(i));
                        mu(i)= asin(1/M(i))*180/pi;
                        mI = tand((mu(i-n+k)-theta(i-n+k)+mu(i))/2);
                        x(i)= (y(i-n+k)/mI)+x(i-n+k);
%here ends 32 element
                        if (t==2)
                            plot([x(i-n+k) x(i)],[y(i-n+k) y(i)]);
                            hold on;
                        end
                    end
                else if (i==h+1)
                        if (i==node_number)
                            mI = tand((theta(i-n+k)+theta(i))/2);
                            mII =tand((mu(i-1)+theta(i-1)+theta(i)+mu(i))/2);
                            x(i) = (y(i-1)-y(i-n+k)+(x(i-n+k)*mI)-(x(i-1)*mII))/(mI-mII);
                            y(i) = y(i-1)+((x(i)-x(i-1))*mII);
                            y1(z) = y(i-1)+((x(i)-x(i-1))*mII);
                            x1(z) = (y(i-1)-y(i-n+k)+(x(i-n+k)*mI)-(x(i-1)*mII))/(mI-mII);
                            theta(i) = theta(i-1);
                            nu(i) = nu(i-1);
                            M(i) = M(i-1);
                            mu(i)= asin(1/M(i))*180/pi;
                            z=z+1;
                            if (t==2)
                                plot([x(i-1) x(i)],[y(i-1) y(i)]);
                                hold on;
                                plot([x(i-n+k) x(i)],[y(i-n+k) y(i)]);
                                hold on;
                            end
                        else
%here is 33 element
                            mI = tand((mu(i-n+k)-theta(i-n+k)+mu(i)-theta(i))/2);
                            mII = tand((mu(i)+theta(i)+mu(i-1))/2);
                            x(i) = (y(i-n+k)+(x(i-1)*mII)+(x(i-n+k)*mI))/(mII+mI);
                            y(i) = (x(i)-x(i-1))*mII;
                            nu(i) = (theta(i-n+k)+nu(i-n+k)+(nu(i-1)/2)+((1/(sqrt((M(i-n+k)^2)-1)-cotd(theta(i-n+k))))*((y(i)-y(i-n+k))/y(i-n+k))))/1.5;
                            theta(i)= (nu(i)-nu(i-1))/2;
                            M(i) = InversePrandtlMeyer(nu(i));
                            mu(i)= asin(1/M(i))*180/pi;
%here end 33 element
                            if (t==2)
                                plot([x(i-1) x(i)],[y(i-1) y(i)]);
                                hold on;
                                plot([x(i-n+k) x(i)],[y(i-n+k) y(i)]);
                                hold on;
                            end
                        end
                    else if (i==h+n-k-1)
                            mI = tand((theta(i-n+k)+theta(i))/2);
                            mII =tand((mu(i-1)+theta(i-1)+theta(i)+mu(i))/2);
                            x(i) = (y(i-1)-y(i-n+k)+(x(i-n+k)*mI)-(x(i-1)*mII))/(mI-mII);
                            y(i) = y(i-1)+((x(i)-x(i-1))*mII);
                            y1(z)=y(i-1)+((x(i)-x(i-1))*mII);
                            x1(z) = (y(i-1)-y(i-n+k)+(x(i-n+k)*mI)-(x(i-1)*mII))/(mI-mII);
                            theta(i) = theta(i-1);
                            nu(i) = nu(i-1);
                            M(i) = M(i-1);
                            mu(i)= asin(1/M(i))*180/pi;
                            z=z+1;
                            if (t==2)
                                plot([x(i-1) x(i)],[y(i-1) y(i)]);
                                hold on;
                                plot([x(i-n+k) x(i)],[y(i-n+k) y(i)]);
                                hold on;
                            end
                        else
                            mI = tand((mu(i-n+k)-theta(i-n+k)+mu(i)-theta(i))/2);
                            mII = tand((mu(i-1)+theta(i-1)+mu(i)+theta(i))/2);
                            x(i) = (y(i-n+k)-y(i-1)+(x(i-n+k)*mI)+(x(i-1)*mII))/(mII+mI);
                            y(i) = y(i-1)+((x(i)-x(i-1))*mII);
                            theta(i) = (theta(i-n+k)+nu(i-n+k)+theta(i-1)-nu(i-1)+((1/(sqrt((M(i-n+k)^2)-1)-cotd(theta(i-n+k))))*((y(i)-y(i-n+k))/y(i-n+k)))-((1/(sqrt((M(i-1)^2)-1)-cotd(theta(i-1))))*((y(i)-y(i-1))/y(i-1))))/2;
                            nu(i) = (theta(i-n+k)+nu(i-n+k)-theta(i-1)+nu(i-1)+((1/(sqrt((M(i-n+k)^2)-1)-cotd(theta(i-n+k))))*((y(i)-y(i-n+k))/y(i-n+k)))+((1/(sqrt((M(i-1)^2)-1)-cotd(theta(i-1))))*((y(i)-y(i-1))/y(i-1))))/2;
                            M(i) = InversePrandtlMeyer(nu(i));
                            mu(i)= asin(1/M(i))*180/pi;
                                if (t==2)
                                    plot([x(i-1) x(i)],[y(i-1) y(i)]);
                                    hold on;
                                    plot([x(i-n+k) x(i)],[y(i-n+k) y(i)]);
                                    hold on;
                                end
                        end
                    end
                end
                i = i+1;
            end
            k=k+1;
            h=i;
        end
end

And first few rows of Matlab "theta" values:

0.4397    0.8793    1.3190    1.7587    2.1983    2.6380    3.0776    3.5173
3.9570    4.3966    4.8363    5.2760    5.7156    6.1553    6.5949    7.0346
7.4743    7.9139    8.3536    8.7933    9.2329    9.6726   10.1122   10.5519
10.9916   11.4312   11.8709   12.3106   12.7502   13.1899   13.1899         0
0.3032    0.7453    1.1877    1.6301    2.0725    2.5149    2.9572    3.3995
3.8418    4.2840    4.7262    5.1684    5.6106    6.0528    6.4951    6.9373
7.3795    7.8218    8.2641    8.7065    9.1488    9.5913   10.0337   10.4762
10.9188   11.3614   11.8041   12.2469   12.2469         0    0.2977    0.7407
1.1839    1.6271    2.0702    2.5133    2.9562    3.3991    3.8420    4.2848
4.7276    5.1704    5.6131    6.0559    6.4986    6.9414    7.3842    7.8270

-------------------------------------Python:-------------------------------------

from collections import OrderedDict
old_settings = np.seterr(all='print')
OrderedDict(np.geterr())
OrderedDict([('divide', 'print'), ('over', 'print'), ('under', 'print'), ('invalid', 'print')])
np.int16(32000) * np.int16(3)

#here is problematic for loop
i = n + 1
h = i
k = 0
z = 1
for j in range(n-1):
    while i <= h + n - k-1:
        if i == h:
            if i == node_number -1:
                    y[i] = 0
                    theta[i] = 0
                    M[i] = Me
                    nu[i] = PM(M[i], gamma)
                    mu[i] = np.arcsin(1 / M[i]) * 180 / np.pi
                    mI = np.tan((mu[i - n + k] - theta[i - n + k] + mu[i]) / 2 * DTOR)
                    x[i] = (y[i - n + k] / mI) + x[i - n + k]
            else:
#here is 31 element
                    y[i] = 0
                    theta[i] = 0
                    nu[i] = theta[i - n + k] + nu[i - n + k] + (1 / (np.sqrt((M[i - n + k]**2) - 1) - np.cos((theta[i - n + k])*DTOR)))
                    M[i] = InversePrandtlMeyer(nu[i])
                    mu[i] = np.arcsin(1 / M[i]) * 180 / np.pi
                    mI = np.tan(((mu[i - n + k] - theta[i - n + k] + mu[i]) / 2) * DTOR)
                    x[i] = (y[i - n + k] / mI) + x[i - n + k]
#here ends 31 element
        elif i == h + 1:
            if i == node_number:
                    mI = np.tan((theta[i - n + k] + theta[i]) / 2 * DTOR)
                    mII = np.tan((mu[i - 1] + theta[i - 1] + theta[i] + mu[i]) / 2 * DTOR)
                    x[i] = (y[i - 1] - y[i - n + k] + (x[i - n + k] * mI) - (x[i - 1] * mII)) / (mI - mII)
                    y[i] = y[i - 1] + ((x[i] - x[i - 1]) * mII)
                    y1[z] = y[i - 1] + ((x[i] - x[i - 1]) * mII)
                    x1[z] = (y[i - 1] - y[i - n + k] + (x[i - n + k] * mI) - (x[i - 1] * mII)) / (mI - mII)
                    theta[i] = theta[i - 1]
                    nu[i] = nu[i - 1]
                    M[i] = M[i - 1]
                    mu[i] = np.arcsin(1 / M[i]) * 180 / np.pi
                    z = z + 1

            else:
#here is 32 element
                    mI = np.tan((mu[i - n + k] - theta[i - n + k] + mu[i] - theta[i]) / 2 * DTOR)
                    mII = np.tan((mu[i] + theta[i] + mu[i - 1]) / 2 * DTOR)
                    x[i] = (y[i - n + k] + (x[i - 1] * mII) + (x[i - n + k] * mI)) / (mII + mI)
                    y[i] = (x[i] - x[i - 1]) * mII
                    nu[i] = (theta[i - n + k] + nu[i - n + k] + (nu[i - 1] / 2) + ((1 / (np.sqrt((M[i - n + k] ** 2) - 1) - np.cos((theta[i - n + k])*DTOR))) * ((y[i] - y[i - n + k]) / y[i - n + k]))) / 1.5
                    theta[i] = (nu[i] - nu[i - 1]) / 2
                    M[i] = InversePrandtlMeyer(nu[i])
                    mu[i] = np.arcsin(1 / M[i]) * 180 / np.pi
#here end 32 element
        elif i == h + n - k-1:
                mI = np.tan((theta[i - n + k] + theta[i]) / 2 * DTOR)
                mII = np.tan((mu[i - 1] + theta[i - 1] + theta[i] + mu[i]) / 2 * DTOR)
                x[i] = (y[i - 1] - y[i - n + k] + (x[i - n + k] * mI) - (x[i - 1] * mII)) / (mI - mII)
                y[i] = y[i - 1] + ((x[i] - x[i - 1]) * mII)
                y1[z] = y[i - 1] + ((x[i] - x[i - 1]) * mII)
                x1[z] = (y[i - 1] - y[i - n + k] + (x[i - n + k] * mI) - (x[i - 1] * mII)) / (mI - mII)
                theta[i] = theta[i - 1]
                nu[i] = nu[i - 1]
                M[i] = M[i - 1]
                mu[i] = np.arcsin(1 / M[i]) * 180 / np.pi
                z = z + 1
        else:
                mI = np.tan((mu[i - n + k] - theta[i - n + k] + mu[i] - theta[i]) / 2 * DTOR)
                mII = np.tan((mu[i - 1] + theta[i - 1] + mu[i] + theta[i]) / 2 * DTOR)
                x[i] = (y[i - n + k] - y[i - 1] + (x[i - n + k] * mI) + (x[i - 1] * mII)) / (mII + mI)
                y[i] = y[i - 1] + ((x[i] - x[i - 1]) * mII)
                theta[i] = (theta[i - n + k] + nu[i - n + k] + theta[i - 1] - nu[i - 1] + ((1 / (np.sqrt((M[i - n + k] ** 2) - 1) - np.cos((theta[i - n + k]) * DTOR))) * ((y[i] - y[i - n + k]) / y[i - n + k])) - ((1 / (np.sqrt((M[i - 1] ** 2) - 1) - np.cos((theta[i - 1])*DTOR))) * ((y[i] - y[i - 1]) / y[i - 1]))) / 2
                nu[i] = (theta[i - n + k] + nu[i - n + k] - theta[i - 1] + nu[i - 1] + ((1 / (np.sqrt((M[i - n + k] ** 2) - 1) - np.cos((theta[i - n + k])*DTOR))) * ((y[i] - y[i - n + k]) / y[i - n + k])) + ((1 / (np.sqrt((M[i - 1] ** 2) - 1) - np.cos((theta[i - 1])*DTOR))) * ((y[i] - y[i - 1]) / y[i - 1]))) / 2
                M[i] = InversePrandtlMeyer(nu[i])
                mu[i] = np.arcsin(1 / M[i]) * 180 / np.pi
        i = i + 1
    k = k + 1
    h = i

And first few rows of Python "theta" values:

[ 0.43966268  0.87932536  1.31898804  1.75865072  2.1983134   2.63797608
  3.07763876  3.51730144  3.95696412  4.3966268   4.83628948  5.27595216
  5.71561484  6.15527752  6.5949402   7.03460288  7.47426556  7.91392824
  8.35359092  8.7932536   9.23291628  9.67257896 10.11224165 10.55190433
 10.99156701 11.43122969 11.87089237 12.31055505 12.75021773 13.18988041
 13.18988041  0.          1.12391865  2.04904177  2.75249233  3.37815396
  3.96686921  4.53558282  5.09310744  5.64490943  6.19500187  6.74690367
  7.30431746  7.87186621  8.45626406  9.06871626  9.73089387 10.49343203
 11.51563061 13.74816848  0.85949953         nan         nan         nan
         nan         nan         nan         nan         nan         nan
         nan  0.          1.48303468  2.8885035   3.74795106  4.46269741
  5.11686874  5.74107193  6.35015955  6.95312861  7.55674416  8.16737947
  8.79250326  9.44272463 10.13594268 10.90840157 11.85472483 13.39395428
 13.37483292         nan         nan

sequences start to differ after position 32. I guess something is wrong in one of the branches. I suggest to attach a tag indicating which branch had produced specific value — tstanisl, Nov 24 '19 at 20:06
@tstanisl i edited it with comments ```%``` in ```matlab``` and ```#``` in ```python``` — macaroni, Nov 24 '19 at 20:24
@hpaulj ```range(n-1)``` means from 0 to 28 in ```python``` so it iterate 29 times. in ```matlab``` ```for j= 1:n-1``` is from 1 to 29 so it also iterate 29 times — macaroni, Nov 24 '19 at 20:46
@JohanC it change nothig in results. But both ```for j in range(1,n)``` and ```for j in range(n-1)``` iterate 29 times — macaroni, Nov 24 '19 at 20:50
A minor speed point, for whole arrays `np.sin` is right, but for single values, `math.sin` is faster. Same for the other trig calcs. — hpaulj, Nov 24 '19 at 21:13
It would be nice if your python code was complete and runnable. For example you don't show how variables like `x` and `y` are initialized. — hpaulj, Nov 24 '19 at 21:14
@hpaulj here is the whole python code https://pastebin.com/PDrEJ5xZ — macaroni, Nov 24 '19 at 22:25
I don't know if it's relevant, but the argument to `InversePrandtlMeyer(nu[i])` is often `nan` (after j=7). The "Warning: invalid value encountered in double_scalars" warnings occur before that in `y=(nu/nu_0)**(2/3)` line. — hpaulj, Nov 25 '19 at 05:45
Also the argument to `np.sqrt((M[i - n + k]**2) - 1)` is often negative. — hpaulj, Nov 25 '19 at 06:00
All the `theta` calculated in `for j in range(n-1):` loop are wrong. — hpaulj, Nov 25 '19 at 06:07

score 0 · Answer 1 · answered Nov 24 '19 at 20:04

0

one of the problems with python calculations is that they are quiet often wrong then can be one of the reasons why it is not equal to eachother you can find more information about it here: Is floating point math broken?

answered Nov 24 '19 at 20:04

Matthijs990

637
3
26

2

is it an answer? – tstanisl Nov 24 '19 at 20:07
but until 32 element in ```matlab``` (31 in ```python```) everything is ok. I guess something is wrong with ```i``` element numbering, because loop starts from here – macaroni Nov 24 '19 at 20:10
@tstanisl it is a sugestion – Matthijs990 Nov 24 '19 at 20:22
My first guess would indeed be rounding error due to floating point arithmetic as the results seem to diverge after each iteration resulting in error propagation. Try increase the precision by using something like https://docs.python.org/2/library/decimal.html. You may be able to get them exact if you use the same precision as matlab. – Ivor Denham-Dyson Nov 24 '19 at 20:40
1

@Matthijs990 @Jake Denham-Dyson if it would be the main problem, there wont be such difference between 33 element in ```matlab``` equal to 0.3032 and 32 element in ```python``` equal to 1.12391865 – macaroni Nov 24 '19 at 20:42

score 0 · Answer 2 · answered Nov 24 '19 at 20:45

In an Octave session:

>> for j=1:10
   disp(j)
   end
 1
 2
 3
 4
 5
 6
 7
 8
 9
 10

In python:

In [20]: for j in range(10):print(j)                                            
0
1
2
3
4
5
6
7
8
9

While MATLAB lets you iterate and still have reasonable speed (it does a jit compile), numpy should use whole-array operations where possible (as in older versions of MATLAB).

I haven't tried to follow the details of your iteration, but expressions like h + n - k-1 and mu[i - n + k] suggest you are trying to simulate 2 dimensional iterations with 1d iteration. If using numpy I expect to see one or more numpy arrays being initialized.

Matlab for loop, different values in Python

2 Answers2