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

Вычислительные технологии Том 12, № 5, 2007

ИССЛЕДОВАНИЕ (m, 2)-МЕТОДОВ РЕШЕНИЯ ЖЕСТКИХ СИСТЕМ*

E. A. Новиков Институт вычислительного моделирования СО РАН,

Красноярск, Россия e-mail: [email protected]

The (m, 2)-methods are investigated. It has been shown that the maximum order of accuracy for the L-stable (m, 2)-method is equal to four. The coefficients for the method of maximum accuracy are obtained.

Введение

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

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 02-01-00523) и Президентской программы "Ведущие научные школы РФ" (грант № НШ-3428.2006.9).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2007.

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

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

1. Методы типа Розенброка

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

у' = /(у), у(го) = Уo, го < г < 4, (1)

где у и / — вещественные N-мерные векторные функции; г — независимая переменная, которая изменяется на заданном конечном интервале. Известно, что неавтономную систему введением дополнительной переменной можно привести к автономному виду. Поэтому рассмотрение (1) не снижает общности. Ниже потребуется представление точного решения у(гп+1) задачи (1) в виде ряда Тейлора в окрестности точки гп, которое имеет вид

У(*п+1) = у(*п) + к/ + ^ /'/ + у[/'2/ + /''/2] +

+ // + /'/''/2 + З/''/'/2 + /'''/3] + 0(к5), (2)

где элементарные дифференциалы вычислены на точном решении у(гп). Методы типа Розенброка применительно к задаче (1) имеют вид

т

уп+1 = уп + Е ргкг

г=1

г-1 \ /г-1

^п,.кг = к/ ( уп + Е вгзк , Бп,г = Е — а.Н/п уп + Е Ъзкз ) , (3)

3=1 ) \ 3=1

где к — шаг интегрирования; Е — единичная матрица; /п = д/(уп)/ду — матрица Якоби векторной функции /(у); аг, рг, вз и 7.3, 1 < г < т, 1 < ] < г — 1, — числовые коэффициенты, определяющие свойства точности и устойчивости (3). В настоящее время методы типа Розенброка трактуются более широко [2].

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

Наиболее эффективные реализации методов (3) возникают при а1 = а2 = • • • = ат = а и = 0) 1 < г < ш, 1 < ] < г — 1. Соответствующие численные схемы имеют вид

т / г— 1 \

Уп+1 = Уп + ЕРгкг, Епк = к/ Уп + ^ вг]к , Еп = Е — а/(Уп). (4) г=1 \ ]=1 /

В этом случае на каждом временном шаге требуется обращение одной матрицы Еп размерности N. Вместо обращения матрицы обычно решается линейная система алгебраических уравнений с применением Ш-разложения матрицы Еп. Декомпозиция Еп на верхнюю и нижнюю треугольные матрицы приводит к порядка N3 арифметических операций, в то время как последующее вычисление стадий стоит порядка N2 операций. Поэтому в случае большой размерности задачи (1) время декомпозиции матрицы Еп практически полностью определяет общие вычислительные затраты. Рассмотрим в качестве примера одностадийный метод типа Розенброка

Уп+1 = Уп + Р1&1, Опк1 = к/(уп), Еп = Е — аН/'п (уп). (5)

Требование второго порядка точности приводит к соотношениям р1 = 1 и ар1 = 0.5, в то время как условие ¿-устойчивости означает а =1 и оно противоречит второму порядку. Поэтому в настоящее время наиболее известен набор коэффициентов р1 = а = 1 ¿-устойчивого метода (5) первого порядка точности.

2. Класс (т, к)-методов

Прежде чем формулировать класс (ш, к)-методов в общем виде, рассмотрим следующую численную формулу:

Уп+1 = Уп + Р1к1 + Р2к2, Епк1 = к/(уп), Епк2 = кь (6)

В случае большой размерности задачи (1) после введения дополнительной стадии к2 не происходит существенного увеличения вычислительных затрат на фоне декомпозиции матрицы Еп, т.е. вычислительные затраты для численных формул (5) и (6) различаются несущественно. В то же время при р1 = а, р2 = 1 — а, где а есть корень уравнения

а2 — 2а + 0.5 = 0, (7)

метод (6) имеет второй порядок точности и является ¿-устойчивым. Это было использовано при описании класса (ш, к)-методов. Уравнение (7) имеет два корня а1 = 1 — 0.5л/2 и а2 = 1 + 0.5у/2. Обычно для практических расчетов применяется корень а = 1 — 0.5л/2) потому что в этом случае меньше коэффициент в локальной ошибке. Заметим, что в отличие от (5) контролировать точность вычислений и выбор величины шага интегрирования численной схемы (6) можно проверкой неравенства ||к2 — к1|| < е, где || • || — некоторая норма в Ки, е - требуемая точность интегрирования.

Класс (ш, к)-методов [4, 5] вводится следующим образом. Пусть заданы целые положительные числа ш и к, к < ш. Обозначим через Мт множество целых чисел г, 1 < г < ш, а через Мк и . — подмножества из Мт вида

Мк = {шг € Мт|1 = ш1 < ш2 < • • • < шк < ш}, . = {ш]—1 € Мт|^' > 1,ш] € Мк,ш] < г}, 2 < г < ш. (8)

Тогда (т, к)-методы можно представить в виде

уп+1 — Уп + Е ргкг

г=1

Е - ак/П-.

Епкг — к/ [уг

г-1

+ Е А? к3 3 = 1

+ Е аг3 к3 , 3^1

г е мк,

(9)

Е,пкг — кг-1 + Е ®гзк, г е мт \ Мк.

3^1

Множество .]г, 2 < г < т, служит для устранения "лишних" коэффициентов агз, за счет которых нельзя повлиять на свойства точности и устойчивости (9) и которые линейно выражаются через другие коэффициенты. Отметим, что рассмотрение Е3=1 азкз вместо Е3<^- аг3к] приводит к существенным трудностям при построении конкретных численных формул именно по причине наличия "лишних" коэффициентов. Заметим также, что в традиционных одношаговых методах для описания вычислительных затрат на шаг интегрирования достаточно одной константы т — числа стадий, ибо в данных методах каждая стадия сопровождается обязательным вычислением правой части задачи (1). В методах (9) есть два вида стадий — для некоторых требуется вычисление правой части, а для других не требуется. В результате в (9) для описания вычислительных затрат на шаг требуются две постоянные т и к. Затраты на шаг следующие: один раз вычисляется матрица Якоби и осуществляется декомпозиция матрицы Еп, к раз вычисляется функция / и т раз выполняется обратный ход в методе Гаусса. В случае к — т и аг3 — 0 численные схемы (9) совпадают с методами типа Розенброка. В остальных случаях это другие методы, обладающие лучшими свойствами по сравнению с (3).

т

3. Максимальный порядок точности (т, 2)-методов

Рассмотрим (т, 2)-методы следующего вида:

т

Уп+1 — Уп + ЕРгкг, Е,п — Е - ак/п,

г=1

Епк1 — к/(уп), Епкг — кг-1, 2 < г < ^ - 1,

( *!-1 \ Епк.31 — к/ I уп + Рз1,33 + 0^1-1^-1, (10)

Епкг — кг—1 + <^1-1^1-1, + 1 < г < т,

где 51 и т, 51 < т, — произвольные целые постоянные. Нетрудно видеть, что (10) описывают всевозможные варианты (т, 2)-методов.

Теорема. Пусть в (9) имеет место к — 2. Тогда при любом выборе множеств (8) и при любом числе стадий т нельзя построить (т, 2)-метод выше четвертого порядка точности.

Без потери общности для простоты доказательство проведем для скалярной задачи (1), точное решение у(£п+^ которой можно записать в виде

у(<*+1) — У(<п) + к/ + к2 /'/ + у [/'2/ + //2] + к4 [/'3 / + 4/'/" /2 + /'''/ 3] +

+ т20[/'4/ + 4/'/'''/3 + 5/'2/''/2 + /''2/3 + /1¥ /4] + °(к6), (11)

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

Рассмотрим (m, 2)-методы (10). Учитывая

D-1 = E + ah/ + a2h2/n + a3h3/f + a4h4/n + O(h5), (12)

имеем, что второе вычисление функции /(y) будет осуществляться в точке

si-1 4

Уп,с = Уп + ^ kj = Уп + ^ с^/(г-1)/ + O(h5), j=1 i=1

где ci, 1 < i < 4, определяются через коэффициенты схемы (10), а элементарные дифференциалы вычислены на приближенном решении yn. Теперь, учитывая (11) и (12), для доказательства теоремы достаточно того, что в представлении функции /(yn>c) в виде ряда Тейлора по степеням h не содержится слагаемого h5/'2/3 Разлагая /(yn>c) в ряд Тейлора в окрестности точки yn до членов с h5 включительно, имеем

c3

/(yn,c) = h/n + C1h2/n/n + h3 [c2/n2/n + 0.5/П'/П] + h4// + wn//2 + ^/П'/П] +

c4

+ h5// + С1С3/П/П / + 0.5c? С2/П/П7П + 24 /n1V /4] + O(h6),

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

4. A-устойчивый (m, 2)-метод четвертого порядка

Выберем сначала множества (8) следующим образом:

Mm = {1, 2,..., m}, Mk = {1, 2}, J = {1}, 2 < i < m, т. е. рассмотрим численные схемы вида

m

yn+1 = yn + ^Piki, Dn = E - ah/, i=1

D„k1 = h/(y„), (13)

D„k2 = h/(yn + ^21k1) + а21&1, Dnki = ki-1 + ai1 k1, 3 < i < m.

Случай 1. Пусть m =2. Подставим разложения в виде рядов Тейлора для k1 и k2 в первую формулу (13). Полагая yn = y(tn) и сравнивая ряды для точного (2) и приближенного yn+1 решений до членов с

h3

включительно, получим условия третьего

порядка точности схемы (13), т.е.

1) Р1 + (1 + «21)Р2 = 1;

2) ap1 + (a + ^21 + 2aa21)p2 = 0.5; (14)

3) a2p1 + (a2 + 2a^21 + 3a2«21)^2 = 1/6;

4) 3в21Р2 = 1.

Исследуем совместность (14). Умножим первое уравнение (14) на 2a, а второе на —1 и сложим. Затем умножим первое соотношение на 3a2, а третье на —1 и сложим. Будем иметь

ap1 + ap2 — в21р2 = 2a — 0.5, 2a2p1 + 2a2p2 — 2a^21p2 = 3a2 — 1/6.

Теперь, умножая первое из полученных равенств на -2а и складывая со вторым, получим

а2 - 6а + 1 — 0. (15)

Пусть в21 свободный коэффициент. Тогда из последнего уравнения (14) выразим р2, из первого и второго соотношений, которые теперь линейны относительно р1 и а21, вычислим оставшиеся коэффициенты. В результате запишем

1

Р1 — 0.5[в21(12ав21 - 3^21 + 2) - 2а], р2 - 2 ,

«21 — 0.5^21(3^21 - 6ав21 - 2)/а, (16)

где а определяется из уравнения (15).

Исследуем устойчивость схемы (13) на линейном скалярном уравнении:

у' — Ау, у(0) — У0, г > 0. (17)

При т — 2, применяя (13) для решения (17), получим

^ ч ч 1 + (1 - 2а)х + (а2 - 2а + 0.5)х2

Уп+1 — д(х)уп, ^(х) —-^-^ 1 ,2-—, х — кА. (18)

(1 - ах)2

При получении (18) использовались соотношения (16). Из (18) следует, что условие ¿-устойчивости схемы (13) имеет вид выражения (7), которое противоречит (15), т.е. при т — 2 метод (13) третьего порядка точности ¿-устойчивым не является. Уравнение (15) имеет два вещественных корня: а1 — (6+\^12)/12 и а1 — (6--\/12)/12. Нетрудно убедиться, что при а1 —

(6 + 1/12)/12 численная схема (13) с коэффициентами (16) является

А-устойчивой.

Свободный коэффициент в21 можно использовать для минимизации главного члена локальной ошибки. При условии 4в21р2 — 1 слагаемое с элементарным дифференциалом к4/'''/3 в локальной ошибке отсутствует. Используя это соотношение, окончательно получим коэффициенты А-устойчивого метода (13) третьего порядка точности с минимальной локальной ошибкой, которые имеют вид

6 + У12 76а - 3 16 „ 3 3 - 54а

а — Р1 — Р2 — 27, в21 — 4, а21 — ^2^

Случай 2. Исследуем схему (13) при т — 3. Подставим разложения в виде рядов Тейлора для к1, к2 и к3 в первую формулу (13). Полагая уп — у(гп) и сравнивая ряды для точного (2) и приближенного уп+1 решений до членов с

к3

включительно, получим

условия третьего порядка точности схемы (13), т.е.

1) Р1 + (1 + «21)Р2 + (1 + «21 + «31)Р3 — 1,

2) аР1 + (а + ^21 + 2аа21)Р2 + (2а + ^21 + 3аа21 + 2аа31)р3 — 0.5, (19)

3) а2Р1 + (а2 + 2ав21 + 3а2а21)р2 + (3а2 + 3ав21 + 6а2а21 + 3а2а31 )р3 — 1/6,

4) 3в221 (Р2 + Р3) — 1.

Применяя (13) для решения (17), получим условие ¿-устойчивости

а(р1 + Р2) - а2 - в21Р2 — 0. (20)

Выражение для функции устойчивости не приводится в силу ее громоздкости.

Исследуем совместность (19) и (20). Умножим первое уравнение (19) на 6а2, а второе на -3а и сложим. Затем умножим первое уравнение (19) на 3а2, а третье на —1 и сложим. В результате получим

3a2pi + 3а2р2 — 3aß2ip2 — 3а2а21р3 — 3а,021р3 = 6а2 — 1.5а, 2a2pi + 2а2р2 — 2а,в21р2 — 3а2а21р3 — 3а,021р3 = 3а2 — 1/6.

Умножим второе из полученных равенств на —1 и сложим с первым. Учитывая (20), получим условие L-устойчивости схемы (13), т.е.

6а3 — 18а2 + 9а — 1 = 0. (21)

Отметим, что (21) есть известное уравнение L-устойчивости трехстадийных схем, т. е. при одинаковом числе стадий функции устойчивости методов типа Розенброка и (m, k)-схем совпадают.

Пусть а21 и ß21 являются свободными коэффициентами. Умножим первое уравнение

(19) на —2а и сложим со вторым, учитывая (20). Из полученного соотношения выразим p3. Затем последовательно определим из четвертого равенства (19) коэффициент p2, из

(20) — pi, из первого уравнения (19) — а31. В результате будем иметь

аСз + С4 — C5 2c2 — 3ß2i Ci Ci

Pi = -, P2 = -T^ä-, P3 = — ,

C3 6ß2i C2 2C2

2C2 [(1 — а)сз — C4 + C5 — 2ас2 (1 + a2i)]

«3i =-,

C2C3

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

ci = 1 — 4а + 2а2, c2 = ß21 + аа21, c3 = 6ас2 ß|i,

C4 = 2c2(ß2i — а), C5 = 3ciß2i(ß2i — а).

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

Случай 3. Пусть m = 4, т. е. рассмотрим численную схему вида

4

Уп+i = Уп + ^ Piki, D = E — а^/П,

i=1

Dnki = h/(yn), D„k2 = h/(y„ + ß2iki) + a2iki, (22)

D„k3 = k2 + a3iki, D„k4 = кз + a4ikb

Подставим разложения ki, 1 < i < 4, в виде рядов Тейлора по степеням h в первую формулу (22). Полагая yn = y(tn) и сравнивая точное решение (2) и полученное пред-

ставление приближенного решения уп+1 до членов с

к4

включительно, получим условия

четвертого порядка точности, т. е.

1) Р1 + (1 + «21)Р2 + (1 + «21 + «31)Р3 + (1 + «21 + «31 + «41)Р4 — 1,

2) ар1 + (а + ^21 + 2аа21)р2 + (2а + Д21 + 3а«21 + 2а«31 )р3 +

+ (3а + в21 + 4а«21 + 3а«31 + 2а«41)р4 — 0.5,

3) а2р1 + (а2 + 2ав21 + 3а2«21 )р2 + (3а2 + 3ав21 + 6а2«21 + 3а2«31)р3 +

Г) Л Л Оч /

+ (6а + 4ар21 + 10а «21 + 6а «31 + 3а «41)р4 — 1/6,

4) в21(Р2 + Р3 + Р4) —1/3, (23)

5) ав21 (Р2 + Р3 + Р4) — 1/8,

6) ав21(0.5р2 + Р3 + 1.5р4) — 1/24,

7) в31(Р2 + Р3 + Р4) — 1/4,

8) а3р1 + а2 (а + 3^21 + 4а«21)р2 + а2 (4а + 6^21 + 10а«21 + 4а«31)р3 +

+а2(10а + 10в21 + 2а«21 + 10а«31 + 4а«41 )р4 — 1/24.

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

7

Уп+1 — Уп + ^ Ыг, Бп — Е - ак/п, г=1

Бп^1 — к/(Уп), Бп^ — к/(Уп + в21кх), (24)

Бп^г — ¿г-1, 2 < г < 4, 6 < г < 7.

Подставляя йг, 1 < г < 7, в (22), получим

Уп+1 — Уп + Р1^1 + («21Р2 + «31Р3 + «41Р4)^2 + («21Р3 + «31Р4)4 + + «21Р4^4 + Р2 4 + Р3^6 + Р4^7.

Отсюда видим, что коэффициенты методов (22) и (24) связаны соотношениями

Р1 — С1, Р2 — С5, Р3 — Сб, Р4 — С7, «21Р4 — С4, «21Р3 + «31Р4 — С3, «21Р2 + «31Р3 + «41Р4 — С2,

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

Р1 — С1, Р2 — С5, Р3 — Сб, Р4 — С7 «21 — С4/С7,

«31 — (С3С7 - С4С6)/С7, «41 — (с2с2 - С4С5С7 - С3С6С7 + С4с2)/с7. (25)

Теперь перейдем к исследованию численной схемы (24). Для этого разложим стадии 1 < г < 7, в ряды Тейлора в окрестности точки уп до членов с к4 включительно и подставим в первую формулу (24). Полагая уп — у(гп) и сравнивая полученное представ-

ление приближенного решения уга+1 с точным (2), получим условия четвертого порядка точности для (24), т.е.

7

1) £ * = 1,

г=1

2) а(с1 + 2*2 + 3сз + 4*4) + (а + $21 )*5 + (2а + + (3а + Д^К = 0.5,

3) а2(с1 + 3*2 + без + 1ОС4) + (а2 + а$21)с5 + (3а2 + Зав21)сб + (6а2 +

+4ав21)с7 =1/6,

4) в21(е5 + ее + С7) = 1/3, (26)

5) а3(е1 + 4*2 + 10ез + 20*4) + (а3 + 3а2$21)*5 + (4а3 + ба2$21 )ее +

+ (10а3 + 10а2в21)*7 = 1/24,

6) а$21 (*5 + е6 + *7) = 1/8,

7) ав21 (0.5*5 + е6 + 1.5*7) = 1/24,

8) $1^5 + + 07) = 1/4.

Исследуем совместность (26). Из четвертого и шестого уравнений получаем а = 3/8, из четвертого и восьмого имеем в21 = 3/4. Подставляя полученные значения а и в21 в шестое и седьмое уравнения (26), запишем

+ + =16/27, 0.5*5 + + 1.5*7 =16/81. (27)

Пусть есть свободный коэффициент. Тогда, подставляя соотношения (27) в первое, второе, третье и пятое уравнения (26), имеем

1 1 1 1 \ *1 ( +11/27 \

1 2 3 4 *2 -20/81

1 3 6 10 *3 -80/81 - *7

1 4 10 20 / *4 -128/81 - 5*7 /

Из данной линейной системы определим 1 < г < 4, а из (27) — коэффициенты и *6. В результате имеем

а = 3/8, $21 = 3/4, = (60 - 81*7)/81, = (18 - 324*7)/81,

= (405*7 - 64)/81, = (19 - 162*7)/81, = (64 + 81*7)/81, (28) = -(16 + 1626^)/81.

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

Из приведенных выше рассуждений следует, что в классе (т, к)-методов с двумя вычислениями правой части задачи (1) на шаге интегрирования можно построить метод четвертого порядка точности. Однако он не обладает требуемыми свойствами устойчивости, что приводит к очевидным трудностям при решении жестких систем. Оказывается, если по-другому организовать вычисления, то при к = 2 можно построить ¿-устойчивую численную схему четвертого порядка.

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

Выберем множества (8) следующим образом:

Мт = {1, 2,..., т}, Мк = {1, 3}, ^ = {2}, 3 < г < т, т. е. рассмотрим численные схемы вида

т

Уп+1 = Уп + рЛ, = Е - а^/П;

¿=1

ЕпЬ = Л/(Уп), Епк2 = &1, Бпкз = Л/(Уп + вз1 ¿1 + вз2 Ы + «32^2, (29) йпк, = ^¿-1 + «¿2^2, 4 < г < т.

Пусть т = 4. Подставим разложения 1 < г < 4, в виде рядов Тейлора в первую формулу (29). Полагая уп = у(£п) и сравнивая полученное представление с рядом Тейлора для точного решения до членов с Л4 включительно, получим условия четвертого порядка точности схемы (29), т.е.

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

2) ар1 + 2ар2 + (а + А31 + в32 + 3а«32)р3 + (2а + ^31 + +в32 + 4аа32 + 3а«42)р4 = 0.5,

3) а2р1 + 3а2р2 + (а2 + 2а^31 + 3а^32 + 6а2«32)р3 + (3а2 + +3а^31 + 4а^32 + 10а2«32 + 6а2 «42 )р4 = 1/6,

4) а3р1 + 4а3р2 + (а3 + 3а2+ 6а2^32 + 10а3а32)р + (30) + (4а3 + 6а2в31 + 10а2 в32 + 20а3 а32 + 10а3 а42 )р4 = 1/24,

5) (в31 + в32)2(р3 + Р4) = 1/3,

6) а(в31 + в32)(в31 + 2^32)(Р3 + Р4) = 1/8,

7) а(в31 + в32)2(0.5р3 + Р4) = 1/24,

8) (в31 + в32)3(р3 + Р4) = 1/4.

Применяя (29) для решения (17), получим условие ¿-устойчивости, которое имеет вид

а(а - Р1) + (^31 - а)р3 = 0.

Алгебраическая система (30) является нелинейной относительно коэффициентов а^, что приводит к определенным трудностям с ее исследованием. Поэтому при изучении методов (29) поступим по аналогии с изучением численной схемы (22). Для этого рассмотрим следующую вспомогательную схему:

6

Уп+1 = Уп + ^ с^, Дп = Е - аЛ/п,

¿=1

= Л/(Уп), = ¿1, = Л/(Уп + в31^1 + $32^), (31)

Подставляя 1 < г < 6, в (29), получим

Уп+1 = Уп + + Р2 ¿2 + Р3^3 + («32^3 + «42^4)^4 + Р4^5 + «32^4^6.

Отсюда видим, что коэффициенты методов (29) и (31) связаны соотношениями

Р1 = Р2 = *2, Р3 = *3, Р4 = *5, а32 = /*5, «42 = (е4е5 - cзc6)/c2. (32)

Теперь перейдем к изучению метода (31). Подставим разложения 1 < г < 6, в виде рядов Тейлора в первую формулу (31). Полагая = у(£га) и сравнивая полученное представление с рядом Тейлора для точного решения до членов с Л,4 включительно, получим условия четвертого порядка точности схемы (31), т.е.

1) £ ^ = 1,

г=1

2) а*1 + 2а*2 + (а + $31 + в32)*3 + 3а*4 + (2а + $31 + А«^ + 4а*6 = 0.5,

3) а2*1 + 3а2*2 + (а2 + 2ав31 + 3ав32)*3 + 6а2+ (3а2 + 3ав31 + +4а$32 )*5 + Юа2^ = 1/6,

4) ($31 + в32)2(*3 + *5) = 1/3,

5) а(в31 + в32)(в31 + 2в32)(*3 + *5) = 1/8,

6) а(в31 + в32)2(0.5*3 + ) = 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.

Исследуем совместность (33). Из четвертого и седьмого уравнений имеем

в31 + в32 = 3/4, + = 16/27.

(33)

(34)

Тогда из пятого и шестого соотношений (33) получим в32 и соответственно. Из (34) выразим в31 и *5. Учитывая первое уравнение (33), подставим полученные выражения для в31, в32, и во второе, третье и восьмое равенства (33). Получим линейную систему алгебраических уравнений относительно *2, и , которая имеет вид

а*2 + 2а*4 + 3а*6 = -(22а + 5)/54, 2а2*2 + 7а2*4 + 16а2 = (64а2 + 29а - 30)/54, а2*2 + 3а2+ 6а2*6 = (32а2 - 11а - 6)/54.

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

_ 76а2 - 29а + 3 _ -146а2 + 89а - 12 32а - 4 = -2-, = --, -

27а2 27а2

_ 72а2 - 59а + 10 _ 4 - 16а

= -п о 2-, = -^-,

в31

18а2

48а - 9 9 - 24а

в32 =

27а

27а

-18а2 + 19а - 4 18а2 ,

32а 32а

Применяя (31) для решения (17), получим условие ¿-устойчивости

6

а(*1 + *3) - а2 - *3в31 = 0.

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

24а4 - 96а3 + 72а2 - 16а +1 = 0. (36)

Заметим снова, что (36) есть известное уравнение ¿-устойчивости четырехстадийных схем, т. е. при одинаковом числе стадий функции устойчивости методов типа Розенброка и (т, к)-схем совпадают. Подставляя теперь (35) в (32), получим коэффициенты схемы (29), т.е.

76а2 - 29а + 3 146а2 + 89а - 12 32а - 4

Р1 =-27а-' Р2 =-27а2-' Р3 = ^7^'

4 - 16а д 48а - 9 д 9 - 24а

Р4 = ^7^' в31 = ^2^' в32 = ^2^ (37)

= -54а2 + 57а - 12 = -864а3 + 828а2 - 288а + 36 а32 = 8а - 32а2 ' "42 = а(4 - 16а)2 '

при которых она имеет четвертый порядок точности. Коэффициент а определяется из условия ¿-устойчивости (36). Данное уравнение имеет корни а1 = 0.10643879214266, а2 = 0.22042841025921, а3 = 0.57281606248213 и а4 = 3.10031673511599. Для расчетов рекомендуется а = 0.57281606248213, а соответствующие коэффициенты имеют вид

а= +0.57281606248213, Р1 = +1.27836939012447,

Р2 = -1.00738680980438, Р3 = +0.92655391093950,

Р4 = -0.33396131834691, в 1 +1.00900469029922,

в32 = -0.25900469029921, «32 = -0.49552206416578,

«42 = -1.28777648233922.

Заключение

Из приведенных выше рассуждений можно сделать следующие выводы.

Во-первых, свойства устойчивости (т, к)-методов зависят от выбора множеств (8), или, что то же самое, от способа реализации численных схем. Это следует из сравнения численных формул (13) и (29).

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

В-третьих, при двух вычислениях функции f (у) задачи (1) можно построить ¿-устойчивый (4, 2)-метод четвертого порядка точности, в то время как соответствующий ¿-устойчивый метод типа Розенброка (3) может быть только второго порядка

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

Список литературы

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

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

[3] Новиков Е.А., Новиков В.А., Юматова Л.А. Замораживание матрицы Якоби в методах типа Розенброка второго порядка точности // ЖВМ и МФ. 1987. Т. 27, № 3. С. 385-390.

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

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

Поступила в редакцию 18 апреля 2007 г.