Адаптивные методы решения сингулярно возмущенных задач при использовании неоднородных разностных схем

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

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

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

Постановка задачи. В области Q = (0,1) х (0,1) с границей Г рассмотрим следующую сингулярно возмущенную задачу с граничным условием Дирихле:

EHd^u Idx² + d²u Idy²)+ L _Q u=f(x,y), если(x,y)е Q. (1)

u(x,y) = G (x,y), если (x,y) е Г. (2)

В зависимости от вида оператора L0 возникают различные ситуации, которые исследуются в данной статье с применением компактных схем. Рассмотрим сначала простейшие ситуации, относящиеся к обыкновенным дифференциальным уравнениям, а затем обобщим их на двумерный случай.

Самосопряженный оператор. Исследуем случай, когда оператор L₀u и граничные условия представлены в следующем виде:

L o u = d ( x ), n(G)=g₀, n(1)=g_r (3)

Относительно коэффициентов задачи предполагается необходимая гладкость, так что решение u(х) обладает достаточным количеством производных, необходимых для обсуждения разностных схем. Кроме того,

d ( x ) >  d ₀ ² > 0. (4)

Хорошо известно, что структуру решения можно представить в виде разложения гладких функций v0 + е2v1 +... и функций типа пограничного слоя в окрестности точки х = 0 с помощью ряда r0 + е2r1 + .... Для функций типа пограничного слоя доказаны оценки на рост производных в окрестности точки х = 0:

d^krJdx^k\ < Cе^kexp(aр Iе), еслиk=0,1,2,..., а>0, (5) с константой С, не зависящей от малого параметра е. Аналогичные представления существуют в окрестности точки х = 1в виде ряда q ₀ + е²q ₁ +.с ростом производных:

\d^kq.I dx^k\k exp(a р I е),еслик = 0,1,2, .,a>0. (6)

Для построения и анализа схем исходная область должна быть представлена в виде декомпозиции подобластей Q = Q₀ и Q j и Q . Соответственно подобласти имеют вид Q ₀ = (0, о ], Q ₀₁ = [ о ,1 -о ], Q 1 = [1 -о ,1].

В математической литературе обсуждается выбор параметра е, который связан с малым параметром е и с характерным параметром N, отождествляющимся с количеством узлов разностной схемы и порядком аппроксимации k. Обычно принимают о = min [113, k е In (N) I d₀].

Рассмотрим схемы повышенного порядка точности, которые при е = 1 хорошо зарекомендовали себя для достаточно гладких функций в области Q = [о; 1 - о]. Введем множество внутренних четных узлов разностной сетки при N = 2 х n четном: Щ_чет = { х , = i h, 1 = 2, 4, 6,., N-2, h = 11 А} и множество внутренних нечетных узлов разностной сетки Ю_нечет = { х , = ih, 1= 1, 3, 5, ...,N-1, h = 1 IN}. Тогда разностная сетка щ = Щ_чет и Ю_нечет. Используя центральную разность, запишем

Lhv = -е2vh о (х,) + d (x,)v(x,), x е ®нечет, xx

L²hv = -е ² v ^{2 h} о ( x, ) + d ( x, ) v ( x, ), x , g ю _чет .

xx

V 0 0 ( Xi ) = ( V i + 1 - 2 V i + V i -1 )/ h h xx

V ² h 0 ( Xi ) = 0,25 ( Vi ₊ 2 - 2 Vi + Vi _-2)/ h ² . xx

В этих обозначениях рассмотрим разностную задачу:

Lu = f ⁽ X ), X ^е m .. , u h (0) = g 0 , u h (1) = g i , (7)

L^hu^h - 0,5 L^hu ^{2 h} = 0,5 f ( x i ), . х еш .      (8)

В итоге приходим к системе линейных алгебраических уравнений с пятидиагональной матрицей^ * xU^F * :

fbl a 0 0 0 0 ..... 0 1 a e2 a s 0 0 ...... 0 0 a b3 a 0 0 . ..... 0 Ah = 0 s a e4 a s ...... 0 , (9) 0 0 0 0 0 0 ...... 0 0 0 0 .....s a eN-2 0 0 \ 0 0 0 0 ..... a bN-1 , F* =[A«g0, f2/2- CM2 "^Д V1«£1], (10) а = е2*-2, ^ = 0,125е2*2, Ъ. = 2е2/ *2 - d., е. = 1,75 е2* 2 +d. /2.

В областях [0, о] и [1 - о, 1] шаг сетки * = о / N, а в области [о, 1- о] этот шаг Н= 1-2о / Ни их отношение

С = А/Я=[ке ln(N)]/[d0-2^e 1n(V)].    (11)

Рассмотрим в узле s пятиточечный шаблон, который можно записать в виде

“ -2 ^U 2 ⁺ « 1 ^U 1 ⁺ а ^{U +} “Д ⁺ а ^U ^ ( ¹² ) где значение U₀ соответствует узлу о. Учитывая формулу (11), выберем коэффициенты а . , пользуясь разложением вряд Тейлора:

а2= [е2(3С-2)] / [2С2Я2(2С2+3С+1)], а. |4е2(13С)| / [С2Я2(С2+3С+2)], а2= [Се2(3-2С)] / [2Я2(С2+ ЗС+2)], а1-[Се2(С-3)]/[Я2(2С2+3С+1)],      (13)

а d.a . а. а . “ 1 .             (14)

Рассмотрение узла 1- о проводится аналогичным образом. При этом коэффициенты а1 отображаются симметричным образом:

а2= [Се2(3-2С)] / [2Я2(С2+3С+2)], а,=[4Се2(С-3)] / [Н2(2С2+ 3С+1)], а2= [е2(3С2)] / [2С2Я2(2С2+3С+1)], а1-[С2е2(1-3С)]/[Я2(С2+3С + 2)],(15)

Ла^-а^-а,-^.(16)

Теорема 1. Пусть для задачи (3) выполнены условия (4), тогда приближенное решение задачи (9), (10) сходится к точному решению с порядком к = 4

II u-u*L и*<С(1п(Я)/Я)4,(17)

||u -и * || ^ _ш * = тах |u *. -u(x . )|, если x. е ш * .

Доказательство. В подобласти [о; 1-о]решение не содержит функций типа пограничного слоя, поэтому, согласно [1], справедлива оценка

||u-u *Lш * < СЯ⁴.                (18)

Доказательство этой теоремы приведено в работе [1] для гладких функций. Далее достаточно рассмотреть отрезки Q₀ и Q₁ независимо друг от друга с более мелким шагом, т.е. * = о / N= ке 1n (V) < е. В этой области погрешность аппроксимации зависит от шестых производных точного решения и в силу малого значения шага разностной сетки *«е может быть оценена следующим образом:

*⁴{||е²и < ⁶⁾|| „,ш * + ||и⁽⁴⁾|| „, _И * } < С,1п I V) Vx

Остается сделать оценку в узлах о, 1- о. По разложению в ряд Тейлора следует, что

||u-u* || „, ш* < С {Я3 е31п (V) + *3(е2/ С2)}х x ||u(5)||<СЯеTn(Я)||u(5)||.

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

Lu = е²u^,, + u =j(x), u(0) = u(1) = 1       (19)

с точным решением u(x) = exp((x2 -x) / е).

Обозначим А = max|u *. -u(x . )| и величину скорости сходимости V = 1n (2)^-11n{||u-u * /2|| / ||u-u * ||}. Для Я = 8, к ~ 4 составим соответствующие таблицы (табл. 1,2).

По табл.1,2 следует, что погрешность, независимо от малого параметра, выходит на стационарный уровень при фиксированном числе узлов. Рассмотрим влияние Япри фиксированном значении малого параметра (табл. 3, 4).

При увеличении числа узлов в два раза точность возрастает соответственно порядку сходимости. Значение величины V говорит, что она близка к четырем.

Несамосопряженный оператор. Исследуем случай, когда оператор L₀ и и граничные условия представлены в следующем виде:

L0u = d(x) du/dx, u(0) = g0, u(1)=gp    (20)

Lu = е d²u / dx² + L₀u =Kx).        (21)

Относительно коэффициентов задачи предполагается необходимая гладкость, так что решение u(x) обладает

Таблица 1

Погрешности задачи (19) в области По, Я = 256

е	1,0	10 2	10 2	10 2	10 5	10 6
N	8,1 ■ 10 lu	5,4 ■ 10 9	2,4- 10 6	2,15 ■ 10 6	2,13 ■ 10 6	2,12 ■ 10 6

Таблица 2

Погрешности задачи (19) в области йо, Я = 8

е 0,1 102 102 104 105 106 А 1,0 ■ 103 1,2 ■ 102 1,1 ■ 102 1,1 ■ 102 1,1 ■ 102 1,1 ■ 102 достаточным количеством производных, необходимых для обсуждения разностных схем. Кроме того,

d(^)>d0>0. (22)

Известно, что структуру решения можно представить в виде разложения гладких функций v₀ + е²v₁ +... и функций типа пограничного слоя в окрестности точки х - 0 с помощью ряда г ₀ + е²г ₁ +.... Для функций типа пограничного слоя доказаны аналогичные оценки на рост производных в окрестности точки х -1.

Для построения и анализа схем исходную область представим в виде декомпозиции подобластей Q = Q₀ и Q₁,а подобласти - в виде Q 0 = (0,1 - о ], О , = [1 - о , 1).

Выбор параметра е связан с малым параметром е и характерным параметром N, отождествляющимся с количеством узлов разностной схемы и порядком аппроксимации к. Обычно выбирают о=min [1 / 2, к е In (N) / 8₀].

Рассмотрим схемы четвертого порядка точности, которые при е = 1, хорошо зарекомендовали себя для достаточно гладких функций в области Q₀- [0; 1- о]. Введем множество внутренних четных узлов разностной сетки при N- 2 х и четном:

to = {^ = /*,7= 1,2,3,4,5,6, ^,Я-1,*-1/У}.

В итоге приходим к системе линейных алгебраических уравнений с пятидиагональной матрицей^ * xU * -F^h:

' bl a2 b2 e2 b3 c2 b4 s2 0 0 0 0 ..... 0 ! 0 t3 a3 e3 c3 s3 0 ...... 0 Ah = 0 t4 a4 e4 c4 s4 ...... 0 , (23) 0 0 0 0 0 0 ...... 0 0 0 0 ..... tN-2 aN-2 eN-2 cN-2 0 V 0 0 ...... . p4 p3 p2 pl J F*-[/;- ^Л"^’/^’ • " ,^'n-2 -s^^-2S1,L-1-PйS1], (24) а = -4е / 3*2- 2d./ 3*, с= -4е / 3*2 + 2d./ 3*, s= е /12*2 - d./12*, е. - 2,5е/ *2, t = е/12*2 + d./12*,^0 = -11е /12*2 + 110Д/ 294*, р.    4е 3*2 2d 3*, р2 - -6е / 12*2 -222d./ 294*,р3 - -4е /12*2 + d./294*, р4-е/12*2+12d./294*, Ь0--11е/12*2-110d./294*, b1 - 20е /12*2 -99d./ 294*, b2 - -6е /12*2 + 222d./ 294*, b3 - -4е /12*2 - d./ 294*, b4 - е /12*2 - 12d./ 294*.

В области [1- о, 1] шаг сетки * - о / N, а в области [0,1-о] шаг сетки Н- 1-о/Яи их отношение

С-*/Н-[ке ln(N)]/[d0-ks ln(N)].      (25)

Необходимо исследовать схемы в узлах 1-о-Н, 1- о, 1-о-*.

Рассмотрим в узле 1- о пятиточечный шаблон, который можно записать в виде a-2U-2 + (Х.^ + а U + «1^1 + «2^2-^

Значение U ₀ соответствует узлу 1-о. Учитывая формулу (24), выберем коэффициенты а . , пользуясь разложением в ряд Тейлора:

а2-е(3С-2)/[2С2Я2(2С2+3С+1)]-

-dw/[2CН(2C2+3C1)],(26)

а₁ - 4е(1-3С) / [2С²Н²(С²+ 3С + 2)] +

+ 4 dw/[СЯ(2С2+3С +2)],(27)

а_-2 - Се(3-2С) / [2Н²(С²+ 3С+2)] +

+ C2dw/[2Н(C2+3C+2)],(28)

а-1 - 4Се(С-3) / [Н2(2С2+3С +1)] -

- 4C2dw/[Н(2С2+3С+1)],(29)

а0--а-2-а-1-а2-а1.

Рассмотрим узел 1-о-Н. Значение U0соответствует узлу 1-о-Н. Учитывая формулу (24), выберем коэффициенты а, пользуясь разложением в ряд Тейлора:

а2-2е/[СН2(С3 + 6С2+11С+6)]-

-2dw/[CН(C3+6C2+11C+6)],       (30)

а₁--е(1 + 3С)/[3СН²] + (1+СН/[СН(2С+1)], (31) а_-2 - -е / [3Н²(С + 3) - (С+ 1)d_w/ [Н(С+2)], (32) а_-1 --е(С+3) / [Н²(С+2)] + (С+ 1)d_N/ [6Н(С+3)], (33) а₀--а_-2-а_-1-а₂-а₁.

Рассмотрим узел 1- о + *. Значение U₀ соответствует узлу 1- о + h. Учитывая формулу (24), выберем коэффициенты а . , пользуясь разложением в ряд Тейлора:

а-2-С2е/[Н2(6С3+11С2+6С+1)] +

+ 2СУЛДН(6С3+11С2+6С+1)],      (34)

а_-1--е(3 + С)/[6С²Н²]-(1+C)d_^/[3CН], (35) а₂--е/ [6СЯ(3С+ 1)] -(С+ 1)d „ Z [6СН(3С+ 1)], (36) а₁ --е(3С+ 1)/ [2С²Н²(2С+1)] +

+ (C+1)dw/[CН(2C+1)],         (37)

Таблица 3

Погрешности задачи (19) в области П_о, е = 10 " ⁵

N	8	16	32	64	128	256
A	1,11 ■ 10 2	3,8 ■ 10 j	8,4 ■ 10 4	1,4 ■ 10 4	1,84 ■ 10 5	2,16 ■ 10 6
V	1,546 5	2,17	2,61	2,9	3,16	-

Таблица 4

Погрешности задачи (19) в области П_о, е = 10 " ¹

N 8 16 32 64 128 256 A 1,0 ■ 10j 8,99 ■ 105 6,52 ■ 106 4,4 ■ 107 4,3 ■ 108 5,4 ■ 109 V 3,0 3,3 3,89 3,79 3,7 - а^-а^-а^-а^-а,.

Теорема 2. Пусть для задачи (18), (19) выполнены условия (20), тогда приближенное решение задачи (21) сходится на отрезке Q₀ к точному решению с порядком к = 4:

II и-и^Ц _m ^h 4.               (38)

Приведем численные расчеты для задачи

Lu ^-ги +и ^: =ф(х), и(0) = ехр(-2/е),и(1) = 1 (39) с точным решением и(х) = ехр((х² + х -2) / г). В области Q₀ решение имеет ограниченные производные, поэтому для задачи (23) рассмотрим область Q₁. Обозначим А = тах|и *. - и(х . )| и для N= 8,к = 4 составим соответствующие таблицы (табл. 5,6).

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

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