Итерации над двумерными массивами с операторами while и for - PullRequest
0 голосов
/ 01 апреля 2019

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

 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.

Ответы [ 2 ]

1 голос
/ 01 апреля 2019

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

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 (фортран), но давайте не будем туда идти. Думая о 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].

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

0 голосов
/ 01 апреля 2019

Относительно ошибок:

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

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

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

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...