使用 while 和 for 语句迭代 2d numpy 数组

Question

在下面提供的代码中，我试图迭代 2D numpy 数组 [i][k] 最初它是用比我祖父还老的 Fortran 77 编写的代码。我正在尝试使其适应 python。 （对于感兴趣的人：它是一个简单的液压瞬态事件求解器） 请记住，所有变量都是在我的代码中引入的，我不会在此处粘贴。

  H = np.zeros((NS,50))
  Q = np.zeros((NS,50))

我在这里分配第一行值：

    for i in range(NS):
       H[0][i] = HR-i*R*Q0**2
       Q[0][i] = Q0

 CVP = .5*Q0**2/H[N]
 T = 0
 k = 0
 TAU = 1
 #Interior points:
 HP = np.zeros((NS,50))
 QP = np.zeros((NS,50))

 while T<=Tmax:
     T += dt
     k += 1
     for i in range(1,N):
         CP = H[k][i-1]+Q[k][i-1]*(B-R*abs(Q[k][i-1]))
         CM = H[k][i+1]-Q[k][i+1]*(B-R*abs(Q[k][i+1]))
         HP[k][i-1] = 0.5*(CP+CM)
         QP[k][i-1] = (HP[k][i-1]-CM)/B
 #Boundary Conditions:
 HP[k][0] = HR
 QP[k][0] = Q[k][1]+(HP[k][0]-H[k][1]-R*Q[k][1]*abs(Q[k][1]))/B
     if T == Tc:
         TAU = 0
         CV = 0
     else:
         TAU = (1.-T/Tc)**Em
         CV = CVP*TAU**2
     CP = H[k][N-1]+Q[k][N-1]*(B-R*abs(Q[k][N-1]))
     QP[k][N] = -CV*B+np.sqrt(CV**2*(B**2)+2*CV*CP)
     HP[k][N] = CP-B*QP[k][N]
     for i in range(NS):
         H[k][i] = HP[k][i]
         Q[k][i] = QP[k][i]

记住i 用于行，k 用于列 我期望的是，对于所有 k 列数，应计算值，直到满足 T

 RuntimeWarning: divide by zero encountered in true_divide
 CVP = .5*Q0**2/H[N]

 RuntimeWarning: invalid value encountered in multiply
 QP[N][k] = -CV*B+np.sqrt(CV**2*(B**2)+2*CV*CP)

 QP[N][k] = -CV*B+np.sqrt(CV**2*(B**2)+2*CV*CP)
 ValueError: setting an array element with a sequence.

Answer 1

查看您的第一次迭代：

H = np.zeros((NS,50))
Q = np.zeros((NS,50))
for i in range(NS):
   H[0][i] = HR-i*R*Q0**2
   Q[0][i] = Q0

H 的形状是 (NS,50)，但是当您迭代 range(NS) 时，您将该索引应用于第二维。为什么？它不应该适用于大小为NS的维度吗？

numpy 中的数组默认为“C”顺序。最后一个维度是最里面的。他们可以有F (fortran) 订单，但我们不要去那里。将二维数组视为一个表，我们通常谈论行和列，尽管它们在 numpy 中没有正式定义。

假设您要将第一列设置为这些值：

for i in range(NS):
    H[i, 0] = HR - i*R*Q0**2
    Q[i, 0] = Q0

但是我们可以一次完成整个行或列的赋值。我相信新版本的 Fortran 也有这些“全数组”功能。

Q[:, 0] = Q0  
H[:, 0] = HR - np.arange(NS) * R * Q0**2

翻译成Python 时要注意一点。索引从 0 开始；范围和np.arange(...) 也是如此。

H[0][i] 在功能上与H[0,i] 相同。但是在使用切片时，您必须使用H[:,i] 格式。

我怀疑你的其他迭代也有类似的问题，但我现在就到此为止。

Answer 2

关于错误：

• 第一个：

RuntimeWarning: divide by zero encountered in true_divide
 CVP = .5*Q0**2/H[N]

您将 H 初始化为零，因此它抱怨被零除是正常的。也许你应该添加一个条件。

• 第三次：

 QP[N][k] = -CV*B+np.sqrt(CV**2*(B**2)+2*CV*CP)
 ValueError: setting an array element with a sequence.

您定义CVP = .5*Q0**2/H[N]，然后定义CV = CVP*TAU**2，这是一个序列。然后您尝试将其派生形式分配给QP[N][K]，这是一个元素。您正在尝试将数组插入到值中。

• 对于第二个错误，我认为它可能与第三个有关。如果您能提供更多信息，我想尝试了解会发生什么。

希望这有帮助。