L-устойчивый (4,2)-метод четвертого порядка для решения жестких задач Текст научной статьи по специальности «Математика»

УДК 519.622

Ь-УСТОЙЧИВЫЙ (4,2)-МЕТОД ЧЕТВЕРТОГО ПОРЯДКА ДЛЯ РЕШЕНИЯ ЖЕСТКИХ ЗАДАЧ

© 2011 Е.А. Новиков1

Исследованы (т,к)-методы решения жестких задач, в которых на каждом шаге два раза вычисляется правая часть системы дифференциальных уравнений. Показано, что максимальный порядок точности Ь-устойчивого (т,2)-метода равен четырем. Построен (4,2)-метод максимального порядка.

Ключевые слова: жесткая задача, методы Розенброка, (т,к)-методы, Ь-устой-чивость.

Введение

При решении задачи Коши для жестких систем обыкновенных дифференциальных уравнений широкое распространение получили методы типа Розенброка [1] благодаря простоте реализации и достаточно хорошим свойствам точности и устойчивости. Данные численные схемы получены из полуявных методов типа Рунге-Кутта, в которых для решения нелинейной системы алгебраических уравнений, возникающей при вычислении каждой стадии, используется одна итерация метода Ньютона [2]. Все остальные проблемы решаются выбором величины шага интегрирования. Наибольшее распространение получили методы типа Розенброка, в которых при вычислении каждой стадии применяется одна и та же матрица Якоби. Известно (см., например, [2]), что в этом случае для т-стадийного метода Розенброка максимальный порядок точности равен т +1, причем схема максимального порядка может быть только А-устойчивой. Если отказаться от максимального порядка, то можно построить Ь-устойчивую численную формулу т-го порядка точности. В практических расчетах, как правило, отказываются от максимального порядка в пользу ¿-устойчивости. Заметим, что на основе методов типа Розенброка нельзя построить схему с замораживанием матрицы Якоби выше второго порядка точности [3-4], что ограничивает применение данных методов расчетами с небольшой точностью или задачами небольшой размерности.

В [5-6] предложен класс (т, к)-методов, в которых нахождение стадий не связывается с обязательным вычислением правой части системы дифференциальных уравнений. Числа т и к означают соответственно число стадий и количество вычислений правой части на шаг интегрирования. Реализация (т, к)-методов так же проста, как и методов Розенброка, однако (т, к)-схемы имеют лучшие свойства

1 Новиков Евгений Александрович ([email protected]), Институт вычислительного моделирования СО РАН, 660036, Российская Федерация, г. Красноярск, Академгородок, ИВМ СО РАН.

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

Здесь исследуются (т, к)-методы решения жестких задач, в которых на каждом шаге два раза (к=2) вычисляется правая часть системы дифференциальных уравнений. Показано, что максимальный порядок точности ¿-устойчивого (т, 2)-метода равен четырем, и построен (4,2)-метод максимального порядка.

1. Схемы типа Розенброка

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

у' = f (y), y(to) = yo, to < t < tk, (1.1)

где y и f — вещественные N-мерные векторные функции, t — независимая переменная. Рассмотрение автономной задачи (1.1) не снижает общности. Введением дополнительной переменной y'N+1 = 1, yn +i(to) = to всегда можно неавтономную задачу

y' = f (y,t), y(to) = yo, to < t < tk, привести к автономному виду.

Для решения задачи (1.1) рассмотрим методы типа Розенброка вида

m

Уп+1 = Уп + Piki,

i=1

Dn = E — ahf, (1.2)

i-i

Dnki = hf (yn + fojkj),

i=1

где E — единичная матрица, f!n = df(yn)/dy — матрица Якоби системы (1.1), ki, 1 ^ г ^ m, — стадии метода, a, pi, eij, 1 ^ i ^ m, 1 ^ j ^ i — 1, — коэффициенты, определяющие свойства точности и устойчивости (1.2). В настоящее время методы типа Розенброка трактуются более широко. Под ними понимаются все численные схемы, в которых матрица Якоби или ее аппроксимация вводятся непосредственно в формулу интегрирования. Численные формулы вида (1.2) можно получить из класса полуявных методов типа Рунге-Кутта, если в них при вычислении каждой стадии ограничиться одной итерацией метода Ньютона. Однако в (1.2) при вычислении стадий необходимо решать только линейные системы алгебраических уравнений, в то время как в неявных или полуявных методах типа Рунге-Кутта требуется использовать итерационный процесс типа ньютоновского.

Рассмотрим в качестве примера одностадийный метод типа Розенброка

Уп+1 = Уп + piki, Dnki = hf (yn), (1.3)

где матрица Dn определена в (1.2). Разлагая приближенное решение yn+i в ряд Тейлора по степеням h до членов с h2 включительно, получим

Уп+1 = Уп + Pihfn + apihfn fn + O(h3),

где fn = f(yn) и fn = df(yn)/dy — матрица Якоби системы (1.1). Представление точного решения y(tn+i) в виде ряда Тейлора в окрестности точки tn имеет вид

y(tn+i) = y (tn) + hf + 0, 5h2f' f + O(h3),

где элементарные дифференциалы f и f'f вычислены на точном решении y(tn). Сравнивая ряды для точного и приближенного решений видим, что одностадийная схема Розенброка (1.3) будет иметь второй порядок точности, если pi = 1 и api = 0, 5, то есть pi = 1 и a = 0, 5.

Теперь исследуем устойчивость данной численной формулы. Для этого применим его для решения скалярного тестового уравнения y' = Ay, где Л есть произвольное комплексное число, Ж(Л) < 0. Параметр Л интерпретируется как некоторое собственное число матрицы Якоби задачи (1.1). Обозначая x = hA, получим yn+i = Q(x)yn, где функция устойчивости Q(x) имеет вид

Q(x) = [1 + (pi — a)x]/(1 — ax).

Подставляя сюда значения коэффициентов pi = 1 и a = 0, 5, имеем

Q(x) = (1+0, 5x)/(1 — 0, 5x),

то есть схема (1.3) второго порядка точности является A-устойчивой. Из вида функции устойчивости следует, что эта схема будет L-устойчивой, если pi = a =1, что противоречит второму порядку точности. Обычно отказываются от второго порядка в пользу L-устойчивости, что приводит к более эффективному методу, хотя и первого порядка.

В случае большой размерности задачи (1.1) основные вычислительные затраты связаны с обращением матрицы Dn. Обычно вместо обращения решается линейная система алгебраических уравнений Dnki = hf(yn) с применением LU-разложения матрицы Dn с выбором главного элемента либо по столбцу, либо по строке, а иногда по всей матрице. Декомпозиция матрицы Dn приводит к порядку N 3 арифметических операций. Обратный ход метода Гаусса стоит порядка N2 арифметических операций. Таким образом, при большой размерности исходной задачи (1.1) общие вычислительные затраты определяются временем декомпозиции матрицы Dn.

Одновременно с численной формулой Розенброка (1.3) рассмотрим метод следующего вида:

Уп+1 = Уп + piki + p2&2, D„ki = hf(y„), D„k2 = ki. (1.4)

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

pi + p2 = 1, a(pi + 2p2) = 0, 5.

Условие L-устойчивости (1.4) записывается в виде квадратичного уравнения относительно коэффициента a, то есть

a2 + 2a + 0, 5 = 0.

Это уравнение имеет два корня ai = 1—0, 5л/2 и ai = 1+0, 5л/2. Обычно в расчетах применяется a = ai , потому что в этом случае меньше коэффициент в локальной ошибке схемы (1.4). В результате имеем коэффициенты

a =1 — 0, 5%/2, pi = a, p2 = 1 — a

L-устойчивой численной формулы (1.4) второго порядка точности.

Отметим, что в случае большой размерности задачи (1.1) методы (1.3) и (1.4) различаются на количество арифметических операций, необходимых для выполнения обратного хода в методе Гаусса, которые на фоне декомпозиции матрицы Dn

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

Применение стадий типа было положено в основу (т, к)-методов.

2. Класс (ш,к)-методов решения жестких задач

Пусть Z есть множество целых чисел, и заданы m, k G Z, k ^ m. Обозначим через Mm, Mk и J¿ множества вида

Mm = {i G Z | 1 < i < m}, Mk = {m¡ G Mm|1 = mi < m2 < ... < mj ^ m}, (2.1)

J¿ = {mj_i G Mm|j > 1, mj G Mj, mj ^ i}, 1 < i ^ m. Рассмотрим следующие численные схемы:

m

Уп+i = Уп + P¿k¿, D„ = E - ah/П,

i=i

i-1 i-1 D„k¿ = h/(yn + ^ e¿j kj aij kj + h/' ^2 kj, i G Mj, (2.2)

j=i jeJ j=i

i-i

D„k¿ = kj-i + ^ a¿jkj + h/' ^ cjkj, i G Mm\Mj,

jeJ j=i

где k¿ стадии метода, a, , e¿j, ®¿j и Cj постоянные коэффициенты, h шаг интегрирования, E — единичная матрица, /' = д/(yn)/ду — матрица Якоби системы (1.1), k — количество вычислений функции /, m — число стадий или количество обратных ходов в методе Гаусса. На каждом шаге интегрирования осуществляются одно вычисление матрицы Якоби и одна декомпозиция матрицы Dn. Допускается аппроксимация матрицы Якоби / матрицей An вида

An = /' + h2B„ + O(h3),

где матрица Bn не зависит от размера шага интегрирования. Данное условие позволяет применять (2.2) с замораживанием как аналитической, так и численной матрицы Якоби [7]. Так как k и m полностью определяют затраты на шаг, а набор чисел mi, •••, mj из множества Mj только распределяет их внутри шага, то методы типа (2.2) названы (m, ^-методами.

Отметим, что при k = m и a¿j = Cj = 0 методы (2.2) совпадают со схемами типа Розенброка (1.2), а при k = m и a¿j =0 — с ROW-методами [2]. В отличие от ROW-методов в численных формулах (2.2) более точно определены затраты на шаг интегрирования и более правильно описана область определения коэффициентов численных формул, что упрощает их исследование и делает их предпочтительнее.

Заметим так же, что при рассмотрении методов такого типа, как правило, все авторы изучали случай k = m, то есть когда число стадий и количество вычислений правой части совпадают. В этом случае k-стадийную схему (2.2) можно поставить в соответствие k-стадийной полуявной формуле типа Рунге-Кутта [2], при реализации которой на каждом шаге используется одна матрица размерности N. Относительно таких численных формул известно, что нельзя построить k-стадийную схему выше (k + 1)-го порядка точности.

При рассмотрении методов (2.2) при т > к можно показать, что при к, равном 1 и 2, можно построить численные формулы порядка точности, равном 2к. Таким образом, среди методов (2.2) имеются такие численные схемы, что в смысле точности они не хуже неявных методов типа Рунге-Кутта и в то же время они требуют существенно меньших вычислительных затрат при реализации.

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

3. Представление стадий в виде рядов Тейлора

При исследовании (т, к)-методов (2.2) потребуются представления в виде рядов Тейлора в окрестности точки у„ до членов с Л4 включительно стадий , 1 ^ г ^ 4, вида

£„к1 = Л/(у„), Д„к2 = к1, £„кз = Л/(у„ + вз1к1 + вз2к2) + «З2к2, (3.1)

П„к4 = кз + «42к2.

Данные ряды можно записать следующим образом:

к1 = / + а/ / + а2 // + а3// + 0(Л5), к2 = / + 2аЛ2/ / + 3а2/2/„ + 4а3//„ + 0(Л5), кз = (1 + «21)Л/„ + (а + ^31 + в32 + 3а«32)Л2/,П /„+ +Л3[(а2 + 2ав31 + 3а^32 + 6а2 «32)// + 0, 5(^31+

+в32)2/П'/П] + Л4[(а3 + 3а2в31 + 6а2^32 + +10а3а32)/П3/п + а(в31 + ^32)(^31 + 2^32)// /+ (3.2)

+0, 5а(в31 + в32)2/,;// + 1(в31 + в32)3//] + 0(Л5), к4 = (1 + «32 + «42)Л/„ + (2а + в31 + в32 + 4а«32 + 3аа42)Л2/„/„+

+Л3[(3а2 + 3ав31 + 4а^32 + 10а2 «32 + 6а2 «42)//+ +0, 5(в31 + в32)2//] + Л4[(4а3 + 6а2в31 + 10а3а2 + +10а3а42 )// + а(в31 + ^32)(^31 + 2^32)///2+

+а(в31 + в32)2/' /^'/2 + 6(в31 + в32)3/Г /3] + 0(Л5).

При записи (3.2) применялось представление П—1 в виде ряда Тейлора

П-1 = Е + а/ + а2/ + а3/ + • • • = £

¿=0

Ниже соотношения (3.2) будут использоваться при построении конкретных методов. Далее потребуется разложением точного решения у(£„+1) в ряд Тейлора в окрестности точки до членов с Л4 включительно, которое имеет вид

у(*п +1) = у(*п) + Л/ + 0, 5Л2/'/+

+1 ^3[/'2/ + / /21 + 1 ^4[/'3/ + (3.3)

6 24

+/7 "/2 + 3/ "/72 + / "7 31 + о(^5),

где элементарные дифференциалы вычислены на точном решении у(£и).

4. Ь-устойчивый метод четвертого порядка

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

Теорема. Пусть в формуле (2.2) имеет место к = 2. Тогда при любом выборе множеств (2.1) и при любом значении параметра т нельзя построить метод (2.2) выше четвертого порядка точности.

Для простоты доказательство проведем для скалярной задачи (1.1), точное решение которой можно записать в виде

Л2 1

У(*п+1) = У(*п) + / + у /'/ + 6 ^3[/'2/ + /"/2] +

+24 ^4[/'3 / + 4/'/"/2 + / '"/3]+ (4.1)

+ Т2о ^ [/'4/ + 4/'/'"/3 + 5/'2/"/2 + / "2/3 + /^ /41 + °(^6),

где элементарные дифференциалы вычислены в точке у(£и).

С использованием (3.2) имеем, что при любом выборе множества М^ второе вычисление функции / будет осуществляться в точке

Уп,с = Уп + ^ + O(h5),

j=1

где cj, 1 ^ i ^ 4, определяются через коэффициенты схемы (2.2). Поэтому для доказательства сформулированного утверждения достаточно показать, что в представлении функции h/(yn,c) в виде ряда Тейлора в окрестности точки yn не содержится слагаемого h5/ПП/2/П. В этом случае в разложении (4.1) имеется слагаемое h5///2/3, а в соответствующих представлениях kj, 1 ^ i ^ m оно отсутствует. Разлагая h/(yn,c) в ряд Тейлора, будем иметь

h/(Уп,с) = h/n + cih2///п + h3// + 0, 5c2//] + +h4[c3/f/„ + ci/ // + 1 /731 +

+h5[c4/,/4/„ + C1C3/,?//+

+0, 5c1c2////«/з + 24//4] + O(h6),

что завершает доказательство.

Перейдем к исследованию (m, к)-методов (2.2) при k = 2. Выберем Mm = {i|1 ^ ^ i ^ m}, M2 = {1, 3}, тогда Jj = {2}, 3 ^ i ^ m. Рассмотрим численные схемы вида

m

Уп+1 = Уп + Pikj, j=1

Dkl = h/(yn), Dk2 = ki,

D„k3 = hf (y„ + esiki + ^32^2) + «32^2, (4.2)

D„k¿ = kj_i + «¿2^2, 4 < i < m,

где матрица Dn определяется по второй формуле (2.2).

Пусть m = 4. Подставим разложения k¿, 1 ^ i ^ 4, из (3.2) в первую формулу (4.2). Полагая yn = y(tn) и сравнивая (3.3) и полученное представление приближенного решения до членов с h4 включительно, получим условия четвертого порядка точности схемы (4.2), которые имеют вид

Р1 + Р2 + (1 + «32^3 + (1 + «32 + «42)Р4 = 1,

api + 2ap2 + (a + вз1 + вз2 + За«з2)рз+ +(2a + вз1 + вз2 + 4ааз2 + 3a«42)p4 = 1/2, a2pi + 3a2p2 + (a2 + 2a^3i + Завз2 + 6a2 «42^3+ +(3a2 + 3a^3i + 4авз2 + 10a2 «32 + 6a2 «42^4 = 1/6, a3pi + 4a3p2 + (a3 + 3a2^3i + 6a2 ^32 + 10a3 «32^3 + (4a3 + 6a2^3i + +10а2в32 + 2a3 «32 + 10a3«42)p4 = 1/24,

(Ají + ^32)2(p3 + p4) = 1/3, a(^3i + вз2)(вз1 + 2^32)(p3 + p4) = 1/8, aOSsi + вз2)2 (0, 5pi + p4) = 1/24, (вз1 + ^32)3(p3 + p4) = 1/4.

Данная система является нелинейной относительно коэффициентов p¿ и a¿j, что вызывает определенные трудности с ее исследованием. Поэтому при изучении методов (4.2) поступим следующим образом. Рассмотрим вспомогательную численную схему

6

Уп+1 = Уп + Е c¿d¿, ¿=1

Dndi = hf(yn), D„¿2 = di, Dn dj = hf (y„ + e3idi + вз2 d2), (4.3)

D„d4 = d2, Dnd5 = d3, Dnd6 = d4.

Подставляя d¿, 1 ^ i ^ 6, из (4.3) в формулу (4.2), получим

Уп+1 = Уп + pidi + p2 d2 + p3d3 + («32p3 + «42p4)d4 + p4d5 + «32p4d6.

Сравнивая полученное соотношение с (4.3), видим, что параметры схем (4.2) и (4.3) связаны соотношениями

pi = Ci, p2 = С2, p3 = Сз, p4 = C5,

«32 = Сб/С5, «42 = (C4C5 - СзСб)/С2 (4.4)

Теперь перейдем к изучению численной формулы (4.3). Для этого требуются разложения стадий d¿, 1 ^ i ^ 6, в ряды Тейлора в окрестности точки yn до членов с

h4

включительно, которые легко получить из (3.2). Подставим их в первую формулу (4.3). Полагая yn = y(tn) и сравнивая полученное представление приближенного решения и (3.3), получим условия четвертого порядка точности схемы (4.3), то есть

6

1) Е c¿ = 1, ¿=i

2) ас1 + 2ас2 + (а + ^31 + в32)с3 + 3ас4+ +(2а + в31 + в32)с5 + 4асб = 1/2, 3) а2С1 + 3а2С2 + (а2 + 2а^31 + 3а^32)с3 + 6а2С4+ + (3а2 + 3ав31 + 4ав32)с5 + 10а2Сб = 1/6,

4) (Ал + в32)2(с3 + с5) = 1/3. (4.5)

5) а(в31 + в32)(в31 + 2в32)(с3 + с5) = 1/8, 6) а(в31 + в32)2(0, 5с3 + с5) = 1/24, 7) (в31 + в32)3(с3 + с5) = 1/4,

8) а3с1 + 4а3с2 + (а3 + 3а^31 + 6а^32)с3 + 10а3с4+ +(4а3 + 6а2в31 + 10а2в32)с5 + 20а3сб = 1/24. Исследуем совместность (4.5). Из четвертого и седьмого уравнений имеем

в31 + в32 = 3/4, с3 + с5 = 16/27. (4.6)

Тогда из пятого и шестого соотношений (4.5) получим ^32 и с3 соответственно. Из (4.6) выразим ^31 и с5. Учитывая первое уравнение (4.5), подставим полученные выражения для ^31, в32, с3 и с5 во второе, третье и восьмое равенства (4.5). Получим линейную систему относительно параметров с2, с4 и с6, которая имеет вид

ас2 + 2ас4 + 3асб = -(22а + 5)/54,

2а2с2 + 7а2с4 + 16а2с6 = (64а2 + 29а - 30)/54,

а2с2 + 3а2 с4 + 6а2 сб = (32а2 — 11а — 6)/54.

Разрешив данную систему, из первого уравнения (4.5) получим с1. В конечном результате будем иметь

с1 = (76 — 29/а + 3/а2)/27, с2 = ( —146 + 89/а — 12/а2)/27,

с3 = (32 — 4/а)/27, с4 = (72 — 59/а +10/а2)/18,

с5 = (4/а — 16)/27, с6 = ( — 18 + 19/а — 4/а2)/18, (4.7)

в31 = (48 — 9/а)/32, ^32 = (9/а — 24)/32.

Применяя (4.3) для решения задачи у' = Ау, Ж(А) < 0, получим условие ¿-устойчивости

а(с1 + с3) — а2 — с3^31 = 0.

Функция устойчивости не приводится в силу ее громоздкости. Подставляя в последнее равенство значение параметров (4.7), будем иметь

24а4 — 96а3 + 72а2 — 16а + 1 = 0. (4.8)

Подставляя теперь (4.7) в (4.4), получим параметры схемы (4.2), то есть

р1 = (76 — 29/а + 3/а2)/27,

р2 = ( — 146 + 89/а — 12/а2)/27,

р3 = (32 — 4/а)/27, р4 = (4/а — 16)/27,

в31 = (48 — 9/а)/32, ^32 = (9/а — 24)/32, (4.9)

а32 = [-54а + 57 - 12/а]/(8 - 32а),

«42 = [-864а2 + 828а - 288 + 36/а]/(4 - 16а)2,

при которых она имеет четвертый порядок точности. Параметр а определяется из условия L-устойчивости (4.8). Заметим, что уравнение (4.8) имеет два вещественных корня

а! = 0, 5728160624821, а2 = 3,100316735116.

Из результатов расчетов ряда задач [2] следует, что наиболее подходящим является корень а = а!, который приводит к более эффективному и, что существеннее, более надежному алгоритму интегрирования.

Заключение

При k = 2 в классе (m, к)-методов можно построить численную схему, которая по свойствам точности не уступает неявному методу типа Рунге-Кутта с двумя вычислениями правой части дифференциальной задачи, а при линейном анализе устойчивости она также не хуже. В то же время при записи (4.2) сразу заложен способ ее реализации, то есть до начала расчетов можно оценить вычислительные затраты на шаг интегрирования. Что касается неявных методов типа Рунге-Кут-та, то для них вычислительные затраты в сильной степени зависят от способа реализации и, в частности, использование двухстадийной схемы вовсе не означает, что на каждом шаге будет два раза вычисляться правая часть дифференциальной задачи. Поэтому на некоторых задачах (m, к)-методы будут выгоднее неявных формул типа Рунге-Кутта.

Здесь не обсуждается вопрос о построении неравенства для контроля точности вычислений и для выбора величины шага интегрирования в методе (4.2). Это неравенство легко можно получить с применением вложенных методов, введя дополнительную стадию Dk = k4 + «52 k2.

Литература

[1] Rosenbrock H.H. Some general implicit processes for the numerical solution of differential equations // Computer. 1963. № 5. P. 329-330.

[2] Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М.: Мир, 1999. 685 c.

[3] Novikov V.A., Novikov E.A., Umatova L.A. Freezing of the Jacobi matrix in the Rosenbrock type method of the second order accuracy // Proc. BAIL-IV Conf.: Bool Press, 1986. P. 380-386.

[4] Новиков Е.А., Двинский А.Л. Замораживание матрицы Якоби в методах типа Розенброка // Вычислительные технологии. 2005. Т. 10. С. 108-114.

[5] Новиков Е.А. Об одном классе одношаговых безытерационных методов решения жестких систем // Актуальные проблемы вычислительной и прикладной математики. 1987. С. 138-139.

[6] Новиков Е.А., Шитов Ю.А., Шокин Ю.И. Одношаговые безытерационные методы решения жестких систем // ДАН СССР. 1988. Т. 301. № 6. С. 1310-1314.

[7] Новиков Е.А., Шитов Ю.А. Алгоритм интегрирования жестких систем на основе (т,к)-метода второго порядка точности с численным вычислением матрицы Якоби // Препринт № 20. Красноярск: ВЦ СО АН СССР, 1988. 23 с.

Поступила в редакцию 18/V/2011;

в окончательном варианте — 19/VT/2011.

L-STABLE (4,2)-METHOD OF THE FOURTH ORDER FOR SOLVING STIFF PROBLEMS

© 2011 E.A. Novikov2

(M,k)-methods for solving stiff problems, in which on each step two times the right-hand side of the system of differential equations is calculated are investigated. It is shown that the maximum order of accuracy of the L-stable (m,2)-method is equal to four. (4,2)-method of maximal order is built.

Key words: stiff problem, methods of Rosenbrock, (m,k)-methods, L-stability.

Paper received 18/V/2011. Paper accepted 19/W/2011.

2Novikov Evgeniy Alexandrovich (novikovaicm.krasn.ru), the Institute of Computational Modelling, SB RAS, Krasnoyarsk, 660036, Russian Federation.