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

© Васильев Е.И., Васильева Т.А., Киселева М.Н., 2013

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

УДК 519.62 ББК 22.19

^-УСТОЙЧИВОСТЬ МУЛЬТИ-НЕЯВНЫХ МЕТОДОВ 8-го ПОРЯДКА ДЛЯ ЖЕСТКИХ СИСТЕМ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

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

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

профессор кафедры фундаментальной информатики

и оптимального управления

Волгоградского государственного университета

[email protected]

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

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

Кандидат физико-математических наук, доцент кафедры фундаментальной информатики и оптимального управления Волгоградского государственного университета tatiana_vas @mail. ru

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

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

Аспирант кафедры фундаментальной информатики и оптимального управления Волгоградского государственного университета kiselevamaria 15.08@gmail. com

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

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

численного решения от величины шага интегрирования, на основе которых проведен анализ свойств устойчивости и точности предлагаемых SISD-схем.

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

Введение

Жесткие задачи встречаются во многих областях науки и техники, включая кинетику химических реакций, гидродинамику многофазных течений, приложения в биологии и другое. При численном решении таких задач возникает ряд проблем, связанных с устойчивостью, точностью и вычислительной сложностью применяемых методов. Наиболее полные обзоры по численным методам для жестких систем содержатся в [3; 6]. В настоящей работе рассматривается новое семейство разностных схем, которое является расширением семейства мульти-неявных схем со второй производной, представленных в [1; 7].

В приложениях численных методов для жестких систем центральной является проблема устойчивости методов. Требование абсолютной устойчивости неизбежно приводит к неявным разностным методам. Однако и для них существуют ограничения, например, порядковый барьер Далквиста [4; 5]: неявные линейные многошаговые методы выше второго порядка не могут быть Л-устойчивыми. Однако для очень жестких задач условие A-устойчивости может оказаться недостаточным. Для них требуются Z-устойчивые методы, простейшим примером которых является неявный метод Эйлера первого порядка.

Одним из видов численных методов, в которых возможно одновременно совместить высокие требования к устойчивости и порядку аппроксимации, являются неявные методы Рунге -Кутта [2]. Однако ресурсоемкость их реализации очень велика, что связано с необходимостью применения итерационного поиска решения на каждом шаге интегрирования. Для этого обычно применяется метод Ньютона, основные затраты которого заключаются в многократном вычислении матриц Якоби правой части системы уравнений. Если эти матрицы Якоби использовать в разностной схеме, то можно повысить порядок аппроксимации схемы без увеличения вычислительной сложности при ее реализации. Такой подход использован в многошаговых методах со второй производной.

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

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

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

Ми(г) = /(и), и(0) = ио, г > 0, (!)

м

где и = {М](г), и2(г), ... , ир(г)}т и/(и) являются векторными функциями размерности р. Введем сетку {гп = хп, п = 0, 1, 2...} с постоянным шагом х> 0. Функции, определенные на сетке и приближенные решения обозначим как гп = и(г„), /п = /(гп).

Для численного решения задачи (1) предлагается т раз неявный разностный метод

Юп+к - Юп+к-1

= £ (ак.Е + ТЪкг'1п+г )п+г , к = 1,2,...5 т :

(2)

где VI и Л - р-векторы. Е и J - единичная матрица и матрица Якоби размера р х р:

їг = /( V X ^ =|Ц-(V ), (і = П П + І-, П + т).

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

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

¥кп =

■+£ (акги"+г + тЪЫи"+ г X

і=0

к = 1,2,..., т

(3)

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

щ(гп) = 0(х2т+2), к = 1,2,..., т .

Ниже приведены коэффициенты для случаев т = 2, 3 и 4 в матричной форме записи.

При т = 3:

(акг ) =

(акг ) = (101 240 11 ^ 240 128 240 128 240 11 I 240 101 , 240 У ( 13 ъ )=[т V 80 -1 -1 I 6 80 1 -13 6 240

( 6893 313 89 397 I ( 1283 -851 -269 ГО 40 т

18144 672 672 18144 30240 3360 3360 30240

3 109 109 3 , ( Ък г ) = 31 113 -113 -31

224 224 224 224 10080 1120 1120 10080

397 89 313 6893 163 269 851 -1283

V18144 672 672 18144 ^30240 3360 3360 30240 У

(4)

Т

ип+к ип+к-1

Т

При m = 4:

К-) =

4354560 272160 630 272160 4354560 725760 90720 320 18144 725760

26081 122341 313 12091 14111 893 6887 -47 -1721 -103

4354560 14111 272160 12091 630 313 272160 122341 4354560 26081 , (Ьк1) = 725760 103 90720 1721 320 47 90720 -6887 145152 -893

4354560 272160 630 272160 4354560 145152 90720 320 90720 725760

59681 38341 103 89371 1539551 2237 1243 81 31207 -26051

V4354560 272160

272160 4354560 J

V725760 18144

90720 725760 J

Разностные схемы с этими коэффициентами имеют соответственно 6-й, 8-й и 10-й порядок аппроксимации. Две первые были представлены в работах [1; 7]. Далее будем их именовать 2ISD, 3ISD и 4ISD соответственно (m-Implicit Second Derivative).

Видно, что построенные разностные схемы обладают центральной симметрией для коэффициентов aki и антисимметрией для коэффициентов bki, а именно

% = aM,m-i , К = bm+1-k,m-i , к = i = °,-m. (5)

Свойство (5) является следствием структуры разностной схемы. Несложно показать, что

оно приводит к взаимной симметрии областей устойчивости и неустойчивости. Разностная

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

dtu(t) = А и, А, и({) eC, Re А < 0 (6)

асимптотически стремится к нулю с ростом t при любой величине шага интегрирования т.

Применяя схему (2) к уравнению (6), получим однородную систему линейных уравнений (7)

относительно vn, v„+i, ... vn+m.

m

vn+k - Vn+k-1 =Z (akrz + V'K+i k = I2— m . (7)

i=0

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

Устойчивость mISD-схемы разумно идентифицировать по крайним значениям решения системы (7). Функцией роста Rm(z) для (7) назовем отношение

v

Rm (Z) = -V+m . (8)

^ n

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

Z » 1 Rm(Z)|<1. (9)

Несложно доказать, что в силу вещественности и симметрии коэффициентов (5) решение системы (7) обладает следующими свойствами.

Свойство 1. Если при z = ц решением является v = (vn, vn+1,..., vn+m), то при z = ц решением будет w = (vn , vn + 1, ..., vn + m ) ^ след0вательн0, Rm (1) = Rm (l) .

Свойство 2. Если при z = ц решением является v = (vn, vn+1,..., vn+m), то при z = -ц решением будет w = (vn+m, vn+m-1,..., vn) и, следовательно, Rm(-ц) = [Rm(i)]-1.

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

переходит в область неустойчивости. Если при этом у области устойчивости нет других границ, то она совпадает с левой комплексной полуплоскостью, и схема (3) будет являться Л-устойчивой. Вопрос о других границах требует более детального рассмотрения функции роста.

При т = 2 функция роста найдена в [7]

1 + z + 104 z2 + — z3 + — z4

Щ (.г) = 24^ ^10^ 90

При т = 3 функция роста приведена в [1]

к (* ) = 1 + 3 * + ! *2 + 7 *3 +Но *4 + 56) *5 + 560 *6

3 1 -1*+1*2 - 7 *3+163) *4 - 56) *5+560*6 ‘

В обоих приведенных случаях числитель и знаменатель отличаются только знаком при нечетных степенях, что приводит к равенству их модулей при чисто мнимых значениях *. Это подтверждает установленный факт, что мнимая ось является границей области устойчивости. Также видно, что для отрицательных действительных * во всех случаях |Л(*)|<1. Известно, что такая дробно рациональная функция роста обеспечивает Л-устойчивость, если все корни знаменателя лежат в правой полуплоскости.

Для всех mISD-схем это выполняется. В самом деле, комплексные корни знаменателя непрерывно зависят от коэффициентов многочлена, которые, в свою очередь, непрерывно зависят от коэффициентов разностной схемы (2). С помощью непрерывного изменения коэффициентов при сохранении их симметрии (5) любую разностную mISD-схему (2) можно трансформировать в более простую ^-устойчивость схему, например, в совокупность однократно неявных схем 4-го порядка. При непрерывном изменении коэффициентов с сохранением их знаков и свойства симметрии (5) Л-устойчивая схема (2) переходит в Л-устойчивую. Следовательно, mISD-схема также будет А -устойчивой.

2. Двухпараметрическое семейство 3ISD-схем 8-го порядка

Рассмотрим семейство разностных 3ISD-схем, в которых только для последнего значения обеспечивается максимальный 8-й порядок точности, а промежуточные значения вычисляются с порядком на 1 меньше. Для этого перепишем разностную схему (2) применительно к 3ISD-схеме в ином виде, который может быть получен путем линейной комбинации отдельных разностных уравнений в (2):

11" =2>1,Е + ). 1„+,,

L 1=0

*Т)"+2Т Т)" = Е(а2,Е + тЬ2Л+г)!"+,, (10)

2 1=0

3

Т)"+3_Т)" = Е (а3,Е+тМ"+г)/"+,,

'>Т 1=0

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

( 6893 313 89 397 I ( 1283 -851 _269 _163 I

18144 672 672 18144 30240 3360 3360 30240

К,) = 223 1134 10 21 13 42 10 567 , Ь) = 43 1890 _8 105 _19 210 _4 945

31 V 224 81 224 81 224 31 224 19 V 1120 _27 1120 27 1120 _19 1120 у

Порядок аппроксимации последнего разностного уравнения в (10) не изменится, если значения V п+1 и Vп+2 в его правой части будут вычисляться лишь с 7-м порядком. Из этих соображений введем в первое и второе разностные уравнения по одному свободному параметру с одновременным уменьшением их порядка аппроксимации на 1, то есть до 7-го. В итоге получим двухпараметрическое семейство трижды неявных схем со следующими коэффициентами

(акг) =

6893 + 11 а 313 + _8^ _

18144 3 612 672 ^

— 11 18144 3 ^

397

_223 +11 в 10 + 9п ц _ дп _1^ _ 11 в

1134 ^ 3 у 21 ^ 42 ^ 567 3 Г

1134

31

V 224

( 1283 30240

(Ь,) =

+ а

81

224

_851 + 9а

3360 ^ ^

81

224

_269 + 9а

3360 ^ ^

31

224

^ + ал

30240

_8 + 9в + 9В

945 + В

19

V 1120

_27

1120

27

1120

_19

1120

(11)

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

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

R( 2 ) =

Р( 2 ) 1 + р12 + р2 2 2 + р323 + р4 24 + р525 + р6 26

Q(2) 1 _ q12 + q222 _ qъz3 + q424 _ q525 + q626

где

р1 = _ 9а + 18В,

п — 193 27 ^ _|_ 36 В

Р4 = 1680 _ Т а+ Т Р,

П1 = -| + 9а _ 18В, п = ^53 + 18а_ 54 В

4.4 1680 ^ 7 ы 7 Г’

Для ^-устойчивости требуется, чтобы

Р( 2 )|

„ _ 29 99 ™ 180 Я

р2 = 28 _ Та +~Г Р,

Р = 11___243 а + 27 В

У5 560 280 35 ^ ’

п = 29 + 99 а 198 В

П2 = 28 + ~ а_— Р,

п = -Л + 27 а _ 243 В 45 560 70 140 ^’

Р3 = | _ ^9 а + Ш В, р = —---------------а

У6 560 280 “ ’

П3 = 1 + 11а_ Ш В

П6 = 5» _ 120 В,

\б( 2 )|

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

(12)

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

10(,у)|2 _\р(,у)|2 ^ ^ Чу, (13)

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

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

Шу )12 -Р0> )12 = 1580(«- 2Д) у10Уо + («+2Д)+§ у и - («+2Д))].

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

а > 2Д

4 1 .

-^т< «+

135 н 27

(14)

В плоскости (а, Р) область (14) изображена на рисунке 1 с границей ВСАЕ. Внутри области имеется семейство разностных схем, для которых обращается в нуль коэффициент р6 в числителе функции роста (а = ^54 ). Это отрезок АВ на рисунке 1. Он содержит семейство Х1-устойчивых разностных схем. На отрезке АВ имеется точка, в которой р5 = 0 (а = ^54, Р = -1216). При этих параметрах получаем Х2-устойчивую разностную схему. На рисунке 1 она отмечена буквой F.

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

и., . ^+з р( ^)

= ехр(3г),

Q( z)

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

ехр(3г)-

Р(г)

Q( ? )

243 31360 V 270

Уо - (а + 20))г9 + ^(^ - (3а + 40))г

10

Q( г )

(15)

Как и ожидалось, при нулевых параметрах метод имеет 8-й порядок точности, так как правая часть в (15) является 0(г9). С помощью параметров можно повысить порядок точности. Семейство схем 9-го порядка задается условием а + 2 Д = 1/210. Часть этой прямой попадает в область А-устойчивости и на рисунке 1 изображена штриховой линией. На ней примечательными являются две точки. Точка G (а = у54, Д = - Ц35), которая соответствует р1-устойчивой схеме 9-го порядка, и точка D (а = 1540, Д = Х080 ), которая соответствует А-устойчивой схеме 10-го порядка. При этих значениях параметров обнуляется и коэффициент при ^10 в (15).

о Р А Е) і і а 54 27

-2 135 > Ь-устойчивость р 9-ый порядок область А-устойчиеости

Рис. 1. Диаграмма состояний 3ЙВ-схемы в плоскости свободных параметров (а, Д)

и

V

п

п

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

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

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

2) А(10) - А-устойчивая схема 10-го порядка (а = 1/540, Р = 1/1080).

3) Ll(9) - Ll-устойчивaя схема 9-го порядка (а = 1/54, Р = -1/135).

4) L2(8) - L2-устойчивaя схема 8-го порядка (а = 1/54, Р = -1/216).

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

Тестирование проводилось на трех тестовых задачах, которые перечислены ниже.

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

—и(і)=

2 9

-8 -3

1 2

2 -12

и(і) на і є [0,1] с начальными условиями и(0) =

)

1

У1)

(16)

1

1

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

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

и'( і) = -( р + 2)и, + ри2

^ Иі, і є [0,2], (17)

и (і) = щ — и2 — и2

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

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

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

є =

\\«Т)-и\г )|| № )11

(18)

и(і)

и(і)

ЩУУ\ чи2Ш

\и2(/)

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

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

Ы?) Ьі(9)

А(10)

У\ Ь2(8)

X А(8)

и -и.Э -1 -1.Э -X

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

б

а

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

а

А(10) А(8)

б

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

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

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

Следующее тестирование проведено для решения с пограничным слоем (рис. 2, б). На рисунке 5 также в логарифмических осях представлены зависимости погрешностей от шага интегрирования для жесткого случая с р = 103.

Видно, что поведение погрешности сильно меняется по сравнению с предыдущими тестами. Отчетливо видны недостатки строго А-устойчивых схем, которые плохо рассчитывают переходный процесс при большом числе жесткости. Для L-устойчивых схем таких проблем не возникает, в этом и заключается их главное преимущество. При данном числе жесткости преимущество L2 над L1 почти не видно, оно проявляется при большей жесткости. Почти во всем диапазоне графика на рисунке 5 шаг по времени т > 5, 5 - размер погра-

ничного слоя. В результате погрешность крайне неравномерно распределена по отрезку интегрирования. Суммарная погрешность складывается из двух частей: £1 - погрешность на пограничном слое и £2 - погрешность на всем остальном участке интегрирования. Для методов А(8) и А(10) £ 1 очень большая. В этом и заключается главный их недостаток по сравнению с Х-устойчивыми методами. Существенное отличие А-устойчивых методов от Х-устойчивых исчезает лишь при т< 105.

lg(£)

Li(9) L2(8)

A(10) t A(8)

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

Заключение

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

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

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

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

3. Ракитский, Ю. В. Численные методы для решения жестких систем / Ю. В. Ракитский, С. М. Устинов, И. Г. Черноруцкий. - М. : Наука, 1979.

4. Dahlquist, G. Error analysis of a class of methods for stiff non-linear initial value problems / G. Dahlquist // Numerical Analysis. Lecture Notes in Mathematics 506. - Berlin : Springer, 1975. -P. 60-74.

5. Dahlquist, G. A Special Stability Problem for Linear Multistep Methods / G. Dahlquist // BIT. - 1963. - № 3. - P. 27-43.

6. Hairer, E. Solving ordinary differential equations. II. Stiff and differential-algebraic problems / E. Hairer, G. Wanner. - Berlin : Springer-Verlag, 1991.

7. Vasilev, E. High order implicit method for ODEs stiff systems / E. Vasilev, T. Vasilyeva // Korean Journal of Computational & Applied Mathematics. - 2001. - V. S, № 1. - P. 165-180.

REFERENCES

1. Vasil'ev E.I., Vasil'eva T.A., Kiseleva M.N. Mul'ti-nejavnye metody so vtoroj proizvodnoj dlja zhestkih sistem differencial'nyh uravnenij [Multi-implicit methods with the second derivative for stiff systems of differential equations]. Vеstnik Vo/GU. Sеr. 1, Matеmatika. Fizika [Journal of Volgograd State University, series 1, Mathematics. Physics], 2012, no. 2 (17), pp. 68-77.

2. Dekker K., Verver Ja. Ustojchivost' metodov Runge - Kutty d/ja zhestkih ne/inejnyh differencia/'nyh uravnenij [Stability of Runge Kutta methods for stiff nonlinear differential equations]. Moscow, Mir Publ., 19SS. 334 p.

3. Rakitskij Ju.V., Ustinov S.M., Chernoruckij I.G. Chis/ennye metody d/ja reshenija zhestkih system [Numerical methods for solving stiff systems]. Moscow, Nauka Publ., 1979.

4. Dahlquist G. Error analysis of a class of methods for stiff non-linear initial value problems. Numerica/ Ana/ysis. Lecture Notes in Mathematics 506. Berlin, Springer Verlag, 1975. Pp. 60-74.

5. Dahlquist G. A Special Stability Problem for Linear Multistep Methods. BIT, 1963, no. 3, pp.27-43.

6. Hairer E, Wanner G. So/ving ordinary differentia/ equations. ll. Stiff and differentia/a/gebraic prob/ems. Berlin, Springer-Verlag, 1991.

7. Vasilev E., Vasilyeva T. High order implicit method for ODEs stiff systems. Korean Journa/ of Computationa/ & App/iedMathematics, 2001, vol. S, no. 1, pp. 165-180.

L-STABILITY OF MULTI-IMPLICIT METHODS OF 8-th ORDER FOR DIFFERENTIAL STIFF SYSTEMS

Vasilyev Evgeniy Ivanovich

Doctor of Physics and Mathematics, Professor,

Department of Computer Science and Optimal Control,

Volgograd State University [email protected]

Prospect Universitetskij, 100, 400062 Volgograd, Russian Federation

Associate Professor,

Department of Computer Science and Optimal Control,

Volgograd State University [email protected]

Prospect Universitetskij, 100, 400062 Volgograd, Russian Federation

Postgraduate student,

Department of Computer Science and Optimal Control,

Volgograd State University [email protected]

Prospect Universitetskij, 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.

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

Vasilyeva Tatiana Anatolievna

Kiseleva Maria Nikolaevna

du(t) = f (u), dt

u(0) = u0, t > 0.

(1)

(2)

(% ) =

( 6893 , 11^ Ц3 , 9^

1814^ ^ 672

.223. + 11 й 10 + 9g

1134^ 3 у 21 ^

31 V 224

81

224

<892 - 9«

■Ц - 9Р

81

224

397 - 11^А 18144 3 ы

_Ж - 11 й

567 3 У

31

224 ,

(3)

( 1283 | 30240 а 85119а 3360 ^ уи- 269 19а 3360 ^ уи- -163 | а1 30240 ^ ы

b)= 1890 ^ й 18 | 9 'со о|чэ | 9 'со 945+й

19 -27 27 -19

V 1120 1120 1120 1120 J

At arbitrary (a, B) parameters last difference equation in (2) has 8-th order of approximating.

We found that the set of absolutely stable 3i5'D-schemes includes two one-parameter families: the set of the Z-stable schemes and the set of the schemes of heightened accuracy for linear problems. For example:

at values a = ^40, B = Xoso we have A-stable scheme with 10-th order of approximating,

at values a = 154, B = -X35 we have Z1-stable scheme with 9-th order of approximating,

at values a = 154, B = -1216 we have Z2-stable scheme with 8-th order of approximating.

The testing of this difference schemes on linear and non-linear 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 3ZSD-schemes.

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