Параллельные методы решения дифференциальных уравнений в частных производных

Рассматриваются вопросы организации параллельных вычислений для решения задач, в которых при математическом моделировании используются дифференциальные уравнения в частных производных. Для численного решения подобных задач обычно применяется метод конечных разностей (метод сеток), обладающий высокой вычислительной трудоемкостью. В лекции последовательно разбираются возможные способы распараллеливания сеточных методов на многопроцессорных вычислительных системах с общей и распределенной памятью. При этом большое внимание уделяется проблемам, возникающим при организации параллельных вычислений, анализу причин появления таких проблем и нахождению путей их преодоления. Для наглядной демонстрации излагаемого материала в качестве учебного примера рассматривается проблема численного решения задачи Дирихле для уравнения Пуассона
Дифференциальные уравнения в частных производных представляют собой широко применяемый математический аппарат при разработке моделей в самых разных областях науки и техники. К сожалению, явное решение этих уравнений в аналитическом виде оказывается возможным только в частных простых случаях, и, как результат, возможность анализа математических моделей, построенных на основе дифференциальных уравнений, обеспечивается при помощи приближенных численных методов решения.
Объем выполняемых при этом вычислений обычно является значительным, и использование высокопроизводительных вычислительных систем традиционно для данной области вычислительной математики.
Проблематика численного решения дифференциальных уравнений в частных производных является областью интенсивных исследований.
Рассмотрим в качестве учебного примера проблему численного решения задачи Дирихле для уравнения Пуассона, которая определяется как задача нахождения функции u=u(x,y), удовлетворяющей в области определения D уравнению

и принимающей значения g(x,y) на границе D0 области D (f и g являются функциями, задаваемыми при постановке задачи). Подобная модель может применяться для описания установившегося течения жидкости, стационарных тепловых полей, процессов теплопередачи с внутренними источниками тепла и деформации упругих пластин. Данный пример часто используется в качестве учебно-практической задачи при изложении возможных способов организации эффективных параллельных вычислений (см. [[12], [60]]).
Для простоты изложения материала в качестве области задания D функции u(x,y) далее будет использоваться единичный квадрат
D={(x,y)D:0 x,y 1}
11.1. Последовательные методы решения задачи Дирихле
Одним из наиболее распространенных подходов к численному решению дифференциальных уравнений является метод конечных разностей (метод сеток) (см., например, [[6], [13], [60]]). Следуя этому подходу, область решения D можно представить в виде дискретного (как правило, равномерного) набора (сетки) точек (узлов). Так, например, прямоугольная сетка в области D может быть задана в виде (рис. 11.1)

где величина N задает количество внутренних узлов по каждой из координат области D.
Обозначим оцениваемую при подобном дискретном представлении аппроксимацию функции u(x,y) в точках (xi, yj) через uij. Тогда, используя пятиточечный шаблон (см. рис. 11.1) для вычисления значений производных, мы можем представить уравнение Пуассона в конечно-разностной форме

Данное уравнение может быть разрешено относительно u_ij:

uij=0,25(ui-1,j+ui+1,j+ui,j-1-h2fij).

по которому очередное k-е приближение значения u_ij вычисляется по последнему k-му приближению значений u_i-1,j и u_i,j-1 и предпоследнему (k-1)-му приближению значений u_i+1,j и u_i,j+1. Выполнение итераций обычно продолжается до тех пор, пока получаемые в результате итераций изменения значений uij не станут меньше некоторой заданной величины (требуемой точности вычислений). Сходимость описанной процедуры (получение решения с любой желаемой точностью) является предметом всестороннего математического анализа (см., например, [[6], [13], [60]]), здесь же отметим, что последовательность решений, получаемых методом сеток, равномерно сходится к решению задачи Дирихле, а погрешность решения имеет порядок h2.

Рис. 11.1.  Прямоугольная сетка в области D (темные точки представляют внутренние узлы сетки, нумерация узлов в строках слева направо, а в столбцах — сверху вниз)
Рассмотренный алгоритм (метод Гаусса – Зейделя) на псевдокоде, приближенном к алгоритмическому языку С++, может быть представлен в виде:
Алгоритм 11.1. Последовательный алгоритм Гаусса – Зейделя
// Алгоритм 11.1
do {
dmax = 0; // максимальное изменение значений u
for ( i=1; i<N+1; i++ )
for ( j=1; j<N+1; j++ ) {
temp = u[i][j];
u[i][j] = 0.25*(u[i-1][j]+u[i+1][j]+
u[i][j-1]+u[i][j+1]–h*h*f[i][j]);
dm = fabs(temp-u[i][j]);
if ( dmax < dm ) dmax = dm;
}
} while ( dmax > eps );
(напомним, что значения u_ij при индексах i,j=0,N+1 являются граничными, задаются при постановке задачи и не изменяются в ходе вычислений).

Рис. 11.2.  Вид функции u(x,y) в примере для задачи Дирихле
Для примера на рис. 11.2 приведен вид функции u(x,y), полученной для задачи Дирихле при следующих граничных условиях:

Общее количество итераций метода Гаусса – Зейделя составило 210 при точности решения eps=0,1 и N=100 (в качестве начального приближения величин uij использовались значения, сгенерированные датчиком случайных чисел из диапазона [-100, 100]).