Блог

В приведенном ниже коде я пытаюсь выполнить итерацию по 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<= Tmax. Я не могу понять, в чем моя ошибка, я получаю следующие ошибки:

  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.
  

Глядя на вашу первую итерацию:

 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) , вы применяете этот индекс ко 2-му измерению. Почему? Разве это не должно применяться к измерению с размером NS ?

По numpy умолчанию массивы имеют порядок ‘C’. Последнее измерение является наиболее внутренним. Они могут иметь F (fortran) порядок, но давайте не будем вдаваться в подробности. Думая о 2d массиве как о таблице, мы обычно говорим о строках и столбцах, хотя у них нет формального определения в 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] формат.

Я подозреваю, что другие ваши итерации имеют аналогичные проблемы, но я пока остановлюсь на этом.

Что касается ошибок:

• Первый:

 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] которая является элементом. Вы пытаетесь вставить массив в значение.

• Что касается второй ошибки, я думаю, это может быть связано с третьей. Если бы вы могли предоставить больше информации, я хотел бы попытаться понять, что происходит.

Надеюсь, это помогло.