Блог

Я пытаюсь создать матрицу из трех индексов, которая содержит 1 значение (V) для каждого узла числовой пространственной сетки (xyz) (проблема реального мира: электростатический потенциал, создаваемый конечными параллельными пластинами в точке пространства). Эта матрица изначально должна быть заполнена нулями, за исключением некоторых определенных точек (где находятся таблички и ограничения пространства), а затем итеративно обновлять значение в каждом узле в соответствии со следующим методом трафарета из 7 точек (индексы jk и l координат x y z соответственно):

V [j, k, l] = (V [j 1, k, l] V [j-1, k, l] V [j, k 1, l] V [j, k-1, l] V [j, k,l 1] V[j, k, l-1]) / 6 (т. е. Замените значение узла средним значением других 6 соседних узлов)

Я пробовал np.zeros и np.meshgrid, но я думаю, что, возможно, у меня просто серьезный концептуальный и базовый пробел в отношении массивов, поскольку, похоже, ничто не делает то, что я хочу. Любая ориентация была бы очень признательна и извините, если я неправильно объяснил себя. Вот некоторый код, который я пробовал:

 V1 = 10
V2 = -5
Mx = 101
My = 151
Mz = 301

V = np.zeros([Mx, My, Mz]).astype(int)
V[46, 51:101, 101:201] = V1   #the values of these nodes should stay fixed throughout iteration
V[56, 51:101, 101:201] = V2   #the values of these nodes should stay fixed throughout iteration
V[1,:,:] =V[100,:,:] =V[:,1,:] =V[:,150,:] =V[:,:,1] =V[:,:,300] = 0     #the values of these nodes should stay fixed throughout iteration

for j  in V:
    for k in j:
        for l in k:
            V[j, k, l] = (V[j 1, k, l]   V[j-1, k, l]   V[j, k 1, l]   V[j, k-1, l]   V[j, k, l 1]  V[j, k, l-1])/6
 

(Обновление после справки от пользователя kcw78)

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

«Будет выполнено столько этих циклов, сколько необходимо, чтобы ошибка упала ниже определенного заданного допуска, rtol. И какова хорошая мера ошибки здесь? Мы будем использовать максимальное значение локального остатка, определяемое как (абсолютное значение) разницы между потенциальным значением в центральном узле и средним арифметическим других значений в шаблоне. В качестве дополнительной меры предосторожности мы также сравним ошибки любых двух последовательных циклов и остановим ослабление, если они станут равными. Лучшее решение больше невозможно «.

Теперь пытаюсь выполнить приведенный ниже код, но не уверен, попал ли он в бесконечный цикл while или просто занимает много времени, так как я должен остановить его через 20 минут без получения какого-либо вывода (также не уверен, может быть, мне следует использовать .all() вместо .any()):

 import numpy as np

V1 = 10
V2 = -5
Mx = 101
My = 151
Mz = 301
rtol = 10**-2

V1_set = { (46,k,l) for k in range(51,101,1) for l in range(101,201,1) }
V2_set = { (56,k,l) for k in range(51,101,1) for l in range(101,201,1) }

V = np.zeros((Mx, My, Mz))
Vnew = np.copy(V)
V[46, 51:101, 101:201] = V1   
V[56, 51:101, 101:201] = V2   
V[1,:,:] =V[100,:,:] =V[:,1,:] =V[:,150,:] =V[:,:,1] =V[:,:,300] = 0



check_set = set().union(V1_set,V2_set)

error = np.zeros((Mx, My, Mz))
errornew = np.zeros((Mx, My, Mz))

while float(errornew.any()) < rtol or error.any() != errornew.any():
 V = Vnew
 error = errornew
 for j in range(1,V.shape[0]-1):
    for k in range(1,V.shape[1]-1):
        for l in range(1,V.shape[2]-1):
            if (j,k,l) not in check_set:
                Vnew[j, k, l] = (V[j 1, k, l]   V[j-1, k, l]   V[j, k 1, l]   V[j, k-1, l]   V[j, k, l 1]  V[j, k, l-1])/6
                errornew[j, k, l] = abs(Vnew[j, k, l]-V[j, k, l])
 
 

Если я понимаю ваш вопрос, вам потребуется 2 изменения:

1. Сначала вам понадобятся дополнительные переменные, чтобы проверить позиции, которые фиксируются на протяжении всей итерации. Для этого я добавил наборы с (j,k,l) кортежами. Чтобы вы могли следовать моей логике, я изначально создал 3 набора; по 1 для каждого из этих индексов: 1) фиксированный V1 ( V1_set ), 2) фиксированный V2 ( V2_set ) и 3) boundary ( zero_set ), затем объединил все 3 набора в один набор (вызывается check_set ). Вы могли бы начать с одного набора и обновлять его по мере добавления. Примечание: в вашем коде есть V[1,:,:] = 0 , но я думаю, что вы действительно хотите V[0,:,:] = 0 . Дайте мне знать, если я неправильно это истолковал.
2. Во-вторых, вам нужно выполнить цикл по длине оси в каждом направлении (атрибуты есть V.shape[0], V.shape[1], V.shape[2] ). Внутри цикла я проверяю каждую (i,j,k) из check_set них и вычисляю новое V1[j, k, l] значение только в том случае, если его НЕТ в наборе.

Смотрите код ниже:

 V1 = 10
V2 = -5
Mx = 101
My = 151
Mz = 301

V1_set = { (46,k,l) for k in range(51,101,1) for l in range(101,201,1) }
V2_set = { (56,k,l) for k in range(51,101,1) for l in range(101,201,1) }

zero_set = set()
zero_set.update( { (0,k,l) for k in range(My) for l in range(Mz) } )
zero_set.update( { (100,k,l) for k in range(My) for l in range(Mz) } )
zero_set.update( { (j,0,l) for j in range(Mx) for l in range(Mz) } )
zero_set.update( { (j,150,l) for j in range(Mx) for l in range(Mz) } )
zero_set.update( { (j,k,0) for j in range(Mx) for k in range(My) } )
zero_set.update( { (j,k,300) for j in range(Mx) for k in range(My) } )

check_set = set().union(V1_set,V2_set,zero_set)

V = np.zeros((Mx, My, Mz)).astype(int)
V[46, 51:101, 101:201] = V1   #the values of these nodes should stay fixed throughout iteration
V[56, 51:101, 101:201] = V2   #the values of these nodes should stay fixed throughout iteration
V[1,:,:] =V[100,:,:] =V[:,1,:] =V[:,150,:] =V[:,:,1] =V[:,:,300] = 0     #the values of these nodes should stay fixed throughout iteration

for j in range(V.shape[0]):
    for k in range(V.shape[1]):
        for l in range(V.shape[2]):
            if (j,k,l) not in check_set:
                V[j, k, l] = (V[j 1, k, l]   V[j-1, k, l]   V[j, k 1, l]   V[j, k-1, l]   V[j, k, l 1]  V[j, k, l-1])/6
 

После публикации приведенного выше решения мне пришло в голову, что диапазоны, используемые в zero_set , предназначены для того, чтобы избежать первого / последнего индексов (границы массива). Если это так, то в этом нет никакой необходимости zero_set . Вы можете справиться с этим, изменив аргументы диапазона, как показано ниже:

 check_set = set().union(V1_set,V2_set)
for j in range(1,V.shape[0]-1):
    for k in range(1,V.shape[1]-1):
        for l in range(1,V.shape[2]-1):
            if (j,k,l) not in check_set:
                V[j, k, l] = (V[j 1, k, l]   V[j-1, k, l]   V[j, k 1, l]   V[j, k-1, l]   V[j, k, l 1]  V[j, k, l-1])/6
 

Дополнительные замечания для рассмотрения:

• Я заметил , что вы создали array V с .astype(int) помощью . Вы уверены, что это то, чего вы хотите (а не поплавки)? Как правило, ваши вычисления не будут возвращать целочисленные значения.
• В зависимости от того, как написан ваш код, вы меняете значения V[j,k,l] по ходу дела. Таким образом, вы используете обновленные значения V[j,k,l] для j,k,l меньших, чем текущее j,k,l , и предыдущие V[j,k,l] значения для j,k,l больших, чем текущее j,k,l .
• Наконец, я предполагаю, что вы собираетесь повторять это вычисление до тех пор, пока изменение между 2 циклами не станет «приемлемо малым». Если это так, вам нужно иметь 2 копии массива («старый» и «новый»), чтобы получить разницу. Позаботьтесь об использовании .copy() при копировании для создания нового / другого объекта np.array.

Это обновленный ответ, основанный на новой информации и коде, добавленных в начальный пост. У вас есть как минимум 1 проблема с вашей логикой. if (j,k,l) not in check_set: Блок пропускает значения (j, k, l), которые вы хотите сохранить постоянными. В результате вы не выполняете вычисления Vnew в этих точках. Это вызовет проблемы с вычислением изменений с каждой итерацией (и даст неверный результат). Кроме того, я думаю, вам нужно V = Vnew.copy() . В противном случае V и Vnew ссылайтесь на тот же объект.

Вот мой простой подход к итерации с жестко запрограммированной погрешностью.

 check_set = set().union(V1_set,V2_set)
Vi = V.copy()
Vn = np.zeros((Mx, My, Mz))
diff = max(abs(V1), abs(V2))
i = 1
print('Start Cycle#',i,'; diff =',diff)
while diff > 0.25:
    for j in range(1,V.shape[0]-1):
        for k in range(1,V.shape[1]-1):
            for l in range(1,V.shape[2]-1):
                if (j,k,l) in check_set:
                    Vn[j, k, l] = Vi[j, k, l]
                else:
                    Vn[j, k, l] = (Vi[j 1, k, l]   Vi[j-1, k, l]   Vi[j, k 1, l]   Vi[j, k-1, l]   Vi[j, k, l 1]  Vi[j, k, l-1])/6      
      
    diff = max(abs(np.amax(Vn-Vi)), abs(np.amin(Vn-Vi)))
    print('Cycle#',i,'completed; diff =',diff)
    i  = 1
    Vi = Vn.copy()
 

Эта реализация будет «сходиться» за 10 итераций. Однако это только проверяет, что ошибка между двумя последовательными циклами меньше жестко заданного допуска (аналогично второй части желаемой проверки ошибок).
Я НЕ реализовал первую проверку ошибок: «используйте максимальное значение локального остатка, определяемое как (абсолютное значение) разницы между потенциальным значением в центральном узле и средним арифметическим других значений в шаблоне». Я не уверен на 100% в намерении. Является ли трафарет 6 точками вокруг [j,k,l] ? Если это так, я думаю, вам понадобится аналогичный расчет ПОСЛЕ того, как вы вычислите новые значения Vn, что-то вроде этого:

 error[j, k, l] = abs(Vn[j, k, l] - (Vn[j 1, k, l]   Vn[j-1, k, l]   Vn[j, k 1, l]   Vn[j, k-1, l]   Vn[j, k, l 1]  Vn[j, k, l-1])/6 )