Matlab for循环，Python中的不同值

【问题标题】：Matlab for loop, different values in PythonMatlab for循环，Python中的不同值
【发布时间】：2020-03-20 03:35:07
【问题描述】：

我在 matlab 中有代码，我正在尝试将其“重写”成 Python。问题在于 for 循环 的值不同。不幸的是，我找不到解决方案。 问题在于每个数组中的值：nu theta mu M 来自 i 元素（matlab 中的 32 和 python 中的 31）。

对于两个代码n=30

matlab 中的 32 个元素 = python 中的 31 个元素

在 Python 中，DTOR 或 RTOD 是弧度的度数，反之亦然。

这是整个python code

----------------------------------Matlab： em>-------------------------


%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

以及Matlab“theta”值的前几行：

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： em>-------------------------

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

以及Python“theta”值的前几行：

[ 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

【问题讨论】：

  • 序列在第 32 位之后开始不同。我猜其中一个分支出了点问题。我建议附上一个标签，表明哪个分支产生了特定的价值
  • @tstanisl 我用 cmets % in matlab 和 # in python 编辑了它
  • @hpaulj range(n-1) 表示 python 中从 0 到 28，因此它迭代 29 次。在 matlab for j= 1:n-1 是从 1 到 29 所以它也迭代 29 次
  • @JohanC 结果没有任何改变。但是for j in range(1,n) 和for j in range(n-1) 都迭代了 29 次
  • 一个较小的速度点，对于整个数组np.sin 是正确的，但对于单个值，math.sin 更快。其他三角计算也一样。

标签： python matlab numpy for-loop


【解决方案1】：

在 Octave 会话中：

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

在python中：

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

虽然 MATLAB 允许您进行迭代并且仍然具有合理的速度（它进行 jit 编译），但numpy 应该尽可能使用整个数组操作（就像在旧版本的 MATLAB 中一样）。

我没有尝试关注您的迭代细节，但像h + n - k-1 和mu[i - n + k] 这样的表达式建议您尝试使用一维迭代来模拟二维迭代。如果使用numpy，我希望看到一个或多个numpy 数组被初始化。

【讨论】：

    【解决方案2】：

    python 计算的一个问题是它们经常出错 那么可能是彼此不相等的原因之一 您可以在此处找到有关它的更多信息： Is floating point math broken?

    【讨论】：

    • 这是一个答案吗？
    • 但直到matlab 中的 32 个元素（python 中的 31 个元素）一切正常。我猜i 元素编号有问题，因为循环从这里开始
    • @tstanisl 这是一个建议
    • 我的第一个猜测确实是由于浮点运算导致的舍入误差，因为每次迭代后结果似乎都会发散，从而导致错误传播。尝试使用docs.python.org/2/library/decimal.html 之类的方法来提高精度。如果您使用与 matlab 相同的精度，您可能能够得到准确的结果。
    • @Matthijs990 @Jake Denham-Dyson 如果这将是主要问题，那么matlab 中的 33 个元素等于 0.3032 和 python 中的 32 个元素等于 1.12391865 之间不会有这种差异跨度>
    猜你喜欢
    • 2021-10-06
    • 1970-01-01
    • 2020-10-14
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2013-06-19
    • 1970-01-01
    相关资源
    最近更新 更多
    热门标签
    Java Python linux javascript Mysql C# Docker 算法 前端 SpringBoot Redis Vue spring 设计模式 .net core .net kubernetes c++ 数据库 数据结构 大数据 js 机器学习 微服务 Android Go 程序员 面试 JVM ASP.net core 云原生 人工智能 后端 PHP git CSS golang k8s Nginx Django mybatis 深度学习 多线程 React 架构 devops 爬虫 云计算 Spring Boot LeetCode