Построено семейство разностных схем для численного решения многомерного уравнения первого порядка с запаздывающим аргументом. Изучены порядок аппроксимации, устойчивость и порядок сходимости. Приведен алгоритм решения, допускающий параллельную реализацию. На тестовых примерах показана хорошая слабая и сильная масштабируемость параллельной версии. Экспериментальные оценки порядков сходимости соответствуют теоретическим оценкам.
Увеличение производительности суперкомпьютеров приводит к повышению интереса к разработке
параллельных вариантов численных методов. Ниже мы приведем обзор некоторых подходов. Первый класс
методов основан на декомпозиции области решения, классическим является альтернирующий метод
Шварца для эллиптических уравнений, последовательный по своей природе. Чтобы распараллелить метод
Шварца, применяются разнообразные подходы: аддитивный метод Шварца, параллельные многосеточные
методы с предобуславливанием, параллельный метод Шварца с весами и т. д. [
Параллельная реализация метода полу-Лагранжевого типа для уравнения переноса рассматривалась
в [
Настоящая работа обобщает результат [
Рассмотрим многомерное уравнение переноса с запаздыванием: (1) (2) (3) здесь x = (x1; : : : ; xp); x 2 [0; X ]; и t 2 [t0; ] пространственная и временная переменные; u(x; t) искомая функция; ut(x; ) = fu(x; t + ); 6 < 0g предыстория искомой функции к моменту t; > 0 – величина запаздывания, коэффициент a > 0:
Пусть G = Qp =1[0; X ] p-мерный прямоугольник с границей : Заданы начальные и граничные условия u(x; t) = '(x; t); x 2 G; t 2 [t0
; t0]; u(x; t) = g(x; t); x 2 ; t 2 [t0; ]: Пусть f; '; g таковы, что (1)–(3) имеет единственное решение, обладающее той степенью гладкости, которая требуется для построения численного метода [4].
Обозначим через Q = Q[ ; 0) множество функций u( ), кусочно-непрерывных на [ ; 0) с конечным числом точек разрыва первого рода, непрерывных справа в точках разрыва.
Введем в Q норму kukQ = sup 2[ ;0) ju( )j. Для построения численного метода и обоснования его сходимости дополнительно предположим, что функционал f (x; t; u; v( )) определен на G [t0; ] R Q и Липшицев по двум последним аргументам: Рассмотрим разбиение отрезка [0; X ] на части с шагом h = X =N ; где N число отрезков разбиения по соответствующей координате. На G введем сетку !h = f(i1h1; : : : ; iphp); i = 0; : : : ; N g: Мы используем обозначение x(i) = (xi1; : : : ; xip) = (i1h1; : : : ; iphp); i = 0; 1; : : : ; N ; здесь индекс i это p-мерный вектор. Обратим внимание, что x 2 R это -я координата вектора x, в то время как x(i) 2 Rp узел в сетке !h.
Введем ! = ftj = t0 + ; j = 0; : : : ; M g, где = ( t0)=M: Для простоты изложения предполагаем, что = = m целое.
Наконец, определим сетку !h = !h ! = f(x(i); tj); x(i) 2 !h; tj 2 ! g: В узлах сетки !h будем искать приближенное решение задачи (1)–(3).
Обозначим через uij приближение функции u в узле: uij u(x(i); tj); где (x(i); tj) 2 !h : Функционал f (x(i); tj; u(x(i); tj); utj (x(i); )) может зависеть от значений u между узами, а потому для его вычисления необходима интерполяция сеточной функции uij: Для каждого момента времени tj и величины запаздывания 2 [ ; 0) либо tj + 6 t0; в этом случае полагаем u(x(i); tj + ) = '(x(i); tj + ), либо tj + > t0; в этом случае используем интерполяцию. fulgij = uli j maxf0; j mg 6 l 6 j : Определение 1 Отображение I; определенное на множестве Aij всех возможных дискретных предысторий и действующее по правилу
I : Aij ! Q[ minf ; tjg; 0] : fulgij 7! vi;j( ) = vi;j(tj + ); называется оператором интерполяции дискретной предыстории.
Для построения численного метода, обладающего нужными свойствами, достаточно применять кусочнолинейную интерполяцию. Именно, для дискретной предыстории fulgij определим
1 vi;j(tj + ) = (tl+1 tj ) uli + (tj + tl) uli+1 ; tl 6 tj + 6 tl+1: (4) Оператор кусочно-линейной интерполяции имеет 2-й порядок на точном решении, т. е. существуют константы C1 и C2 такие, что для всех i = (i1; : : : ; ip); i = 0; 1; : : : ; N ; j = 1; : : : ; M; и t 2 [maxf0; tj g; tj] выполняется: vi;j(t) u(x(i); t) 6 C1 maxf0;mjamxg6l6j uli u(x(i); tl) + C2 2: Рассмотрим разностные операторы, которые аппроксимируют первые производные по пространственным переменным, 4uij[ 1 ] 2h a fji[ 1 ] @g(x(i)[ 1 ];tj) @g(x(i)[ 1 ];tj) @x @t
+ 4uij ! 8 > > > > > > uij = < > > i[ 2 ] >> u > j > : 4uij[ 1 ] + 3uij 2h p
P =1; 6=
2h ; для i > 2: ; для i = 1; Здесь x(i)[ 1 ] сосед узла x(i); сдвинутый влево на один шаг h по соответствующей координате, т. е. узел (i1h1; : : : ; i 1h 1; (i 1)h ; i +1h +1; : : : ; iphp); uij[ 1 ] сеточное приближение функции u(x; t); вычисленное в узле (x(i)[ 1 ]; tj): Рассмотрим семейство разностных схем (0 6 s 6 1) для j = 0; : : : ; M 1 и i = (i1; : : : ; ip); i = 1; : : : ; N : i uj+1 uij + a Xp s =1
i uj+1 + (1 s) uij = fji (5) с начальными и граничными условиями ui0 = '(x(i); t0); vi;0(t) = '(x(i); t); t < t0;
uij = g(x(i); tj); x(i) 2 : Здесь fji = f x(i); tj; uij; vi;j( ) значение функционала f; вычисленное на приближенном решении, vi;j( ) результат интерполяции. Для построения численного метода мы дополнительно предполагаем, что g(x; t) дифференцируема. 3
Исследование сходимости
для всех i и j:
Определение 2 Обозначим "ij = u(xi; tj) uij; i = (i1; :::; ip); i = 1; :::; N ; j = 0; : : : ; M: Будем говорить,
что метод (5) сходится, если "ij ! 0 при h ! 0 и ! 0 для всех i и j: Будем говорить, что метод
p p
сходится с пордяком P hq1 + q2; если существует константа C такая, что k"ijk 6 C P hq1 + q2
=1 =1
=1
u(x(i); tj+1) + (1 s)
u(x(i); tj)) fji;
(6)
где fji = f (x(i); tj; u(x(i); tj); utj (x(i); )) – значение функционала f , вычисленное на точном решении.
Теорема 1 Пусть точное решение u(x; t) задачи (1) (3) трижды непрерывно дифференцируемо по x ;
дважды непрерывно дифференцируемо по времени t и все первые производные решения по x непрерывно
p
дифференцируемы по t, тогда невязка метода (5) имеет порядок P h2 + :
=1
Доказательство теоремы сводится к разложению функции u в ряд Тейлора в окрестности точки (x(i); tj);
для одномерного случая см. [
В силу нелинейного характера зависимости функционала f от решения и его предыстории, мы
применим аппарат общих схем с последействием, разработанный в [
Без ограничения общности мы рассмотрим случай однородных граничных условий u(x; t) = 0, x 2 ; t 2 [t0; ]. Действительно, замена вида u~(x; t) = u(x; t) g^(x; t); где g^(x; t) произвольная функция, обладающая нужной степенью гладкости, совпадающая с g(x; t) при x 2 ; t 2 [t0; ]; сводит изначально поставленную задачу к задаче с однородными граничными условиями.
Доказательство сходимости будем проводить методом вложения семейства (5) в общую разностную
схему с последействем [
Пусть p фиксировано. Для каждого tj 2 ! определим значения дискретной предыстории вектором ~yj; для этого упорядочим компоненты вектора uij в лексикографическом порядке:
~yj = u(j1;:::;1); u(j2;:::;1); : : : ; u(jN1;:::;1); : : : ; u(jN1;:::;Np) > 2 Yp; здесь j = 0; : : : ; M; знак > означает транспонирование, Yp линейное пространство, dim Yp = N1
N2 : : : Np: Например, при p = 1; имеем ~yj = uj1; uj2; : : : ; ujN1 > : последовательность матриц D ; определяется через предыдущую.
Рассмотрим N1 N1 матрицу D1; которая соответствует разностному оператору 1 : p Построим разностный оператор A : Yp ! Yp; который соответствует P . Для этого рассмотрим =1 = 1; : : : ; p; где каждая следующая (начиная со второй) рекуррентно D1 = 0 4 0 B 4 3
BB 1 4 a BB 0 1
B . 2h1 BB ..
BB 0 B@ 0 0 Определим матрицу E1 = 2ha2 I1; где I1 единичная N1 N1-матрица. Далее, определим матрицу D2; которая имеет блочно-трехдиагональную форму, каждый блок имеет размер N1 N1; и в каждой строке N2 блоков: B B
B D2 = BBB
B B B @ C C C C C : C C C C
A Доказательство. Матрица D1 имеет положительные собственые числа 1(D1) = 2a=h; 2(D1) = : : : = = N1(D1) = 3a=2h; следовательно, положительно определена. Для всякого p матрица Dp нижнетреугольная с положительными диагональными элементами. Все ее главные миноры равны произведению соответствующих диагональных элементов, следовательно, положительны. Матрица Dp положительно определена по критерию Сильвестра. Лемма доказана.
Перепишем систему (5) в форме
~yj+1 ~yj + sA ~yj+1 + (1 s)A ~yj = F~j:
Определим линейный оператор B = I + s A (I единичный оператор соответствующей размерности)
и запишем (7) в виде двуслойной разностной схемы в каноническом виде [
B ~yj+1 ~yj + A ~yj = F~j: ~yj+1 = S ~yj + B 1 F~j; В силу леммы 1, оператор A положительно определен, следовательно, B положительно определен. Таким образом, B обратим, и можно записать (8) в виде где S = (I B 1 A) так называемый оператор продвижения на шаг.
В пространстве Yp определим скалярное произведение и энергетическую норму
N1
Np X i=1 (y; u) = h1 h2 hp yiui; kykYp = p(Ay; y); далее определим соответствующую индуцированую операторную норму. Определение 4 Разностную схему (8) будем называть устойчивой, если kSkYp < 1: Отметим, что в [13, стр. 324–330] дана эквивалентная формализация устойчивости двуслойных разностных схем. Теорема 2 Если s > 1=2, разностная схема (8) устойчива. 0 0 0 . . .
0
Dp 1 + 3Ep 1
1
C
C
C
C ;
C
C
C
A
(7)
(8)
Доказательство. Рассмотрим (8) с позиции операторно-разностных уравнений и применим методы
проверки устойчивости двуслойной разностной схемы [
B^ ~yj+1 ~yj + A^ ~yj = F~^j: Метод (9) усточив в энергетической норме тогда и только тогда, когда 2B^ > A^ [13, стр. 333, теорема 1]. Это требование эквивалентно тому, что A 1+ E(s 0:5) > 0: Поскольку A 1 положительно определенный, последнее неравенство выполняется для любых , как только s > 0:5. Теорема доказана.
Важной особенностью предлагаемого метода является то, что условие устойчивости s > 1=2 не накладывает ограничений на шаги разбиения, как это имеет место быть в условиях Куранта–Леви.
Методом вложения в общую схему [
Параллельная реализация и численные эксперименты Для простоты изложения все конструкции и обозначения в этом разделе вводятся для случая p = 2; по этой причине результаты численных экспериментов также приводятся для случая p = 2: Для нахождения uij1;i2 при i1 > 2; i2 > 2; требуется найти решение в 9 узлах сетки (рис. 1а). Такой тип рекуррентных информационных зависимостей позволяет последовательно находить решения. Соответственно, метод (5) может быть запрограммирован в виде трех вложенных циклов: внешний по времени, два внутренних по пространству. (9) Рис. 1: а) сетка и информационные зависимости; б) расположение блоков в пространстве и информационная зависимость между блоками.
Реорганизуем вычисления таким образом, чтобы решения в некоторых узлах сетки можно было находить параллельно. Разрежем сетку на K1 K2 пространственных колодцев. Назовем (k1; k2)-м, k1 = 1; :::; K1; k2 = 1; :::; K2; пространственным колодцем множество узлов сетки (i1; i2; j) i1 = (k1 1)N1=K1; :::; k1N1=K1 1; i2 = (k2 1)N2=K2; :::; k2N2=K2 1; j = 1; :::; M : Задачу нахождения сеточного решения разделим между процессами; именно, пусть (k1; k2)-й процесс вычисляет решение в (k1; k2)-м пространственном колодце, k1 = 1; :::; K1; k2 = 1; :::; K2; т. е. требуется организовать K1 K2 потоков.
Поток (1; 1) не имеет информационных зависимостей от других потоков. Для нахождения решения на (j +1)-м временном слое, процесс (k1; k2); k1 > 2; k2 > 2; должен получить от процесса (k1 1; k2) значения uij1;i2; где i1 = (k1 1)N1=K1 2; (k1 1)N1=K1 1; i2 = (k2 1)N2=K2; :::; k2N2=K2 1; и от процесса (k1; k2 1) значения uij1;i2; где i1 = (k1 1)N1=K1; :::; k1N1=K1 1; i2 = (k2 1)N2=K2 2; (k2 1)N2=K2 1; в этом состоит информационная зависимость между потоками.
Для каждого фиксированного набора k1; k2; j назовем блоком Bk1;k2 cледующее множество узлов: j
= 2u(x; y; t) u(x; y; t 1) x y t + 2; где 0 6 x 6 2; 0 6 y 6 2; 0 6 t 6 4; с начальными и граничными условиями u(x; y; t) = x + y + t; 0 6 x 6 2; 0 6 y 6 2; 1 6 t 6 0; u(0; y; t) = y + t; u(x; 0; t) = x + t; 0 6 x 6 2; 0 6 y 6 2; 0 6 t 6 4: Легко проверить, что поставленная начально-краевая задача имеет точное решение u(x; y; t) = x + y + t: В силу того, что решение линейно зависит от всех аргументов, а метод (5) имеет порядок второй по пространству и первый по времени, погрешность численного решения во всех экспериментах будет мала, порядка 10 14:
Таблица 1 содержит результаты численных экспериментов. Число отрезков разбиения по пространственным переменным выбрано одинаковым и изменяется от 64 до 1024, число отрезков разбиения по времени постоянно и равно 100. Число процессов указано в первой строке таблицы и меняется от 1 до 64. Ячейки таблицы содержат медианное время счета в минутах. Ускорение рассчитывалось по отношению к последовательному, т. е. однопроцессорному, варианту.
Для оценки качества распараллеливания сделана оценка сильной масштабируемости. Сильная
масштабируемость характеризует изменение времени решения задачи с увеличением количества процессоров (или
вычислительных узлов) при неизменном общем объеме задачи [
Чтобы проиллюстрировать качество численного метода, рассмотрим уравнение где 2 < x < 2; 5 < y < 5; 0 < t 6 4 ; с начальными и граничными условиями
1 u(x; y; t) = 1 + x2 + y2 + sin t; 2 6 x 6 2; 5 6 y 6 5; 6 t 6 0; u( 2; y; t) = 5 +1 y2 + sin t; u(x; 5; t) = 26 +1 x2 + sin t; 2 6 x 6 2; 5 6 y 6 5; 0 6 t 6 4 : Таблица 2: Погрешность численного решения в зависимости от шагов Вариант 1 2 3 4 5 6 7 h 1=5 1=10 1=10 1=20 1=5 1=10 1=20
=20 =20 =40 =40 =200 =200 =200 diff 0.0902 0.0919 0.0481 0.0495 0.0530 0.0135 0.0038 В экспериментах № 5 7 погрешность, обусловленная дискретизацией по времени, мала в сравнении с погрешностью, обусловленной дискретизацией по координате; анализ поведения погрешности показывает квадратичную сходимость по пространственным переменным: уменьшение шага в два раза приводит к уменьшению погрешности примерно в два раза. Анализ таблицы показывает, что только согласованное уменьшение шагов приводит к уменьшению погрешности. Так, в экспериментах № 7 8 уменьшение h в два раза не приводит к соответствующему уменьшению погрешности, так как общая погрешность уже в большей степени обусловлена дискретизацией по времени.
Согласно теореме 2, при s = 0:8 разностная схема (5) является устойчивой при любом соотношении шагов, однако в силу некорректности операции численного дифференцирования при уменьшении h аппроксимация производных по пространству в (5) становится более чувствительной к ошибкам машинного округления, что приводит к росту погрешности. Согласованное с h уменьшение является своего рода регуляризатором, который не позволяет вычислительным ошибкам расти и накапливаться. Эксперименты № 1 4 служат тому иллюстрацией.
Отметим, что порядок численного интегрирования для вычисления операторов распределенного запаздывания должен быть согласован с порядком метода (5). [2] S.A. Gourley, R. Liu, J. H. Wu. Spatiotemporal Distributions of Migratory Birds: Patchy Models with
Delay. SIAM J. Appl. Dyn. Syst., 9(2): 589–610, 2010. [3] T. Luzyanina, J. Cupovic, B. Ludewig and G. Bocharov. Mathematical models for CFSE labelled lymphocyte dynamics: asymmetry and time-lag in division. J. Math. Biol., 69(6-7): 1547–83, 2013. [4] J. Wu. Theory and applications of partial functional differential equations. New York: Springer–Verlag, 1996. [5] L. Tavernini. Finite difference approximation for a class of semilinear Volterra evolution problems. SIAM
J. Numer. Anal., 14(5): 931–949, 1977. [6] B. Zubik-Kowal. Stability in the Numerical Solution of Linear Parabolic Equations with a Delay Term. BIT
Numerical Mathematics, 41(1): 191–206, 2001. [7] Z. Kamont and M. Netka. Numerical method of lines for evolution functional differential equations. J.
Numer. Math., 19(1): 63–89, 2011. [8] Z. Kamont and K. Przadka. Difference methods for first order partial differential-functional equations with initial-boundary conditions. Computational Mathematics and Mathematical Physics, 31(10): 1476–1488, 1991. [9] Z. Kamont and K. Kropielnicka. Implicit difference methods for evolution functional differential equations.
Siberian Journal of Numerical Mathematics, 14(4): 361-379, 2011. [10] W. Czernous and Z. Kamont. Comparison of explicit and implicit difference methods for quasilinear functional differential equations. Appl. Math., 38(3): 315-340, 2011.
Parallel variant of numerical algorithm for solving a multidimensional advection equation with time delay
Svyatoslav I. Solodushkin1;2, Arsen A. Sagoyan1, Irina F. Iumanova1 1 – Ural Federal University (Yekaterinburg, Russia) 2 – Krasovskii Institute of Mathematics and Mechanics (Yekaterinburg, Russia)
Keywords: parallel numerical methods, difference scheme, advection equation, time delay, parallel programming.