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

www.volsu.ru

DOI: http://dx.doi.org/10.15688/jvolsu1.2015.3.4

УДК 519.62 ББК 22.19

РАСШИРЕННОЕ СЕМЕЙСТВО ДВАЖДЫ НЕЯВНЫХ МЕТОДОВ ДЛЯ ЖЕСТКИХ СИСТЕМ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

Евгений Иванович Васильев

Доктор физико-математических наук,

профессор кафедры фундаментальной информатики и оптимального управления, Волгоградский государственный университет [email protected], [email protected]

просп. Университетский, 100, 400062 г Волгоград, Российская Федерация

Татьяна Анатольевна Васильева

Кандидат физико-математических наук,

доцент кафедры фундаментальной информатики и оптимального управления, Волгоградский государственный университет [email protected], [email protected]

просп. Университетский, 100, 400062 г. Волгоград, Российская Федерация

Мария Николаевна Киселева

Аспирант,

кафедра фундаментальной информатики и оптимального управления, Волгоградский государственный университет [email protected], [email protected]

просп. Университетский, 100, 400062 г. Волгоград, Российская Федерация

Аннотация. Проведено исследование свойств устойчивости расширенного трех-параметрического семейства дважды неявных разностных схем со второй производной. Показано, что семейство имеет 5-й порядок точности. Также показано, что среди множества Л-устойчивых 2ISD-схем существует два двухпараметрических семейства: семейство ¿-устойчивых схем и семейство схем повышенной точности для линейных задач. Проведено тестирование таких разностных схем на задачах с различной степенью жесткости. Представлены зависимости интегральной погрешности численного решения от величины шага интегрирования, на основе которых проведен анализ свойств устойчивости и точности предлагаемых 2ISD-схем.

Ключевые слова: ¿-устойчивость, Л-устойчивость, жесткие системы, неявные методы, мульти-неявные методы, методы со второй производной.

Введение

В приложениях численных методов для жестких систем центральной является проблема © устойчивости методов. Требование абсолютной устойчивости неизбежно приводит к неявным

о <м

л «

<и

ч

<и

Н'

л «

<и

«

о л

т

м «

<и Л

разностным методам. Однако для очень жестких задач условие A-устойчивости может оказаться недостаточным. Для них требуются Z-устойчивые методы, простейшим примером которых является неявный метод Эйлера первого порядка.

В настоящей работе рассматривается новое семейство разностных схем для жестких систем дифференциальных уравнений, которое является расширением мульти-неявных схем со второй производной, представленных в публикациях [1; 2; 4].

В классических многошаговых методах используется информация в нескольких предшествующих точках. В противоположность этому в данной работе рассматриваются методы, которые используют несколько последующих точек, для каждой из которых записывается отдельное разностное уравнение. Такой подход приводит к построению семейства так называемых мульти-неявных методов со второй производной (m-Implicit Second Derivative). mISD-схе-ма представляет собой систему разностных уравнений. Порядок аппроксимации отдельных уравнений можно варьировать с сохранением порядка схемы в целом, что позволяет расширить семейство mISD-схем.

В работе [2] рассмотрено обобщение SISD-схем (3-Implicit with Second Derivative) 8-го порядка, проведенное за счет снижения порядка аппроксимации промежуточных разностных уравнений. Такой подход позволяет ввести свободные параметры, за счет которых можно обеспечить улучшенные свойства устойчивости. В данной работе аналогичная процедура обобщения выполнена для 2ISD-схем 6-го порядка [1; 4]. Построено трехпараметрическое семейство 2ISD-схем. Показано, что среди множества A-устойчивых схем существуют семейство Z-устойчивых схем и семейство схем повышенной точности для линейных задач. Проведено тестирование таких разностных схем на различных задачах с различной степенью жесткости. Представленные зависимости интегральной погрешности численного решения от величины шага интегрирования хорошо демонстрируют качество устойчивости и точности предлагаемых дважды неявных схем.

1. Мульти-неявные методы со второй производной

Рассмотрим задачу Коши для автономной системы обыкновенных дифференциальных уравнений

d-uit) = f (u), u(0) = Uo, t > 0, (1)

где u (Uj(t), u2(t), ... , up(t)}T и f(u) являются векторными функциями размерности р. Введем сетку {t = in, n = 0, 1, 2, ...} с постоянным шагом т > 0. Функции, определенные на сетке, и приближенные решения обозначим как vn = v(tn), fn = f (vn).

Основная идея мульти-неявных методов для задачи (1) заключается в совместном использовании m разностных уравнений, включающих в неявном виде m искомых значений численного решения

V"+k _ = £(akE + $,kiJn+i)fn+i, k = 1, 2, ..., m, (2)

где vt и f - р-векторы; E и J - единичная матрица и матрица Якоби размера p х p:

f = f (х), J = f (х), i = n, n+1, ..., n+m.

ou

При известном vn = v(tn) разностная схема (2) представляет собой систему взаимосвязанных нелинейных разностных уравнений относительно неизвестных vn+1, vn+2, ..., vn+m в m последующих точках tn+1, tn+2, ..., tn+m . Такую схему будем называть m-неявной.

Для нахождения коэффициентов разностной схемы (2) запишем ошибку аппроксимации, используя, что для решения системы (1) и выполнено u" = fuu' = f f:

u _u m

шП =_ n+k u+k1 +£(atiu'n+i + фкк+,), k = 1,2, ...,m. (3)

ф i=0

Разложение функций, первых и вторых производных в формуле (3) в ряд по степеням т с центром в точке tn позволяет выписать систему линейных уравнений на коэффициенты, которые обеспечивают

ш,) = О^2), к = 1,2, ..., т.

В работах [1; 2; 4] приведены и исследованы схемы 6-го, 8-го и 10-го порядка аппроксимации.

Разностная схема называется абсолютно устойчивой, если численное решение модельного уравнения

^и^) = ли, л, и(}) е С, Re л < 0, (4)

асимптотически стремится к нулю с ростом t при любой величине шага интегрирования т. Применяя схему (2) к уравнению (4), получим однородную систему линейных уравнений (5) относи-

тельно V , V ,, ..., V .

п п+1 ' п+т

^^п^к — ^+к-1 =Х (аи2 + , к =1,2, ..., т. (5)

I=0

Коэффициенты системы и, следовательно, само решение системы зависят от г = А,т.

Устойчивость рассматриваемой схемы разумно идентифицировать по крайним значениям решения системы (5). Функцией роста R (г) для системы линейных уравнений (5) назовем отношение

V

Rm (г) = . (6)

п

Область в комплексной плоскости, в которой модуль функции роста меньше единицы, называется областью устойчивости метода:

г е5 » ^т(г)|< 1. (7)

В работе [2] показано, что все mISD-схемы максимального порядка (то есть 2т + 2) являются А-устойчивыми, причем область устойчивости в точности совпадает с левой полуплоскостью комплексной плоскости.

2. Трехпараметрическое семейство 2ISD-схем 5-го порядка

В работе [2] рассмотрено обобщение 3ISD-схем 8-го порядка за счет снижения порядка аппроксимации двух первых разностных уравнений. Такой подход сохраняет суммарный 8-й порядок аппроксимации, но в то же время позволяет ввести два свободных параметра, за счет которых можно обеспечить улучшенные свойства устойчивости. Аналогичный подход для 2ISD-схем не приводит к успеху. Поэтому расширим количество свободных параметров за счет снижения суммарного порядка аппроксимации с 6-го до 5-го.

Запишем 2ISD-схему в ином виде, который может быть получен путем линейной комбинации отдельных разностных уравнений:

V , — V

п+1 п = (0,0 Е + #10 0п ) /п + (апЕ + фп Уп+1) 1п+1 + К Е + ф^ 2) /п+2

ф

V о — V (8)

"+* п = (020Е + #20 Jn ) 1п + (021Е + ф^п+О /п+1 + (022Е + <ф>22 Л+2) /п+ 2 2ф

Коэффициенты для метода, обеспечивающие 6-й порядок аппроксимации всех входящих уравнений, в этом несимметричном виде были получены в работе [4]:

(% ) =

( 101 128 11 Л

240 240 240

56 128 56

V 240 240 240 )

Ъ ) =

( 13 40 3 Л

240 240 240

8 V 240 0 8 240 )

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

( % ) =

^ + 3б - 2в

-56- - 3г

V 240

128 + 4в

240 + 4в

128 240

^40 -36 - 2в

140 + 3г

Л

( 13

( ъи ) =

20 + б - в

V 240

- г

-Ж. + 4б

240 ^

-4г

- 240+б+в

-----г

240

(9)

В частном случае а = р = у = 0 все разностные уравнения в (9) имеют 6-й порядок аппроксимации.

Для аналитического вычисления функции роста схемы (8)-(9) требуется решить систему, коэффициентами которой являются квадратичные многочлены с неопределенными коэффициентами. Объем выкладок очень большой, и для их выполнения использовались пакеты символьных вычислений. Функция роста будет иметь следующий вид:

R( г) =

Р(г) _ 1 + р1г + р2г2 + р3г3 + р4г4 Q(г) 1 + q1г + q2г2 + q3г3 + q4г4 '

где

р1 = 1 - 4в - 6г,

Р2 = 30-|б - 4в -34 г + 24гв,

Р3 = 110 -1б -?в -§г + 24гв,

Р4 = ^ - Т1б - 1г + 8гв, Для А -устойчивости требуется, чтобы

дх = -1 - 4в - 6г,

д2 = -30 - |б + 4в + ^г + 24гв,

=-110 + Тб-^в - 17г -24гв, = 1- Т1б + 30г + 8гв.

|Р( г )| IQ( г )|

< 1 для всех г: Re г < 0.

(10)

Для этого достаточно, чтобы аналогичное неравенство выполнялось для мнимой оси, то есть

|0Оу)|2-|Р(/>)|2 > 0, Уу, (11)

и чтобы все корни многочлена 0(г) лежали в правой полуплоскости. Так как в нашем случае условие для корней выполнено при а = р = у = 0 и корни непрерывно зависят от параметров, то неравенство (11) будет описывать все множество А-устойчивых схем. Выражение для разности квадратов знаменателя имеет вид (выкладки опускаем)

| Q(iy) |2 -1 РЦу) |2 = -1 у [Цв -| г + (15 -|б -8г + 48гв)гу2].

Отсюда видно, что неравенство (11) выполняется при параметрах, удовлетворяющих неравенствам

г > 0

в > 1г (12) б < £ - г(1 - 30в)

Заметим, что первое неравенство следует и из более простых соображений, так как для А -устойчивости необходимо (0 < p4 < q4), что приводит к следствию у> 0. Семейство ¿1-устой-чивых разностных схем получим, если в дополнение к неравенствам (12) потребуем выполнение условия p4 = 0. В результате получим следующее соотношение:

6 = 114 - г(т - 30в). (13)

В трехмерном пространстве (а, р, у) неравенства (12) образуют область, сечения которой плоскостями у = const (для значений г : %5, 14315, 6з15 и 0 ) изображены на рисунке 1. Внутренность области А-устойчивости закрашена серым цветом. Внутри нее изображено семейство параметров (13), при которых схема становится L1-устойчивой. Заметим, что для всех сечений боковая граница области А-устойчивости и семейство L1-устойчивости являются параллельными прямыми, всегда проходящими через фиксированные точки C (б = 124,в = %0) и D (б = У24, в = %„) соответственно. Как видно на рисунка 1, г, при у = 0 обе прямые сливаются в одну. В диапазоне 0 < г < 14^15 на луче L1-устойчивости имеется точка, в которой в которой p3 = 0:

б = 8 - }в - Kf - 30в). (14)

При этих параметрах получаем L2-устойчивую разностную схему. На рисунке 1, б и в, она отмечена буквой F.

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

= exp(2z),

Piz) Q (z).

Их разность будет характеризовать погрешность метода для модельной задачи. Также, опуская выкладки, приведем аналитическое выражение для этой разности

ч Р(г) С6г6 + С7г7 + С8г еХР(2г) — Ш - 6 150(г) ' ^ (15)

где

С6 = г — fв, С7 = 315 — Т5б — |в + 15г + 16^ С8 = 31? — Т5б — ^Цв + Цг + 16гв.

Как и ожидалось, в общем случае внутри области А-устойчивости схема имеет 5-й порядок точности, так как правая часть в (15) является О(г6). Однако на нижней границе схема имеет 6-й порядок точности. Помимо этого, на нижней границе существуют точки с 7-м и 8-м порядком аппроксимации. Для 8-го порядка существует единственная точка Е (б = }Г68, в = 0, г = 0), которая помечена на рисунке 1, г. Там же помечена точка О (а = р = у = 0), соответствующая базовой схеме 6-го порядка. Для 7-го порядка существует целое семейство. На рисунке 1, в, помече-

un+ 2

Vn+2

Un

vn

на точка G (б = -535880, в = }Г48, г = 6315), которая лежит на пересечении этого семейства с семейством ¿^устойчивости. На рисунке 1, б, помечена точка F(б = -23360,в = 160, г = 14315), которая лежит на пересечении семейства 6-го порядка, то есть нижней границы области А-устойчивости, с семейством ¿2-устойчивых схем.

Р 0.09 и/ /Ь

0.06

0.03 //<* ч:

о/

-0.09 -0.06 -0.03 0.03 а

У= 6/315 -0.03

Рис. 1. Диаграмма четырех сечений области А-устойчивости 2ISD-схемы в пространстве свободных параметров (а, р,у)

в

г

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

Таким образом, для последующего тестирования отберем из семейства (8)-(9) четыре разностные схемы, которые соответствуют точкам О, Е, G и F:

1) А(6) - базовая А-устойчивая схема 6-го порядка (а = р = у = 0);

2) А(8) - А-устойчивая схема 8-го порядка (б = )168, в = 0, г = 0);

3) L1(7) - ¿^устойчивая схема 7-го порядка (б = -535880,в = 1148, г = 2Го5);

4) L2(6) - ¿2-устойчивая схема 6-го порядка (б = -23360, в = г = %5).

3. Результаты тестирования 2ISD-схем

Тестирование проводилось на трех тестовых задачах, которые перечислены ниже. Тест 1. Линейная система с постоянными коэффициентами. Первая стадия тестирования проводилась на примере численного решения линейной системы (16) с тремя переменными, для которой несложно выписать точное решение.

—u(t) =

Г-2 9 -1^ -8 -3 1

1

2 -12

u (^ на t е [0,1] с начальными условиями u(0) =

ГП

V1/

(16)

Матрица системы имеет как вещественные, так и комплексные собственные числа.

Тест 2. Жесткая нелинейная система без пограничного слоя. Дальнейшее тестирование проводилось на частном случае тестовой задачи Капса [3, с. 238], которая представляет собой жесткую (при p >> 1) автономную систему двух нелинейных дифференциальных уравнений

ц'( t) = -(р + 2)п + рп22

2 , t е[0,2] , u2 (t) = u1 - u2 - П2

(17)

с начальными условиями нДО) = 1, и2(0) = 1. При этих условиях реализуется плавное решение щ(() = е-,ы2(() = е~'/2, не зависящее отp (см. рис. 2, а).

Тест 3. Жесткая нелинейная система с пограничным слоем. Та же система Капса, что и в предыдущем случае, но при иных начальных условиях мх(0) = 0, и2(0) = 1. В этом случае реализуется решение с переходной зоной (пограничный слой, см. рис. 2, б), толщина которой 5 « 4/р. График для решения (рис. 2, б) имеет две типичные зоны поведения решения, характерные для жесткой компоненты. Первая зона включает быстрое изменение функции п^), когда она переходит к равновесному значению. Затем эта функция медленно меняется вблизи равновесного значения. Поведение второй компоненты мало зависит от параметра р. Она относится к так называемым медленным компонентам.

б

Ш)

\и2(0

11(1)

Рис. 2. Решение тестовой задачи (17) с параметром р = 100 при различных начальных условиях:

а - плавное решение; б - решение с пограничным слоем

1

а

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

е = || и(Т) - и '(Т )|| (18)

II и*(Т)|| '

На рисунке 3 представлена зависимость в (0 в логарифмических осях для четырех вышеупомянутых схем. Контрольное решение предварительно рассчитывалось с максимально возможной точностью. Для лучшей идентификации кривых дополнительно черными кружками отмечены результаты для базовой схемы А(6).

ш А(8) /

Li(7) / / L2(6)

/ // / У А(6)

г

Рис. 3. Зависимость погрешности в численного решения от величины шага интегрирования т для различных вариантов 2ISD-схем в линейной задаче (16)

По углам наклона кривых в логарифмических осях хорошо видно, что порядки точности всех схем соответствуют теоретическим за исключением начального и конечного участков. Хаотическое поведение погрешности в районе в = 10-15-10-16 объясняется пределом машинной точности. Данные расчеты проводились со стандартной двойной точностью формата вещественных чисел с мантиссой 52 бит.

На рисунке 4 в логарифмических осях представлены зависимости погрешностей (18) от величины шага интегрирования для нежесткого (р = 1) и жесткого (р = 104) случаев на плавном решении (рис. 2, а). По рисунку 4, а, видно, что на нелинейной задаче преимущество схемы А(8) пропадает, она работает с 6-м порядком точности. ¿-устойчивые схемы демонстрируют 5-й порядок точности.

На рисунке 4, б, представлены аналогичные результаты для жесткого случая. Видно, что качество расчета не ухудшилось, все схемы устойчивы, порядок точности у всех схем близок к шести. Проявляется преимущество схем А(8) и L1(7). Это связано с тем, что при больших значениях параметра р система приближается к линейной системе, и при этом начинает отчетливо проявляться преимущество схем повышенного порядка.

Следующее тестирование проведено для решения с пограничным слоем (рис. 2, б). На рисунке 5 также в логарифмических осях представлены зависимости погрешностей от шага интегрирования для жесткого случая с р = 103 и р = 104. Видно, что поведение погрешности сильно меняется по сравнению с предыдущими тестами. Отчетливо видны недостатки строго А-устойчивых схем, которые плохо рассчитывают переходный процесс при большом числе жесткости. Для ¿-устойчивых схем таких проблем не возникает, в этом и заключается их главное преимущество. При данном числе жесткости преимущество ¿2 над ¿1 почти не видно, оно проявляется при большей жесткости. Почти во всем диапазоне графика на рисунке 5 шаг по времени т > 5, 5 - размер пограничного слоя. В ре-

ПРИКЛАДНАЯ МАТЕМАТИКА ^^^^^^^^^^^^^^^^^^^^^^^

зультате погрешность крайне неравномерно распределена по отрезку интегрирования. Суммарная погрешность складывается из двух частей: в1 - погрешность на пограничном слое и в2 - погрешность на всем остальном участке интегрирования. Для методов А(6) и А(8) в1 очень большая. В этом и заключается главный их недостаток по сравнению с ¿-устойчивыми методами. Существенное отличие Л-устойчивых методов от ¿-устойчивых исчезает лишь при т < 105. Немонотонное поведение погрешностей для ¿-устойчивых схем связано с немонотонным поведением разности функций роста и экспоненты в выражении (15) при больших отрицательных значениях переменной г.

а б

Рис. 4. Зависимость погрешности е от величины шага интегрирования т для различных вариантов 2ISD-схем на плавном решении теста Капса (17):

а) р = 1; б) р = 104

а б

Рис. 5. Зависимость погрешности е численного решения от величины шага т для различных вариантов 2ISD-схем для решения с пограничным слоем, тест № 3:

а) р = 103; б) р = 104

Заключение

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

СПИСОК ЛИТЕРА ТУРЫ

1. Васильев, Е. И. Мульти-неявные методы со второй производной для жестких систем дифференциальных уравнений / Е. И. Васильев, Т. А. Васильева, М. Н. Киселева // Вестник Волгоградского государственного университета. Серия 1, Математика. Физика. - 2012. - №> 2 (17). - С. 68-77. - DOI: 10.15688/jvolsu1.2012.2.8.

2. Васильев, Е. И. L-устойчивость мульти-неявных методов 8-го порядка для жестких систем дифференциальных уравнений / Е. И. Васильев, Т. А. Васильева, М. Н. Киселева // Вестник Волгоградского государственного университета. Серия 1, Математика. Физика. - 2013. - №> 1 (18). - С. 70-83. - DOI: 10.15688/jvolsu1.2013.1.6.

3. Деккер, К. Устойчивость методов Рунге - Кутты для жестких нелинейных дифференциальных уравнений / К. Деккер, Я. Вервер. - М. : Мир, 1988. - 334 с.

4. Vasilev, E. High Order Implicit Method for ODEs Stiff Systems / T. Vasilyeva, E. Vasilev // Korean Journal of Computational & Applied Mathematics. - 2001. - Vol. 8, №№ 1. - P. 165-180.

REFERENCES

1. Vasilyev E.I., Vasilyeva T., Kiseleva M.N. Multi-neyavnye metody so vtoroy proizvodnoy dljya zhestkikh sistem differentsialnykh uravneniy [Multi-Implicit SD-Methods for Stiff Systems of Differential Equations]. Vestnik Volgogradskogo gosudarstvennogo universiteta. Seriya 1, Matematika. Fizika [Science Journal of Volgograd State University. Mathematics. Physics], 2012. no. 2 (17), pp. 68-77. DOI: 10.15688/jvolsu1.2012.2.8.

2. Vasilyev E.I., Vasilyeva T.A., Kiseleva M.N. L-ustoychivost multi-neyavnykh metodov 8-go poryadka dlya zhestkikh sistem differentsialnykh uravneniy [L-Stability of Multi-Implicit Methods of 8th Order Differential Stiff Systems]. Vestnik Volgogradskogo gosudarstvennogo universiteta. Seriya 1, Matematika. Fizika [Science Journal of Volgograd State University. Mathematics. Physics], 2013, no. 1 (18), pp. 70-83. DOI: 10.15688/jvolsu1.2013.1.6.

3. Dekker K., Verwer J.G. Ustoychivost metodov Runge-Kutty dlya zhestkikh nelineynykh differentsialnykh uravneniy [Stability ofRunge-Kutta Methods for Stiff Nonlinear Differential Equations]. Moscow, Mir Publ., 1988. 334 p.

4. Vasilyeva T.A., Vasilyev E.I. High Order Implicit Method for ODEs Stiff Systems. Korean Journal of Computational & Applied Mathematics, 2001, vol. 8, no. 1, pp. 165-180.

THE EXTENDED FAMILY OF 2ISD-METHODS FOR DIFFERENTIAL STIFF SYSTEMS

Evgeniy Ivanovich Vasilyev

Doctor of Physical and Mathematical Sciences, Professor, Department of Fundamental Informatics and Optimal Control, Volgograd State University [email protected], [email protected]

Prosp. Universitetsky, 100, 400062 Volgograd, Russian Federation

Tatyana Anatolyevna Vasilyeva

Candidate of Physical and Mathematical Sciences, Associate Professor, Department of Fundamental Informatics and Optimal Control, Volgograd State University [email protected], [email protected]

Prosp. Universitetsky, 100, 400062 Volgograd, Russian Federation

Mariya Nikolaevna Kiseleva

Postgraduate Student,

Department of Fundamental Informatics and Optimal Control, Volgograd State University [email protected], [email protected] Prosp. Universitetsky, 100, 400062 Volgograd, Russian Federation

Abstract. The new set of absolutely stable difference schemes for a numerical solution of ODEs stiff systems (1) is submitted:

du(t) = f (u), dt

t > 0, u(0) = u0.

(1)

The main feature of the set is the multi-implicit finite differences with the second derivatives of the desired solution. The expanded three-parameter (a, p, y) set of 2ISD-schemes (2)-(3) is studied in more details in this paper.

ф

- = X (aliE + «U) fn+

v„^ - v„

i=0 2

2ф

= X (a2lE + ««J) fn+l

(2)

( aki ) =

^ + 3б - 2в - 3г

V 240 i 13

(PU ) =

# + б - в

V 240

128 + 4в

240 + 4B

128 240

- + 4б

240 U

-4г

150 - 36 - 2в

140 + 3г . - 240 + б + ^

(3)

At arbitrary (a, p, y) parameters last difference equation in system (2) has 5 th order of accuracy.

We found that the set of absolutely stable 2ISD-schemes includes two families: the set of the L-stable schemes and the set of the schemes of heightened accuracy for linear problems. For example:

at 6 = Xes, b = 0, r = 0 we have A-stable scheme with 8th order of approximation,

at 6 = -%80,b = X48,r = 6315 we have L 1-stable scheme with 7th order of approximation,

at 6 = - 2;336(), b = %o, r = 143i5 we have L2-stable scheme with 6th order of approximation.

The testing of this difference schemes on linear and nonlinear problems with a different stiff power is conducted. The errors of a numerical solution as functions of integration step size are computed in numerical experiments. These results demonstrate high quality of stability and accuracy of the suggested 2ISD-schemes.

Key words: L-stability, A-stability, stiff systems, implicit methods, multi-implicit methods, methods with second derivative.

vv

n+1 n

-----г

240