Исследование двухсеточного метода повышенной точности для эллиптического уравнения реакции - диффузии с пограничными слоями Текст научной статьи по специальности «Математика»

____________УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА

Том 157, кн. 1 Физико-математические науки

2015

УДК 519.63

ИССЛЕДОВАНИЕ ДВУХСЕТОЧНОГО МЕТОДА ПОВЫШЕННОЙ ТОЧНОСТИ ДЛЯ ЭЛЛИПТИЧЕСКОГО УРАВНЕНИЯ РЕАКЦИИ-ДИФФУЗИИ С ПОГРАНИЧНЫМИ СЛОЯМИ

С.В. Тиховская

Аннотация

Исследуется двухсеточный метод для линейного эллиптического уравнения с малым параметром е при старшей производной. Рассматривается е-равномерно сходящаяся разностная схема на сетке Шишкина. Для разрешения разностной схемы используется двухсеточный метод с е-равномерной интерполяционной формулой. Для повышения точности схемы применяется экстраполяция Ричардсона. Приводятся результаты численных экспериментов. При реализации двухсеточного алгоритма предлагаются различные итерационные методы.

Ключевые слова: эллиптическое уравнение реакции-диффузии, сингулярное возмущение, сетка Шишкина, двухсеточный метод, экстраполяция Ричардсона, равномерная сходимость.

Введение

Рассматривается сингулярно возмущенная задача типа реакции-диффузии для линейного эллиптического уравнения в единичном квадрате. Такого рода задачи возникают в различных приложениях, и интерес к ним высок (см., например, [1-5] и приведенную там библиографию). Известно [6], что при выполнении условий согласования [2, 7] решение такого типа задач имеет регулярные экспоненциальные пограничные слои ширины порядка O(y/s) вблизи границы. Известно также [7, 8], что решение может иметь сингулярности в угловых точках области. При исследовании сингулярно возмущенных задач основным аспектом остается вопрос е-равномерной сходимости разностных схем. Один из способов достижения равномерной сходимости является использование сгущающихся сеток [1, 4, 6, 9]. В настоящей работе рассматривается классическая пятиточечная разностная схема на сетке Шишкина, имеющая в случае выполнения условий согласования е-равномерную точность порядка O(ln2 N/N2). Имеются также оценки и в случае более слабых ограничений на гладкость решения [7, 10].

Рассматриваемая разностная схема представляет собой систему линейных уравнений, которую можно решить на основе итераций. Двухсеточный метод является быстрым и эффективным методом для решения краевых задач, в том числе и нелинейных (см., например, [5, 11-14] и приведенную там библиографию). Идея двухсеточного метода заключается в предварительном решении задачи на более грубой вспомогательной сетке с последующей интерполяцией найденного решения на исходную сетку. Найденное на основе интерполяции сеточное решение далее принимается за начальное приближение для итераций на исходной сетке, что и приводит к уменьшению числа арифметических действий. Таким образом, возникает необходимость в е-равномерной интерполяции разностного решения с сохранением

60

ИССЛЕДОВАНИЕ ДВУХСЕТОЧНОГО МЕТОДА...

61

полученной точности схемы на вспомогательной сетке [14-16]. При использовании двухсеточного алгоритма также можно практически без дополнительных вычислительных затрат повысить точность решения разностной схемы, если использовать решение на вспомогательной сетке для построения численного решения с помощью метода экстраполяции Ричардсона [14, 15, 17].

1. Постановка задачи

Рассмотрим линейную сингулярно возмущенную эллиптическую задачу

euxx + euyy - c(x,y)u = f (x,y), (x,y) € i; u(x,y) = g(x,y), (x,y) € T,

где i = (0,1)2 , Г = di = i \ i, c, f, g - достаточно гладкие функции,

c(x,y) > 2y > 0, £ > 0. (2)

Известно [1, 2, 6], что при выполнении условий (2) решение задачи (1) является равномерно ограниченным и имеет пограничные слои у границ x = 0, x = 1 и y = 0, y = 1. _

В соответствии с [2], если f,c € C4,A(Q), g e C4’X(di) и выполнены условия согласования [8], то решение задачи (1) имеет регулярные экспоненциальные пограничные слои ширины 0(т/£) вблизи дQ и может быть представлено в виде

4 4

u = v + Wi + XI Zi,

i=1 i=1

где

|| v(kj)|| < C (l + £1-k/2-j/2^ , 0 < k + j < 4,

\w1(x,y)\ < Ce-^, \w2(x,y)\ < Св-^х,

\w3(xy)\ < Ce—fY/E(1-y), \w4(xy)\ < Ce^Y/e(1-x), max{||w.k’j)||, ||z(k’j)||} < C£-k/2-j/2, 0 < k + j < 4,

IIw^II < C (l + £1-k/2) , * = 1, 3,

К(0’к)|| < C (l + £1-k/2) , * = 2, 4, 0 < k + j < 4,

IIzi(x,y)II < Ce^Y/eye-/Y/ex, ||z2(k, j)I| < Ce-VY/~e(1-y)e-VY/ex,

IIz3(k,j)II < Ce-VY^(1-y)e-VY^(1-x), IIz4(k,j)II < Ce-/Y/eye-VY^(1-x),

f (k,f) = 9(k+j)f.

dxkdyj

Здесь и далее через C (возможно с индексами) обозначаем положительные постоянные, не зависящие от £ и числа шагов.

62

С.В. ТИХОВСКАЯ

Для достижения е -равномерной сходимости разностной схемы для задачи (1) используем схему центральных разностей на сетке Шишкина.

Зададим в области Q кусочно-равномерную сетку [6]:

Qn = {(xi,yj), i, j = 0, N, hi = xi-Xi-i,Tj = yj-yj-1,xo = 0,xN = 1,y0 = 0,yN = 1},

где

hi =

4ax N ’ 1 < i < N ~4 ’ 3N ~T < i < N; , 2(1 - 2Ox) h = N ’ N ~4 < i < 3N ~T ’

4oy N ’ 1 < j < N ~4 ’ 3N ~T <j < N; _ 2(1 - 2oy) j N ’ N ~4 <j < 3N T ’ (3)

j

12 re

ox = min < -, ln N

|1,^lnnJ, о <у< /Y

, mi, , . oy = min < -, —:— ln N 4 у, J’ y \4’ ix

На сетке Qn выпишем схему с использованием центральных разностей

^ / N N N N

2е I ui+1j - uid uid - ui_1,j

hi + hi+i \ hi

+

+

n / N N N N

2e t uij+i - uij ui j - Ui j_i

Tj + Tj+i V Tj+i

- c(xi, yj)ui j = f (xi,yj), (xi,yj) € Q

N,

u4 = g(xi,yj), (xi,yj) € Tn =^|Qn. (4)

Определим норму непрерывной функции z и норму сеточной функции zN соответственно по формулам

lzll = ™х_\z(x,y)\, \\z ||n = max \zij\.

(x,y)en O^ij^N-

h

j

Пусть [u]n^ - проекция функции непрерывного аргумента u(x, y) на сетку Qn . В соответствии с [2] выполнено неравенство

||uN - [u]nN\\n ^ ^i^N, An = ln2 n/n2. (5)

В настоящей работе исследуем реализацию схемы (4) двухсеточным методом для задачи (1) с повышением точности на основе экстраполяции Ричардсона и использованием е -равномерной интерполяционной формулы.

2. Двухсеточный метод

Решение схемы (4) можно находить на основе итераций

u(m+i) = G(u(m)), (6)

где u(0) задано. Матрица системы (4) является M-матрицей, что обеспечивает сходимость многих итерационных методов [18, 19]. Пусть для решения системы (4) применяется сходящийся итерационный метод, для которого справедлива оценка сходимости итераций

||u(m+i) - uNIn < q\u(m) - uNIn, q< 1. (7)

ИССЛЕДОВАНИЕ ДВУХСЕТОЧНОГО МЕТОДА...

63

Пусть ||u(0) — uN||N ^ 5, тогда, рассуждая аналогично [15], получим, что для разрешения схемы (4) потребуется

Mn « dN2 logg(An/5) (8)

арифметических действий, если предположить, что для реализации одной итерации необходимо выполнить dN2 арифметических действий (пропорционально N2 ).

Теперь рассмотрим двухсеточный метод для нахождения решения разностной схемы (4). Введем сетку Qn такую же по структуре и с такими же параметрами ах и ау, как Qn , только с намного меньшим количеством узлов n, n ^ N, получив таким образом вложенные сетки. Считаем, что N = kn.

В соответствии с идеей двухсеточного метода предварительно с использованием итерационного метода (6) решаем задачу (1) на сетке Qn.

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

Пусть Ij(vn,x,y) - интерполянт сеточной функции vn в области Q и справедлива оценка погрешности

II1J ([u]^^,ctw ,x,y) ull ^ C2Aint,n, Aint,n ^ An, (9)

где u(x,y) - решение задачи (1).

Вычислим, рассуждая аналогично [15], необходимое количество арифметических действий двухсеточного метода:

MnN ~ dn2 logg + dN2 logg AAN + Jn, (10)

5 An

где Jn - количество арифметических действий, необходимых для интерполяции.

Сравнивая (8), (10), оценим выигрыш в числе арифметических действий при использовании двухсеточного метода:

Mn — MnN « d(N — n2) logg (An/5) — Jn.

Теперь оценим относительный выигрыш в количестве арифметических действий. Считая, что для реализации метода интерполяции необходимо Jn = dj (N — - n) действий (пропорционально N - n), и учитывая (8), (10), получим

MnN

Mn

n2\ d logq (AJ5) — 1

N 2 d logq (AN/5)

(11)

Покажем, что ty(N) из представления (11) удовлетворяет неравенству

0 < V(N) < 1.

Представим ty(N) как

*(N )= (log, f Ajn

dj

d

logq[ AN

(12)

Тогда

3No : yN ^ No An < 5,

64

С.В. ТИХОВСКАЯ

следовательно, знаменатель в (12) положителен. А так как

Зпо : Уп ^ n0 Дп < Sqdl/d+1,

что означает выполнение неравенства

logg (Дп/д) > di/d + 1,

следовательно, числитель в (12) также положителен. Оценка ^(N) < 1 следует из того, что ДN < Дп. Таким образом, получаем, что

. MnN < п2

MN ^ N2'

Найдем теперь предел ty(N) при N . Учитывая, что N = C2/N2 , где £ = (ра)/(2^ё), получим

kn и ДN =

lim ty(N)

N —

lim

N —— ж

± ± Дп ln q Дп д

_д_ Дл

ln q ДN д

1.

Следовательно, относительный выигрыш в количестве арифметических действий от использования двухсеточного метода ограничен величиной 1 - п2 /N2 , но при увеличении числа узлов N стремится к этой же величине.

3. Экстраполяция Ричардсона

Исследуем экстраполяцию Ричардсона для повышения точности разностной схемы при использовании двухсеточного метода. Метод экстраполяции Ричардсона для сингулярно возмущенных эллиптических задач исследовался, например, в работах [4, 17, 20]. Рассмотрим случай, когда вспомогательная сетка Шишкина имеет вид (3) и содержит в k раз меньше сеточных интервалов по каждому направлению, чем исходная сетка, где k > 1 - целое число.

Решение uN схемы центральных разностей (4) вычисляется на сетке Qn . Для повышения точности будем также использовать решение этой схемы на сетке Qn, которая содержит п сеточных интервалов и имеет те же параметры ax, ау, что и сетка Qn . Эти сетки вложены так, что Qn = {(Xi,Ym)} С Qn = {(xi ,Vj)}. Очевидно, что сетку Qn можно получить из Qn делением каждой ячейки на k равных частей по каждому направлению.

Обозначим решение схемы (4) на Qn через un, и пусть

п2 1 N2 k2

kn — -»,л ~ — m т 1 kN

N2 п2

k2 1

N2 п2 k2 1

В соответствии с методом экстраполяции Ричардсона зададим функцию unN на сетке Qn , приближающую решение u(x, у) с более высоким порядком точности, чем uN . Для этого сначала в узлах вспомогательной сетки Qn определим сеточную функцию unN следующим образом:

unN(Xl,Ym) = knun(Xt,Ym) + kNuN(Xl,Ym), (Xi,Ym) e Qn'

В узлах исходной сетки Qn , не совпадающих с узлами сетки Qn, зададим сеточную функцию unN (xi ,yj), используя интерполяцию. Тогда для каждого узла (xi,yj) e Qn определим

unN (xi,yj )= In([unN ]n„ ,xi,yj ), (13)

где I^(vn,x,y) - интерполянт сеточной функции vn в области Q.

Отметим, что точность интерполяционной формулы в (13) должна быть е-рав-номерной, порядка, не ниже порядка точности сеточной функции unN в узлах сетки (Xi,Ym) e Qn, полученного в результате применения метода Ричардсона.

ИССЛЕДОВАНИЕ ДВУХСЕТОЧНОГО МЕТОДА...

65

4. Результаты численных экспериментов

Рассмотрим краевую задачу

£ПХХ + £Uyy - 2u = f (x,y), (x,y) € Q, u(x,y) = 0, (x,y) € Г,

где f соответствует решению

(14)

u(x,y)

1 — e

Ej (l - e-(i-x)/^j (l - e-y^ij (l - e-(1-y)^ij

1-

1

1-

1

+

+ sin(^x) sin(^y).

(15)

В соответствии с (14) зададим аХ и ay в (3), используя значение у, = 1. Решение задачи (14) находим на основе схемы (4). Начальное приближение для используемых итерационных методов задаем на основе известных граничных условий дифференциальной задачи (14) следующим образом:

u(0)(xi,yj) = 0, (xi,yj) € Qn.

Итерационный метод на сетке Qn завершаем, если выполнено условие

\\LN ) - fN||n < 7An,

тогда будет выполнена оценка ||u(mN) - uN||n ^ An .

Остановимся подробнее на выборе интерполяционных формул в (9), (13). Для этого обобщим интерполяционную формулу из [16] на двумерный случай и исследуем ее при к = 3,5 и различном выборе функций Ф^) и 0(y). Пусть Lk(v,Zp,... ,zp+k-i,z) - многочлен Лагранжа с узлами zp,..., zp+k-i и [zp,..., zp+k-i]v - разделенная разность для функции v(z), тогда

!ф,е,к (u, x, y) = Lk-i(u^,k ,yi,..., yi+k-2,y)+

[yi,...,yi+k-i]^,k f^( ) r (t^ )\

--------------ГУТ еЫ - Lk-i(®,yi,... ,yi+k-2,y) ,

[yi ,...,yi+k-i]® V )

+

(16)

где

UФ!k (u,x) Lk-i(u: xm> . . . , xm+k-2 ,x^ +

. [xm, . . . , xm+k-i]u (ф( ) у

+ 7--------------TTfT Ф(x) - Lk

[^mi . .., xm+k-i]Ф '•

i (ф, xm, . . . , xm+k-2, x^j ,

(x,y) € [xm,xm+k-i] X [xi ,xi + k-i].

Отметим, что в случае, когда Ф^) = xk—i и 0(y) = yk—i, получаем двумерную формулу полиномиальной интерполяции.

В работе [21] показано, что при к = 2 одномерная формула (16) имеет е-равномерную точность 0(ln N/N), а при к = 3 имеет точность 0(ln2 N/N2). Проводя рассуждения, аналогичные [14], можно обобщить полученные оценки на двумерный случай интерполяционной формулы. В соответствии с (5) и (9) случай к = 2 не рассматриваем, так как погрешность метода интерполяции больше погрешности схемы (4).

В (13) необходима интерполяционная формула повышенного порядка точности, поэтому остановимся также на случае к = 5. В работе [21] показано, что при этом

66

С.В. ТИХОВСКАЯ

формула будет иметь точность O(ln4 N/N4). Результаты численных экспериментов показывают, что данная оценка является е-равномерной.

Выпишем формулу (16) при к = 5 для каждого узла (xi, yj) G &n,o-n , принадлежащего некоторой ячейке Si,m = [Xi_2,Xi+2] х [Ym_2,Ym+2], где (Xp,Yq) G &n,aN . Считаем, что n кратно четырем и не меньше 16. Имеем

т ([ nN ] ) U1N (xi,Ym_1) - u1N (xi,Ym_2)

Ir (u К ,xi,yj ) = --------^-------Г7----)-------

(Ym_1 Ym_2 )

ulN (Xi, Ym) - 2u^N (Xi, Ym-i) + U^N (xh Ym_2)

(yj Ym_2) +

+

(yj Ym_2)(yj Ym_ 1 ) +

(Y^ Y^_2) (Ym Ym_1)

+ (u^N(xi,Ym+i) - 3ulN(xi,Ym)+3ulN(xi,Ym_i)-

(yj Ym_2)(yj Ym_1)(yj Ym)

- unN (xi,Ym_2)

j - m_2 j - m_1 j - m (Ym+1 Ym_2) (Ym+1 Ym_1) (Ym+1 Ym)

+

i Qm_1 Qm_2 / \ Qm 2Qm_1 + Qm_2 ,

+ Gs[ - Y---XV----(yj - Ym_2 )--—----------^2--y -

Ym 1 - Ym 2

2 (Ym - Ym_ 1)

■лг w -лг \ Qm+1 3^m + 3^m_1 Q m_2 / -*r w

- Ym_2)(yj - Ym_ 1) XX? V \3 (yj - Ym_2)(yj-

6(Ym - Ym_1)

- Ym_1)(yj - Ym) + Q(yj ) - Qm_2^j + UnN ^, Ym_2 ), (17)

где

unnN (xi,yj )= unN (Xi_2,yj) +

,nN

(Xi_1,yj) - unN(Xi_2,yj)

Xi 1 - Xi 2

(xi - Xi_2 ) +

+ (unN (Xi, yj) - 2unN (Xi_1,yj)+ unN (Xi_2,yj )

(xi - Xi_2) (xi - Xi_1)

(Xi - Xi_2)(Xi - Xi_1) + (<unN (Xi+1,yj) - 3unN (Xi, yj ) + 3unN (Xi_1,yj)-

(xi Xi_2')(xi Xi_1)(xi Xi)

+

- unN(Xi_2,yj))

(Xi+1 - Xi_2) (Xi+1 - Xi_1) (Xi+1 - Xi)

+

Xl Ф_1 - Ф_2 , s ^i - 2^i_1 + Ф_2

+ GX -

Ч Xi_1 - Xi_2

(„ v ) ф - 2Ф_1 + Ф_2 („ v S(„ v )

(xi Xi_2) ъ /-*r -гг \2 (xi Xi_2)(xi Xi_1)

2(Xi - Xi_1)2

Ф1 + 1 - ЗФ1 + ЗФ1_1 - Ф1_2 6(Xi - Xi_1 f

(xi - Xi_2)(xi - Xi_1)(xi - Xi) + Ф(xi) - Ф_2^),

G5

GX

urjN(xj,Ym+2) - 4urjN(Xj,Ym+1) + 6ulN(xi, Ym) - 4urjN(Xi,Ym-1)+ urjN(Xj,Ym-2) ©m + 2 4©m+l + 6©m 4©m-l + ©m —2

= unN(Xi+2,yj) - 4unN(Xi+!,yj) + 6unN(Xi,yj) - 4unN(Xi-i,yj) + unN(Xi-2,yj) Ф;+2 — 4Ф1+1 + 6Ф1 — 4Ф—1 + Ф1-2 ’

Ф1 = Ф^1), Qm = e(Ym).

Так как сетка (3) является кусочно-равномерной, то условия, налагаемые на количество узлов вспомогательной сетки в (17), гарантируют равномерный шаг внутри каждой ячейки. А сгущение сетки (3) в областях больших градиентов интерполируемой функции обеспечивает оценки точности равномерные по е.

ИССЛЕДОВАНИЕ ДВУХСЕТОЧНОГО МЕТОДА...

67

Табл. 1

Нормы погрешности 7ф,е,з при Ф(х) = х2 и ©(у) = у2

£ N

32 64 128 256 512 1024

1 4.98e-4 6.08e-5 7.48e-6 9.27e-7 1.15e-7 1.44e-8

10-1 3.00e-4 3.78e-5 4.65e-6 5.77e-7 7.18e-8 8.96e-9

10-2 9.18e-3 1.38e-3 1.98e-4 2.65e-5 3.44e-6 4.37e-7

10-3 1.05e-1 3.28e-2 4.95e-3 7.61e-4 1.06e-4 1.39e-5

10-4 1.05e-1 3.72e-2 8.65e-3 1.93e-3 3.91e-4 7.23e-5

10-5 1.05e-1 3.72e-2 8.65e-3 1.93e-3 3.91e-4 7.23e-5

10-6 1.05e-1 3.72e-2 8.65e-3 1.93e-3 3.91e-4 7.23e-5

Табл. 2

Нормы погрешности 1ф,о,ь при Ф(х) = х4 и ©(у) = у4

£ N

32 64 128 256 512 1024

1 9.11e-6 2.74e-7 8.23e-9 2.51e-10 7.73e-12 2.40e-13

10-1 1.39e-5 4.08e-7 1.27e-8 4.01e-10 1.27e-11 3.98e-13

10-2 1.21e-3 5.46e-5 1.93e-6 6.88e-8 2.30e-9 7.43e-11

10-3 3.43e-2 6.24e-3 4.52e-4 1.69e-5 6.35e-7 2.21e-8

10-4 3.43e-2 7.32e-3 1.06e-3 8.96e-5 5.50e-6 3.38e-7

10-5 3.43e-2 7.32e-3 1.06e-3 8.96e-5 5.52e-6 3.38e-7

10-6 3.43e-2 7.32e-3 1.06e-3 8.96e-5 5.59e-6 3.45e-7

В табл. 1 при различных значениях N и £ приведены нормы погрешности метода интерполяции (16) при к = 3 в случае Ф(х) = х2, 0(у) = у2 для функции (15). Отметим, что в случае

Ф(х) = e-x/^E + e-(l-x)/^ и 0(у) = e-y/^E + e~(1~y)/^i

погрешность получается того же порядка.

Из табл. 1 следует, что точность формулы (16) при к = 3 в случае Ф(х) = х2 и 0(у) = у2 даже выше, чем теоретически обоснованная оценка 0(ln2 N/N2), и близка к 0(ln3 N/N3). Следовательно, если применить данную интерполяционную формулу, то будет выполнено неравенство (9).

В табл. 2 при различных значениях N и £ приведены нормы погрешности метода интерполяции (16) при к = 5 в случае Ф(х) = х4, 0(у) = у4 для функции (15).

В табл. 3 при различных значениях N и £ приведены нормы погрешности метода интерполяции (16) при к = 5 в случае

Ф(х) = e-x^E + e~(1~x)^/l и 0(у) = e-y^E + e-(1-y)^E для функции (15).

Из табл. 2 и 3 следует, что точность формулы (16) при к = 5 в случае выбора Ф(х) и 0(у) с учетом вида областей больших градиентов интерполируемой функции выше, чем 0(ln4 N/N4), и близка к 0(ln5 N/N5). Поэтому в (17) возьмем

Ф(х) = e-x/^ + e-(1-x)/^E и 0(у) = e-y/^E + e-(1-y)/^E.

Исследуем реализацию схемы (4) на основе явного метода Зейделя [19]. Пятиточечную схему (4) можно в общем виде представить как

N I г N

ai,j Ui-1,j + bi,j Ui,j-

+ c- uN + d uN — e- uN

1 ' ci,j ui+1,j ' di,j ui,j+1 ei,j ui,j

f

N

i,j5

0 <i, j < N.

68

С.В. ТИХОВСКАЯ

Табл. 3

Нормы погрешности /Ф,е,Б при Ф(ж) = e-x/^ + e-(1-x)A/£ и 0(y) = = eTv^1 + е~(1~у)/л/ш

е N

32 64 128 256 512 1024

i 9.52e-6 2.93e-7 8.86e-9 2.70e-10 8.38e-12 3.90e-13

10-1 1.27e-5 4.37e-7 1.38e-8 4.22e-10 1.31e-11 4.08e-13

10-2 4.14e-5 1.29e-6 4.16e-8 1.31e-9 4.10e-11 1.28e-12

10-3 2.82e-4 5.68e-6 1.61e-7 4.90e-9 1.50e-10 4.62e-12

10-4 2.78e-3 2.41e-4 7.76e-6 1.79e-7 4.17e-9 1.01e-10

10-5 4.03e-3 4.86e-4 3.55e-5 1.64e-6 4.80e-8 1.23e-9

10-6 4.47e-3 5.51e-4 4.35e-5 2.84e-6 1.65e-7 6.56e-9

Табл. 4

Количество итераций односеточного и двухсеточного методов Зейделя, е = 10-4

n N

32 64 128 256 512

16 13(7) 36(7) 113(7) 385(7) 1357(7)

32 30(15) 98(14) 341(13) 1219(13)

64 82(39) 296(35) 1082(33)

128 247(121) 942(108)

256 788(407)

17 45 139 463 1605

Тогда векторно-матричная запись метода Зейделя будет следующей:

u(m = d-1 f + Lu(m) + Uu(m-1^ ,

где

(Lv)i,j ai,j Vi — l,j + bi,j Vi,j-1, (Dv)i,j ei,j Vi,j, (Uv)i,j Ci,j Vi+1,j + di,jvi,j+1'

Результаты численных экспериментов для метода Зейделя при е = 10-4 приведены в табл. 4 и 5.

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

В табл. 5 при различных значениях N и n приведены нормы погрешности двухсеточного метода с экстраполяцией Ричардсона (17). В нижней строке для сравнения приведены нормы погрешности односеточного метода в зависимости от N. Из таблицы следует, что экстраполяция Ричардсона (17) повышает точность схемы до порядка не ниже 0(ln3 N/N3).

Помимо метода Зейделя в качестве итерационного метода использовался метод последовательной верхней релаксации [19]:

u(m = uD-1 (f + Lu(m) + Uu(m-1)) + (1 - u)u(m-1)

ИССЛЕДОВАНИЕ ДВУХСЕТОЧНОГО МЕТОДА...

69

Табл. 5

Нормы погрешности односеточного метода Зейделя и двухсеточного метода Зейделя с экстраполяцией Ричардсона, е = 10-4

n N

32 64 128 256 512

16 1.17e-3 4.23e-4 1.22e-4 2.98e-5 6.20e-6

32 2.74e-4 1.21e-4 4.62e-5 1.59e-5

64 4.16e-5 1.75e-5 6.71e-6

128 4.85e-6 1.87e-6

256 5.27e-7

8.15e-3 3.15e-3 1.10e-3 3.62e-4 1.15e-4

где ш - итерационный параметр, а также метод Писмана-Речфорда (продольнопоперечных прогонок) [19]:

U(m-1/2) = u(m-l) + т (_A1U(m-1/2) + A2U(m-1) — f) ,

u(m) = u(m-1/2) + т (Aiu(m-1/‘2) + A2U(m) - f) ,

где т - итерационный параметр,

(A1 v)i,j = ai,jvi-1,j ei,jvi,j + ci,jvi+1,j ,

(A2 v)i,j bi,j vi,j-1 ei,j vi,j + di,j vi,j+1'

Здесь ei,j = eij + e'ij, где ei j соответствует аппроксимации производных по направлению x, а e'i j - по направлению у.

Использованы также итерации в подпространствах Крылова [22-24], а именно стабилизированный метод бисопряженных градиентов (BiCGSTAB) и стабилизированный метод бисопряженных невязок (BiCRSTAB):

r0 = f — Au°, p0 = r0, n = 0,1,...

s

n

П — an)AP

n

v(«)

n

(rn, (AT)qr0) (Apn, (AT)qr0),

u

n+1

un + anq)p(nq) + Unsn,

шп

(Asn, sn) (Asn,Asn),

p(q)

pn+1

r(q)

'n+1

+ Aq)

e(q) = [an (rn+\ (AT)qr0)]

Pn [шп (rn, (AT)qr0)]

Здесь q = 0 и 1 для методов BiCGSTAB и BiCRSTAB соответственно, номера итераций указаны в данном случае нижними индексами [22].

Остановимся подробнее на выборе итерационного параметра ш в методе последовательной верхней релаксации. Согласно [19], оптимальное значение релаксационного параметра есть

Ш0

2

1 + V1 — Ps

(18)

где ps - коэффициент подавления ошибки метода Зейделя.

В этом случае можно апостериорно выбрать значение Ш0, используя следующий алгоритм автоматического определения итерационного параметра [19]. Сначала

70

С.В. ТИХОВСКАЯ

Табл. 6

Сравнение количества итераций, е = 10-4

Метод N

32 64 128 256 512

Зейделя 13(7) 17 30(15) 45 82(39) 139 247(121) 463 788(407) 1605

верхней релаксации 12(8) 15 23(15) 35 58(31) 95 106(81) 186 215(159) 330

Писмана-Речфорда 7(5) 11 14(12) 24 28(27) 54 58(60) 119 118(130) 259

BiCRSTAB 6(3) 7 10(6) 15 20(13) 31 41(25) 57 65(60) 98

BiCGSTAB 6(3) 7 11(6) 13 20(13) 28 36(26) 58 58(54) 93

в методе последовательной верхней релаксации проводятся итерации при значении ш = 1 и на каждой итерации вычисляется значение

Р

П

S

\\un _ un-1\\

\\un-l _ un-2\\ >

а также проверяется условие

\рП _ Pns-l \ < £1 (19)

при некоторой заданной малой величине £i. Когда неравенство (19) выполняется, величина рП принимается за приближенное значение ps, после чего по формуле (18) определяется значение шо, используемое на последующих итерациях. При этом предварительные итерации по методу Зейделя служат для получения достаточно приемлемого начального приближения.

В [25] показано, что предпочтительнее всего по количеству итераций выглядит автоматический выбор итерационного параметра ш.

В табл. 6 при £ = 10-4 проведено сравнение исследуемых итерационных методов. Двухсеточный метод использовался в случае n = N/2. Для метода последовательной верхней релаксации итерационный параметр ш выбирался в соответствии с описанным выше алгоритмом автоматического определения значения, для метода Писмана - Речфорда итерационный параметр т = 8/N.

Применение методов последовательной верхней релаксации, Писмана -Реч-форда, стабилизированных методов бисопряженных направлений и градиентов привело к уменьшению количества итераций по сравнению с методом Зейделя. Повышение точности на основе экстраполяции Ричардсона не зависит от конкретного итерационного метода.

В табл. 7 при различных значениях £ для стабилизированного метода бисопряженных градиентов приведено количество итераций и время выполнения двухсеточного (сверху) и односеточного (снизу) методов.

В табл. 8 при различных значениях £ и N для стабилизированного метода бисопряженных градиентов приведено преимущество по времени выполнения двухсеточного метода перед односеточным в процентном соотношении:

AT

100 1

Ttgm \

T

OGM

)

где Ttgm - время выполнения двухсеточным методом, Togm - время выполнения односеточным методом.

ИССЛЕДОВАНИЕ ДВУХСЕТОЧНОГО МЕТОДА...

71

Табл. 7

Количество итераций (Н) и время выполнения (Т) двухсеточного и односеточного методов

е = 10—1 £ = 10—2

N 256 512 1024 2048

Н 112(81) 207(185) 499(442) 1212(901)

185 442 901 2251

T 1.79 19.37 206.24 2244.99

2.56 36.55 317.42 3521.70

N 256 512 1024 2048

Н 121(121) 306(259) 584(570) 1443(1437)

259 570 1437 4623

T 2.01 28.18 244.97 2790.42

3.47 45.92 492.94 7194.04

£ = 10—3 е = 10—4

N 256 512 1024 2048

Н 36(26) 58(54) 108(97) 189(179)

58 93 207 434

T 0.69 5.94 46.79 371.31

0.88 7.94 74.27 684.25

N 256 512 1024 2048

Н 47(37) 72 87(72) 148 177(148) 430 318(430) 795

T 0.85 1.04 8.44 12.40 73.98 154.29 661.48 1246.48

Табл. 8

Преимущество ДТ двухсеточного алгоритма по времени выполнения в процентах

£ N

256 512 1024 2048

10—1 42.04 38.62 50.30 61.21

10—2 30.07 46.99 35.03 36.25

10—3 18.04 31.93 52.05 46.93

10—4 21.93 25.10 37.00 45.73

10—5 20.81 32.83 37.19 44.04

Табл. 9

Скорость сходимости в случае односеточного метода Писмана-Речфорда

£ N

32 64 128 256

10—1 2.001 2.000 2.000 2.000

10—2 1.988 1.997 2.000 2.001

10—3 1.508 1.965 1.991 1.998

10—4 1.370 1.517 1.604 1.656

10—5 1.370 1.517 1.604 1.656

CRt 1.356 1.474 1.555 1.615

Из табл. 7, 8 следует, что с увеличением N выигрыш от использования двухсеточного метода растет, что соответствует полученной оценке (11).

В табл. 9 приведена скорость сходимости схемы (4) в зависимости от е и N при ее реализации односеточным методом Писмана-Речфорда:

CR = log2 Dn, Dn = ||mn - [w]n*||-

D2N

В последней строке табл. 9 приведена теоретическая оценка скорости сходимости CRt схемы (4) в зависимости от N, соответствующая оценке погрешности (5).

В табл. 10 приведены значения скорости сходимости схемы (4) в зависимости от е и N при ее реализации двухсеточным методом Писмана - Речфорда с использованием экстраполяции Ричардсона при n = N/2. В последней строке табл. 10

72

С.В. ТИХОВСКАЯ

Табл. 10

Скорость сходимости в случае двухсеточного метода Писмана-Речфорда при использовании экстраполяции Ричардсона

e N

32 64 128 256

10-1 4.727 4.737 4.538 4.158

10-2 2.123 2.090 2.297 2.354

10-3 2.939 4.261 3.694 2.590

10-4 2.639 3.368 3.139 3.200

10-5 2.638 3.362 3.117 3.225

CRt 2.034 2.211 2.333 2.422

2.712 2.948 3.110 3.229

приведены теоретические оценки скорости сходимости CRt схемы (4) с экстраполяцией Ричардсона в зависимости от N, соответствующие оценке 0(ln3 N/N3) и 0(ln4 N/N4).

Итак, применение двухсеточного метода на сетке Шишкина приводит к выигрышу в количестве арифметических действий, а использование экстраполяции Ричардсона повышает точность разностной схемы на порядок. Исследовано обобщение одномерной £-равномерной интерполяции на двумерный случай. С увеличением числа узлов преимущество двухсеточного метода увеличивается. Исследованы методы Зейделя, последовательной верхней релаксации, Писмана-Речфорда и стабилизированные методы бисопряженных градиентов и направлений, наибольшую скорость сходимости среди которых показали два последних.

Автор благодарит профессора А.И. Задорина и профессора В.П. Ильина за полезные советы и интересные обсуждения при подготовке данной работы.

Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (проекты № 13-01-00618, 15-01-06584).

Summary

S.V. Tikhovskaya. Investigation of a Two-Grid Method of Improved Accuracy for the Elliptic Reaction-Diffusion Equation with Boundary Layers.

A two-grid method for the elliptic equation with a small parameter e multiplying the highest derivative is investigated. The e-uniformly convergent difference scheme on the Shishkin mesh is considered. To resolve the difference scheme, a two-grid method with e-uniform interpolation formula is used. To increase the accuracy of the scheme, the Richardson extrapolation in the two-grid method is applied. The results of numerical experiments are discussed. Various iterative methods for implementation of the two-grid algorithm are suggested.

Keywords: elliptic reaction-diffusion equation, singular perturbation, Shishkin mesh, two-grid method, Richardson extrapolation, uniform convergence.

Литература

1. Miller J.J.H., O’Riordan E., Shishkin G.I. Fitted Numerical Methods for Singular Perturbation Problems. - Singapure: World Scientific, 1996. - 163 p.

2. Clavero C., Gracia J.L., O’Riordan E. A parameter robust numerical method for a two dimensional reaction-diffusion problem // Math. Comput. - 2005. - V. 74, No 252. -P. 1743-1758. - doi: 10.1090/S0025-5718-05-01762-X.

ИССЛЕДОВАНИЕ ДВУХСЕТОЧНОГО МЕТОДА...

73

3. Roos H.G., Stynes M., Tobiska L. Robust Numerical Methods for Singularly Perturbed Differential Equations. - Berlin: Springer-Verlag, 2008. - 604 p.

4. Shishkin G.I., Shishkina L.P. Difference Methods for Singular Perturbation Problems. -Boca Raton: Chapman & Hall/CRC, 2009. - 408 p.

5. Angelova I.T., Vulkov L.G. A two-grid method on layer-adapted meshes for a semilinear 2-D reaction-diffusion problem // Lecture Notes in Computer Science. - 2010. - V. 5910. -P. 703-710. - doi: 10.1007/978-3-642-12535-5_84.

6. Шишкин Г.И. Сеточные аппроксимации сингулярно возмущенных эллиптических и параболических уравнений. - Екатеринбург: УрО РАН, 1992. - 232 с.

7. Андреев В.Б. О точности сеточных аппроксимаций негладких решений сингулярно возмущенного уравнения реакции-диффузии в квадрате // Дифференц. уравнения. -2006. - Т. 42, № 7. - С. 954-966.

8. Han H., Kellogg R.B. Differentiability properties of solutions of the equation —e2 Au + + ru = f (x, y) in a square // SIAM J. Math. Anal. - 1990. - V. 21, No 2. - P. 394-408. -doi:10.1137/0521022.

9. Бахвалов Н.С. К оптимизации методов решения краевых задач при наличии пограничного слоя // Журн. вычисл. матем. и матем. физики. - 1969. - Т. 9, № 4. -С. 841-859. - doi: 10.1016/0041-5553(69)90038-X.

10. Ершова Т.Я. О решении задачи Дирихле для сингулярно возмущенного уравнения реакции-диффузии в квадрате на сетке Бахвалова // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и кибернетика. - 2009. - № 4. - С. 7-14. - doi: 10.3103/S0278641909040013.

11. Axelsson O., Layton W. A two-level discretization of nonlinear boundary value problems // SIAM J. Numer. Anal. - 1996. - V. 33 - P. 2359-2374. - doi: 10.1137/S0036142993247104.

12. Xu J. A novel two-grid method for semilinear elliptic equation // SIAM J. Sci. Comput. -1994. - V. 15. - P. 231-237. - doi: 10.1137/0915016.

13. Vulkov L.G., Zadorin A.I. Two-Grid Algorithms for the Solution of 2D Semilinear Singularly perturbed convection-diffusion equations using an exponential finite difference scheme // AIP Conf. Proc. - 2009. - V. 1186. - P. 371-379. - doi: 10.1063/1.3265351.

14. Zadorin A.I., Tikhovskaya S.V., Zadorin N.A. A two-grid method for elliptic problem with boundary layers // Appl. Numer. Math. - 2015. - V. 93. - P. 270-278. - doi: 10.1016/j.apnum.2014.06.003.

15. Тиховская С.В. Двухсеточный метод для эллиптического уравнения с пограничными слоями на сетке Шишкина // Учен. зап. Казан. ун-та. Серия Физ.-матем. науки. -2012. - Т. 154, кн. 4. - С. 49-56.

16. Zadorin A.I., Zadorin N.A. Interpolation formula for functions with a boundary layer component and its application to derivatives calculation // Сиб. электр. матем. изв. -2012. - Т. 9. - С. 1-11.

17. Шишкин Г.И., Шишкина Л.П. Метод Ричардсона высокого порядка точности для квазилинейного сингулярно возмущенного эллиптического уравнения реакции-диффузии // Дифференц. уравнения. - 2005. - Т. 41, № 7 - С. 980-989.

18. Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. - М.: Наука, 1984. - 320 с.

19. Ильин В.П. Методы конечных разностей и конечных объемов для эллиптических уравнений. - Новосибирск: Изд-во ИВМиМГ СО РАН, 2001. - 318 с.

20. Шишкин Г.И., Шишкина Л.П. Улучшенная разностная схема метода декомпозиции решения для сингулярно возмущенного уравнения реакции-диффузии // Труды ИММ УрО РАН. - 2010. - Т. 16, № 1. - C. 255-271.

74

С.В. ТИХОВСКАЯ

21. Задорин А.И., Задорин Н.А. Интерполяция функций с погранслойными составляющими и ее применение в двухсеточном методе // Сиб. электр. матем. изв. - 2011. -Т. 8. - С. 247-267.

22. Ильин В.П. Методы и технологии конечных элементов. - Новосибирск: Изд-во ИВ-МиМГ СО РАН, 2007. - 371 с.

23. Saad Y. Iterative Methods for Sparse Linear Systems. - N. Y.: PWS Publ. Co, 1996. -448 p.

24. Ильин В.П. Методы бисопряженных направлений в подпространствах Крылова // Сиб. журн. индустр. матем. - 2008. - Т. 11, № 4 (36). - С. 47-60. - doi: 10.1134/ S1990478910010102.

25. Тиховская С.В. Разработка разностных схем на сгущающихся сетках для краевых задач с пограничным слоем: Дис. ... канд. физ.-матем. наук. - Омск, 2013. - 105 с.

Поступила в редакцию 23.01.15

Тиховская Светлана Валерьевна - кандидат физико-математических наук, научный сотрудник, Институт математики им. С.Л. Соболева СО РАН, Омский филиал, г. Омск, Россия.

E-mail: s.tihovskaya@yandex.ru