Короче говоря, мне нужно выполнять векторные операции над 2D-матрицей со значениями самой матрицы для тысяч итераций, но по причинам, которые я объясню ниже, мне нужно сделать это в нескольких разделах, и я хочу знать, лучший способ сделать это, по-прежнему получая максимально возможную производительность и удобочитаемость.
Я решаю уравнение Лапласа, чтобы сгенерировать сетки для компьютерного моделирования аэродинамики.
Для этого, скажем, у меня есть двумерная матрица с именем X формы (M, N), где M и N - количество строк и столбцов соответственно, и мне нужно получить значение каждого внутреннего узла с помощью «координаты» (i, j), на которые влияют его соседи, точки (i+1, j) (i-1, j) (i, j+1) (i, j-1). Возьмем, к примеру, следующее уравнение:
X[i, j] = (X[i+1, j] - X[i-1, j] + X[i, j+1] - X[i, j-1]) / 4
Код выполняется за несколько итераций, порядка сотен тысяч, и на каждой итерации мне нужно пройти по всей матрице, вычисляя каждый внутренний узел. Приведенное выше уравнение гласит, что вычисления выполняются в самой матрице, и что значения X[i-1, j] и X[i, j-1] - это те значения, которые уже вычислены в текущей итерации.
Итак, это фон проблемы под вопросом, теперь к коду я пишу. Как новичок ie я начал с очевидного, не оптимального подхода двух вложенных циклов, одного для строк и одного для столбцов, которые уже находятся внутри некоторого времени l oop (количество итераций):
while current_it < it_max:
for i in range(1, M-1):
for j in range(1, N-1):
X[i, j] = (X[i+1, j] - X[i-1, j] + X[i, j+1] - X[i, j-1]) / 4
Это сработало, и для небольших размеров матрицы она выполнялась за относительно короткое время, примерно через 5 минут, я знаю, что это уже огромное время выполнения, но это не было проблемой. Но мне нужны большие сетки, например, сетка размером 1200 x 400, и в этом случае время выполнения возрастает экспоненциально, и для решения сетки требуются ДНИ, которые более недоступны.
Благодаря этот вопрос , я понял, что могу векторизовать уравнение и избавиться от вложенных циклов for, поэтому теперь мой код выглядит примерно так:
while current_it < it_max:
# replacements of i and j
# j or i --> 1:-1
# (j or i) + 1 --> 2:
# (j or i) - 1 --> :-2
X[1:-1, 1:-1] = (X[2:, 1:-1] - X[:-2, 1:-1] + X[1:-1, 2:] - X[1:-1, :-2]) / 4
Это представляет ОГРОМНОЕ улучшение времени выполнения, сетки что в подходе classi c для его генерации потребуется 3 дня, а теперь, возможно, 5 минут.
У меня сейчас проблема в том, что я потерял способность получать значения (i-1) и (j-1). для текущей итерации, и это заставляет код выполнять гораздо больше итераций, что, как я подозреваю, необходимо.
Мое решение - разделить матрицу на секции и вычислять каждый фрагмент за раз.
while current_it < it_max:
# 1st piece [i, 1 : lim_1]
# 2nd piece [i, lim_1 :]
X[1:-1, 1:lim_1] = (X[2:, 1:lim_1] - X[:-2, 1:lim_1] \
+ X[1:-1, 2:lim_1 + 1] - X[1:-1, :lim_1 - 1]) / 4
X[1:-1, lim_1:-1] = (X[2:, lim_1:-1] - X[:-2, lim_1:-1] \
+ X[1:-1, lim_1 + 1:] - X[1:-1, lim_1 - 1:-2]) / 4
Но я знаю, что копирование - плохая практика, а также строки кода быстро растут, так как мне нужно несколько разделов как в i, так и j указания.
Каков наилучший способ переупорядочения этого окончательного кода, чтобы получить наилучшую производительность и удобочитаемость.