【问题标题】:How to simplify for loops in Python and Matlab如何在 Python 和 Matlab 中简化 for 循环
【发布时间】:2012-02-11 07:11:54
【问题描述】:

我是 C 和 MATLAB 用户。当我(一周前)开始学习 Python 时,我注意到我没有充分发挥 MATLAB 的潜力,尤其是数组运算。我经常使用 for 循环,可能是因为我学过 C 语言。

在之前的一个技巧中,我学会了使用cumsum和其他高效的数组操作,例如:

alpha = [1e-4,1e-3,1e-4,1e-1,1e-2,1e-3,1e-6,1e-3];
zeta = alpha / (dz*dz)
nz = 101
l=[0.3,0.1,0.2,0.1,0.1,0.1,0.2];
wz = cumsum(l*(nz-1));
nl = lenght(l);   

是否可以在 Python (Numpy) 或 MATLAB 中简化以下代码?

      A = zeros(nz,nz);
      i=1;
      for j = 2:wz(i)-1
        A(j,j-1) = zeta(1,1);
        A(j,j) = -2*zeta(1,1);
        A(j,j+1) = zeta(1,1); % layer 1 nodes 
      end

      %cicle to n-layers
      for i=2:nl
          for j=wz(i-1):wz(i-1)
              A(j,j-1) = zeta(1,i-1);
              A(j,j) = -zeta(1,i-1)-zeta(1,i);
              A(j,j+1) = zeta(1,i); 
          end

          for j=wz(i-1)+1:wz(i)
              A(j,j-1) = zeta(1,i);
              A(j,j) = -2*zeta(1,i);
              A(j,j+1) = zeta(1,i);
          end

      end
end

【问题讨论】:

  • 如果你能用一两句话解释代码应该做什么,而不只是给我们简单的代码,对我们来说会更容易......
  • @Hans 此代码属于应用于多层的一维方程热求解器。 alpha 是具有每层扩散率的数组,l 是具有每层高度的数组,wz 是聚合点(离散点)的累积和的数组,A 是“状态矩阵”。计算状态矩阵后,我将实现 ode 求解器。
  • @marco:如果您没有意识到,MATLAB 内置了各种 ODE solvers,而对于 Python,您可以在 Scipy 中找到 ODE 求解器。
  • 我将使用内置的 ODE 求解器。但矩阵 A 是求解器的输入。这段代码是正确的(解决了问题),但效率不高。

标签: python arrays matlab for-loop numpy


【解决方案1】:

为了简化你的循环,你可以使用函数spdiags

http://www.mathworks.fr/help/techdoc/ref/spdiags.html

例如你的第一个循环可以写成:

A=full(spdiags(repmat([zeta(1,1),-2*zeta(1,1),zeta(1,1)],wz(i),1),[-1 0 1],wz(i),wz(i)))

【讨论】:

  • 感谢您的参考,但这并不能解决我的问题 =) 我正在尝试找到一种方法来简化问题,换句话说,将迭代问题转换为向量问题
  • 我编辑了答案,向您展示如何在您的情况下使用spdiags
  • 没看懂,把代码粘贴到matlab中,返回???索引超出矩阵维度。 ==> spdiags 在 114 处出错 a((len(k)+1):len(k+1),:) = [i i+d(k) B(i+(m>=n)*d(k ),k)];
  • 有效,但结果不同。 “我的矩阵”的对角线(有 3 个值)为 100,-200,100。 “你的”总是 100。你的矩阵有 30,30 大小,这是正常的,因为它只是循环的第一部分。
  • 我在最后一条评论中忘了提到@Oli =)
【解决方案2】:

在有机会在我的机器上与您的机器并排运行后,我修改了下面的代码。还有几个问题(假设 A 在最终循环中变大了吗?,什么是 dz?)。你在运行它之前遇到的问题是我忘记了 idx_matrix 必须是合乎逻辑的。

dz=0.1;
alpha = [1e-4,1e-3,1e-4,1e-1,1e-2,1e-3,1e-6,1e-3];
zeta = alpha / (dz*dz);
nz = 101;
l=[0.3,0.1,0.2,0.1,0.1,0.1,0.2];
wz = cumsum(l*(nz-1));
nl = length(l);

A = zeros(nz);
i=1;

%replaces 1st loop
j_start = 2;
j_end = wz(i)-1;

idx_matrix = false(size(A));
idx_matrix(j_start:j_end,j_start:j_end) = eye(j_end-j_start+1);
A(idx_matrix) = -2*zeta(1,1);

idx_matrix(idx_matrix) = false;
idx_matrix(j_start:j_end,j_start-1:j_end-1) = eye(j_end-j_start+1);
A(idx_matrix) = zeta(1,1);

idx_matrix(idx_matrix) = false;
idx_matrix(j_start:j_end,j_start+1:j_end+1) = eye(j_end-j_start+1);
A(idx_matrix) = zeta(1,1);

%cicle to n-layers
for i=2:nl

    %replaces 3rd loop
    j_start = wz(i-1);
    A(j_start,j_start) = -zeta(1,i-1)-zeta(1,i);
    A(j_start,j_start-1) = zeta(1,i-1);
    A(j_start,j_start+1) = zeta(1,i);

    %replaces 4th loop
    j_start = wz(i-1)+1;
    j_end = min(wz(i),size(A,2)-1);
    idx_matrix = false(size(A));
    idx_matrix(j_start:j_end,j_start:j_end) = eye(j_end-j_start+1);
    A(idx_matrix) = -2*zeta(1,i);

    idx_matrix(idx_matrix) = false;
    idx_matrix(j_start:j_end,j_start-1:j_end-1) = eye(j_end-j_start+1);
    A(idx_matrix) = zeta(1,i);

    idx_matrix(idx_matrix) = false;
    idx_matrix(j_start:j_end,j_start+1:j_end+1) = eye(j_end-j_start+1);
    A(idx_matrix) = zeta(1,i);

end

【讨论】:

  • 感谢您的回答。你的代码有错误???下标索引必须是实数正整数或逻辑数。 A(idx_matrix) = -2*zeta(1,1);代码太早了,我不知道你想做什么。是的,第三个循环,可以是一个作业,我没有找到更好的方法来这条指令。 (这代表了层之间的边界,如果你把我的代码描述读给@Hans)
  • AHh,“you're A matrix”是一个数组 (A=zeros(nz)) 是不是错误?
  • A = zeros(nz) 与 A = zeros(nz,nz) 相同
  • 您是否打算让 A 在此代码继续运行时变得更大?最后的 for 循环向 A 添加一列,因为 wz(i) == size(A,2) 和 A(j,j+1) 索引超出范围,因此添加了另一列。另外,你从来没有定义过 dz。
  • 你是对的。你的代码有效。然而,在最后一个循环中,“zeta”因子不会对每个系数进行叠加。可以保持通用吗?如果 l 的大小在变化?
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 2021-07-17
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2021-10-06
  • 2012-06-23
  • 1970-01-01
相关资源
最近更新 更多