Алгоритм продольно-поперечной прогонки

Основные авторы описания: Дружкина А.И., Лиходед Н.А.

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Продольно–поперечная схема, которая также носит название метода переменных направлений (ADI – alternate directions implicit method), получила широкое применение для решения многомерных задач, приводящих к уравнениям в частных производных параболического типа (уравнение диффузии, уравнение теплопроводности). Эта схема была предложена в 1955 году Писменом и Рэкфордом.

1.2 Математическое описание алгоритма

Рассмотрим в области

[math]\begin{align} D_{T}=G\times [0\leq t\leq T], \\ G=[0\leq x\leq l_{x}]\times[0\leq y\leq l_{y}]\end{align}[/math]

двумерное нелинейное уравнение теплопроводности

[math]\frac{\partial u}{\partial t} = \frac{\partial }{\partial x_{1}}\left ( k_{1}(x_{1}, x_{2}, t)\frac{\partial }{\partial x_{1}} \right ) + \frac{\partial }{\partial x_{2}}\left ( k_{1}(x_{1}, x_{2}, t)\frac{\partial }{\partial x_{2}} \right ) + f\left ( x_{1}, x_{2}, t \right ) [/math]

с начальными

[math]T(x,y,0) = T_{0}(x,y), (x,y)\in G,[/math]

и граничными условиями

[math]T(x,y,t)|_{S_{T}} = \mu (x,y,t), (x,y,t)\in S_{T}\equiv \partial D_{T}[/math]

Введем в области [math]G[/math] сетку узлов [math]\bar{\omega_{h}} = \left \{ x^{i}=i\cdot h_{x}, i=0,1,...,N_{x}, N_{x}\cdot h_{x}=l_{x},y^{j}=j\cdot h_{y}, j=0,1,...,N_{y}, N_{y}\cdot h_{y}=l_{y} \right \} [/math], а на отрезке [math][0\leq t\leq T][/math] сетку узлов [math]\omega_{\tau } = \left \{ t_{j}=n\tau, n=0,1,...,k, k\tau=T \right \}. [/math]. Запишем разностную схему для заданной задачи ([math]k_{1} = k_{2} = \lambda [/math]):

[math]\left\{\begin{matrix} \frac{y_{i,j}^{n+\frac{1}{2}} - y_{i,j}^{n}}{\tau} = \frac{1}{2}\left ( \frac{\lambda_{i+\frac{1}{2}, j}(y_{i+1,j}^{n+\frac{1}{2}} - y_{i,j}^{n+\frac{1}{2}}) - \lambda_{i-\frac{1}{2}, j}(y_{i,j}^{n+\frac{1}{2}} - y_{i-1,j}^{n+\frac{1}{2}})}{h_{x}^{2}} \right ) + \frac{1}{2}\left ( \frac{\lambda_{i, j+\frac{1}{2}}(y_{i,j+1}^{n} - y_{i,j}^{n}) - \lambda_{i, j-\frac{1}{2}}(y_{i,j}^{n} - y_{i,j-1}^{n})}{h_{y}^{2}} \right ), \\ \frac{y_{i,j}^{n+1} - y_{i,j}^{n+\frac{1}{2}}}{\tau} = \frac{1}{2}\left ( \frac{\lambda_{i+\frac{1}{2}, j}(y_{i+1,j}^{n+\frac{1}{2}} - y_{i,j}^{n+\frac{1}{2}}) - \lambda_{i-\frac{1}{2}, j}(y_{i,j}^{n+\frac{1}{2}} - y_{i-1,j}^{n+\frac{1}{2}})}{h_{x}^{2}} \right ) + \frac{1}{2}\left ( \frac{\lambda_{i, j+\frac{1}{2}}(y_{i,j+1}^{n+1} - y_{i,j}^{n+1}) - \lambda_{i, j-\frac{1}{2}}(y_{i,j}^{n+1} - y_{i,j-1}^{n+1})}{h_{y}^{2}} \right ) \end{matrix}\right.[/math]

где [math]\lambda _{i\pm \frac{1}{2}, j} = \frac{\lambda _{i\pm 1, j} + \lambda _{i, j}}{2}, [/math] [math]\lambda _{i, j\pm \frac{1}{2}} = \frac{\lambda _{i, j\pm 1} + \lambda _{i, j}}{2}. [/math]

Зафиксировав в первом из уравнений системы [math]j[/math], получим систему уравнений относительно значений [math]y_{i,j}^{n+\frac{1}{2}}[/math], где [math]j = 1,...,N_{x}-1[/math], состоящую из ([math]N_{x}-1[/math]) линейного уравнения, которую можно решить методом прогонки. В целом систему на каждом половинном временном слое можно представить как ([math]N_{x}-1[/math]) независимую задачу (для каждого фиксированного [math]j[/math]), решаемую методом прогонки. Аналогично решение второго из уравнений системы на каждом слое [math]t_{n+1}[/math] представляет собой решение ([math]N_{x}-1[/math]) независимой задачи при фиксированном [math]i[/math]. Каждая из указанных задач является системой линейных уравнений относительно значений сеточной функции по неявному направлению и решается методом прогонки. Сеточная функция , является приближенным решением задачи. По каждому из неявных направлений разностная схема является линейной и может быть записана в следующем виде:

[math]A_{i,j}y_{i-1,j}^{n+\frac{1}{2}} + C_{i,j}y_{i,j}^{n+\frac{1}{2}} + B_{i,j}y_{i+1,j}^{n+\frac{1}{2}} = F_{i,j}^{n}[/math]

[math]A_{i,j} = \frac{\lambda_{i-\frac{1}{2},j}}{-2h_{y}^{2}}, B_{i,j} = \frac{\lambda_{i+\frac{1}{2},j}}{-2h_{y}^{2}}, C_{i,j} = \frac{1}{\tau} - A_{i,j} - B_{i,j},[/math]

[math]F_{i,j}^{n+\frac{1}{2}} = \frac{y_{i,j}^{n+\frac{1}{2}}}{\tau} + \frac{\lambda _{i+\frac{1}{2}, j}(y_{i+1, j}^{n+\frac{1}{2}} - y_{i, j}^{n+\frac{1}{2}}) - \lambda _{i-\frac{1}{2}, j}(y_{i,j}^{n+\frac{1}{2}} - y_{i-1, j}^{n+\frac{1}{2}})}{h_{y}^{2}}[/math]

и соответственно

[math]A_{i,j}y_{i,j-1}^{n+1} + C_{i,j}y_{i,j}^{n+1} + B_{i,j}y_{i,j+1}^{n+1} = F_{i,j}^{n+\frac{1}{2}}. [/math]

Для решения уравнений воспользуемся формулами прогонки. Значения прогоночных коэффициентов находятся по рекуррентным формулам:

где [math]k[/math] – индекс неявного направления. Из граничных условий при [math]k=0[/math] и [math]k=N[/math] определяются значения прогоночных коэффициентов. При этом [math]\alpha _{0}[/math] и [math]\alpha _{N}[/math] равны нулю, а значения [math]\beta _{0}[/math] и [math]\beta _{N}[/math] определяются из соответствующих краевых условий.

1.3 Вычислительное ядро алгоритма

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

1.4 Макроструктура алгоритма

Алгоритм представляет собой совокупность продольной и поперечной прогонки, а также прямого и обратного хода.

1.5 Схема реализации последовательного алгоритма

Последовательность исполнения метода следующая: Осуществляется прогонка вдоль строк, как это изображено на рисунке 1.1 [4]:

Рисунок 1.1– Прогонка вдоль строк

Затем осуществляется прогонка вдоль столбцов, как это представлено на рисунке 1.2 [4]:

Рисунок 1.2- Прогонка вдоль столбцов

1.6 Последовательная сложность алгоритма

Таким образом, при классификации по последовательной сложности продольно–поперечная прогонка относится к алгоритмам с линейной сложностью. При переходе от слоя [math]j[/math] к слою [math]j+1[/math] требуется [math]O(\frac{1}{h^{2}})[/math] арифметических действий. Чтобы найти [math]y^{j_{0}}[/math] при [math]t_{0} = j_{0}\tau[/math] по начальным данным требуется, очевидно, [math]O(\frac{1}{h^{2}})j_{0} = O(\frac{1}{\tau h^{2}})[/math] операций, то есть число операций пропорционально числу используемых узлов пространственно–временной сетки [math]w_{h\tau}[/math]. Наряду с основными значениями [math]u_{ij}^{k}[/math] и [math]u_{i,j}^{k+1}[/math] вводится значение на промежуточном слое – [math]u_{ij}^{k+\frac{1}{2}}[/math], что фактически является значением [math]u[/math] при [math]t = t_{k+\frac{1}{2}}=t+\frac{\tau}{2}[/math]. Благодаря этому, переход на следующий слой осуществляется в два шага.

1.8 Описание ресурса параллелизма алгоритма

Алгоритм метода переменных направлений обладает естественным параллелизмом: можно организовать независимые вычислительные процессы на каждом временном слое.

1.9 Входные и выходные данные алгоритма

Входные данные: матрица [math]y[/math] (элементы [math]y_{i,j}^{1}, i = 0,...,N_{x}, j = 0,...,N_{y}[/math]).

Конечные значения [math]y_{(2)}(i_{1}, i_{2})[/math] являются приближенным решением задачи

Выходные данные: обновленная матрица [math]y[/math] (элементы [math]y_{i,j}^{n+1}, i = 0,...,N_{x}, j = 0,...,N_{y}[/math]).

Объём выходных данных: [math](N_{x} + 1) * (N_{y} + 1)[/math]

1.10 Свойства алгоритма

Продольно–поперечная схема является одной из первых экономичных схем. Она сочетает в себе лучшее качество явной схемы – экономичность и неявной – устойчивость. Основной идеей экономичных разностных схем является сведение многомерной задачи к цепочке одномерных задач. Продольно–поперечная схема равномерно и безусловно устойчива по начальным данным, так как при переходе с одного целого слоя на следующий целый слой ошибки начальных данных не нарастают. При переходе с целого слоя на целый погрешность локальной аппроксимации на равномерных сетках есть [math]O(\tau^{2} + h_{x}^{2} + h_{y}^{2})[/math] т.е. продольно–поперечная схема имеет второй порядок аппроксимации по всем переменным.

2 Программная реализация алгоритма

2.1 Особенности реализации последовательного алгоритма

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

2.3 Возможные способы и особенности параллельной реализации алгоритма

Вычислительные процессы реализаций, основанных на естественном па-раллелизме, на каждом временном слое могут выполняться независимо, но при переходе к новому слою требуют групповых коммуникаций “каждый с каждым”. С ростом числа используемых процессов это приводит к большим накладным расходам. Локальность алгоритма – это вычислительное свойство, отражающее степень использования при реализации алгоритма памяти с быстрым доступом. При многопроцессорной обработке памятью с быстрым доступом считается локальная память процессора. В этой работе вместо части естественного параллелизма предлагается использовать конвейерный параллелизм. Такой параллелизм приводит к разгону и торможению вычислительного конвейера на каждом временном слое, но позволяет построить параллельный алгоритм с улучшенной локальностью (с существенно меньшим объемом коммуникационных операций).

2.3.1 Вычислительный алгоритм продольно-поперечной прогонки

Основную вычислительную часть алгоритма, реализующего двухшаговую разностную схему продольно–поперечной прогонки, т.е. метода переменных направлений, можно представить в виде следующего псевдокода (циклы, итерации которых заведомо можно выполнять независимо, запишем как dopar):

do  j =1, j0  
     dopar  i2= 1, N2 –1 
          do  i1= 1, N1 –1
S1:            
S2:            
          enddo
          do  i1= 1, N1 –1
S3:           
          enddo
     enddopar
     dopar  i1= 1, N1 –1
          do  i2= 1, N2 –1
S4:            
S5:            
          enddo
          do  i2= 1, N2 –1
S6:          
          enddo
     enddopar
enddo(j)

Здесь (i) и (i) – коэффициенты прогонки, возникающие при решении промежуточных систем линейных алгебраических уравнений. Конечные значения являются приближенным решением задачи. Алгоритм метода переменных направлений обладает естественным параллелизмом: можно организовать независимые вычислительные процессы на каждом временном слое. Так как число вычислительных ядер обычно гораздо меньше N1 и N2, то при реализации на параллельных компьютерах с распределенной памятью итерации циклов, отмеченных как dopar, следует сгруппировать и получить макрооперации, называемые зернами вычислений или тайлами. Процесс получения макроопераций–тайлов называется тайлингом.

2.7 Существующие реализации алгоритма

В языке Matlab существует реализация алгоритма продольно-поперечной прогонки. Также метод продольно–поперечной прогонки реализован в библиотеке для решения дифференциальных уравнений языка Python.