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

Вычислительные технологии

Том 9, № 2, 2004

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

И. В. ОЛЕМСКОЙ Санкт-Петербургский государственный университет, Россия

e-mail: OLIV@pobox.spbu.ru

An explicit one-step method for the numerical integration of a special type of ordinary differential equations is considered. Some statements on the relation between the order of method and the number of stages are proved. Economical numerical schemes up to the fourth order of accuracy are constructed for the integration of systems and equations of the second order.

Введение

При решении задачи Коши

d: = i = 1,..., (1)

У:(Хо) = y:o, i = 1,... , n; (2)

x G [X0,X1] С R, y: : [X0,X1] —► R, i = 1,...,n; /: :[Xo,Xi] x Rn —► R, i = 1,...,n,

общая схема явного метода Рунге — Кутты [1-3] численного интегрирования системы (1) имеет вид

mi

У:(х + h) ^ Z: = y:(x) + ^ bjkj(h), i = 1,..., n, (3)

j=i

где функции k:j (h) вычисляются по следующей схеме:

j-i j-i kij (h) = h/:(x + Cij h,yi(x) + ^ ^ aijipkip) . . . ) Уп (x) + ^ ^ aijnp knp) ; (4)

p=i P=i

c:i = 0, a:iso = 0, s = 1,... ,n.

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

Здесь шг — число этапов и — порядок точности по г-й компоненте искомой функции в общем случае могут быть различны по каждой из компонент. Поэтому при построении расчетных схем метода их характеристикой могут выступать как векторы числа этапов М = (ш1 ,... , шп) и порядка точности Q = (д1,... , qn), если хотя бы одна из компонент соответствующих векторов отлична от остальных (и тогда метод называется разноэтапным и разнопорядковым соответственно), так и скаляры ш и д, если компоненты соответствующих векторов равны между собой: Ш1 = ... = шп = ш, = ... = дп = q.

Под порядком точности метода (3), (4) при векторе порядка точности Q = (д1,... , дп) будем понимать д = ш1п(д1,... , дп}.

Введем еще одно понятие, необходимое для сравнения методов интегрирования систем обыкновенных дифференциальных уравнений.

Определение. Пусть для некоторого вектора числа этапов в рамках одношагово-го метода существует расчетная схема с вектором порядка точности Q = (д1,... , дп). Будем называть такой вектор числа этапов минимальным и обозначать М= (ш1(д1),... ,шп(дп)), если для любой М = (ш1, ... ,шп)-этапной расчетной схемы с вектором порядка точности Q = (д1,... ,дп) этого же метода справедливы неравенства Ш > шг(дг), г = 1,... ,п.

В скалярном случае минимальное число этапов метода порядка точности д будем обозначать ш(д). При построении расчетных схем [1-4] в рамках метода (3), (4) традиционно в силу равноправности уравнений системы параметры метода полагают не зависящими от номера компоненты г и в:

Такой способ распространения метода интегрирования (3), (4) уравнений на системы, безусловно, сужает его возможности, однако значительно упрощает задачу конструирования методов интегрирования систем. В этом случае методы строятся для скалярных уравнений, а распространение на системы осуществляется простой заменой скалярных функций у, /, к^ на соответствующие векторные. При этом все утверждения о соотношении порядка точности метода д и числа этапов ш, справедливые для скалярного случая, верны и для векторного. Такой способ распространения и методы, построенные с его использованием, будем называть формальными, причем известно [4], что для явного метода Рунге — Кутты д-го порядка точности численного интегрирования дифференциального уравнения минимальное число этапов удовлетворяет равенствам

Заметим, что равенства (5) справедливы не только для формального метода интегрирования системы (1), но и для общей схемы метода (3), (4).

Известная теория условий порядка для скалярных уравнений была обобщена Хайрером [5, 6] на случай разделяющихся систем вида

ш(д) = д, д < 4, ш(5) = 6.

(5)

= /1(У1,У2), % = /2(У1,У2)

для неявного одношагового метода интегрирования

т

Уг(ж + Л) ^ уг(ж) + ^ Ьу к^, % = 1, 2, (7)

¿=1

тт

к^ = Л/г(у1(ж) + ^ «¿1Рк1Р, У2(х) + ^ а^2Рк2Р), % = 1, 2, ^ = 1,..., т.

р=1 р=1

Параметры метода интегрирования (7) разделяющихся систем зависят от номера компоненты решения 5, но все они одинаковы для всех компонент % правой части:

Ьщ = Ьщ , Cj = Су, = aijsp, т' = т,г.

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

В [7, 8] выделен класс систем обыкновенных дифференциальных уравнений:

^ = Л^уо^...^; (8)

аж

= /г(ж,Уо,.. . ,Уг-1,Уд+1; •••,Уn), % = (9)

оу

ож

7

у,-(ж,уо,...,У-1), = £ + 1,...,п (10)

где у., /s — функции размерности г.. Две группы уравнений (9), (10) структурно тождественны. Каждое уравнение одной из групп уравнений (9), (10) занимает определенное место в последовательности уравнений своей группы. Его правая часть не зависит от искомых функций, поведение которых описывается этим и всеми последующими уравнениями этой же группы. Группа уравнений (8), в которую вошли все уравнения, не имеющие структурных особенностей указанного выше типа, называется общей. Она, как и группа уравнений (9), может отсутствовать.

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

Для численного интегрирования систем (8)-(10) в [7, 8] предложен явный одношаговый структурный метод. В отличие от метода (3), (4) его вычислительная схема использует структурные особенности на уровне алгоритма:

т3

у.(ж + Л) ^ у.(ж) + ^ к.,(Л-), 5 = 0,..., п, (11)

¿=1

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

к0Ъ к1Ъ . . . , kn1, к02, . . . ; (12)

к.,(Л) = ,^¿ь ... ,У7п); 5 = 0,... , п, (13)

где

0, если в < д, > 0, если в > д;

Уя

если ] = 1 и в < д,

7-1

sjd

Уd(x) + Е asjdpkdp(h), если 1 и в <

р=1

j

Уd(x) + Е asjdpkdp(h), если = 1 и в > д или > 1 и в > р=1

(14)

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

Благодаря этому обстоятельству компоненты минимального вектора числа этапов М структурного метода ^-го порядка численного интегрирования системы (8)-(10) удовлетворяют равенствам [8]

то(<?о) = <?0, Шг(^г) = ^ - 1, г = 1,...,п, & = 3, 4, в = 0,...,п.

Отсутствие общей группы уравнений (8) (г0 = 0) в рассматриваемой системе может еще больше повысить эффективность метода (11 )-(14). Так, например, в [9-11] для систем

«Х = /1(х,у2), = /2(х,у1)

(15)

(у, /г — функции размерности гг, г = 1, 2) рассмотрен следующий частный случай структурного метода (11)—(14). Приближенное решение ищем в виде

Уг(х + Л,) ^ ^ = Уг(ж) + ^ Ь. kij (Л,), г = 1, 2, Ш1 > Ш2,

j=l

где функции к. = к. (h) вычисляем в строгой последовательности

к1Ъ к21, к12, к22, . . .

(16)

(17)

по схеме

к. =

^ = 1,

.7-1

Л./1(ж + с.h, У2(х) + Е aljnк2Ч), > 1;

П=1

(18)

k2j = Л/2(ж + С27h, У1(х) + ^ ^ a2jn к1п ) , С21 = 0,

П=1

(19)

и в его рамках построены [11] четырехэтапные расчетные схемы пятого порядка, которые на треть экономичнее формального метода Рунге — Кутты.

т

Т аб л и ц а 1

М = (т1, Ш2)-этапный структурный метод (16)—(19) интегрирования системы (15)

0

С12 <121 С21 <211

с1тх <1тх1 . . . <11т1т1 — 1 с2т2 <<2т21 • • • <2т2т2

Ь11 ••• Ь1тх-1 Ь1тх Ь21 • • • Ь2т2 Таблица2

Четырехэтапный метод пятого порядка

0 <1г3

4^6 15 4^6 15

8 36^9^6 128 28^7^6 128

14±^6 20 9954^4419^6 18800 26556±9809^6 112800 -1464±4469\/6 22560

С13/Ь13 328±77^6 1140 -1188^2457^6 5348 2432±448^6 -18184±51676^6 2415 250401

30

10

10 1

<2г3

4^л/б 30

4^6 12^3^6

40 40

5348±1947-Уб -18204^7581^6 8448±3472^6

5000 7000 4375

-206^83^6 40614±11721^6 -504^272^6 -62874±49236^6

76 5348 161 83467

С23/Ь23

36

16±\/6 36

С использованием таблиц Батчера [4,5] по каждой компоненте метод (16)-(19) можно представить в компактной форме (табл. 1).

Так, полученная в [11] в рамках метода (16)—(19) четырехэтапная расчетная схема пятого порядка имеет табличное представление (табл. 2).

Известно [5], что формальный метод Рунге — Кутты при интегрировании дифференциальных уравнений второго порядка

У = /УЫ = У0, У Ы = у0 (20)

уступает по своим характеристикам методу Нюстрема (прямому методу Рунге — Кутты):

1-V

3=0

+ й) « ^ = ]>] й"у(з+^(х) + Д^, V = 0,1,

(21)

г=1

г-1

к = й/(X + Сгй, у(х) + Сгйу0 + й ^ Ац к ), С = 0.

3=1

Для компактной записи явного т-этапного метода Нюстрема используется [5,12] табличное представление (табл. 3). После приведения структурного метода (16)—(19) к виду

0

т

экономичному и алгоритмически простому для непосредственного интегрирования дифференциальных уравнений (20) его расчетная схема примет вид

т1 — 1

у(ж + Л) ^ = у(ж) + Лу' (ж) + Л В0гк,

г=1

т 2

у' (ж + Л) « = у' (ж) + ^ Вцкг, (22)

г=1

г—1

к2г = кг = /ж + СгЛ, у(ж) + СгЛу'(ж) + Л ^ Агзк,) ,

3=1

т1 г

В0г = ^^ Ь13а 13г) В1г = Ь2г, Аг3 = ^^ а2г4а13, Сг = с2г-

3=г+1 ¿=3+1

Так, метод пятого порядка (см. табл. 2), приведенный к виду (22) для прямого интегрирования уравнения (20), имеет вид, представленный в табл. 4.

Надо отметить, что, во-первых, упомянутый выше (см. табл. 4) метод согласно классификации в [5] принадлежит классу методов Нюстрема, причем по своим характеристикам (числу этапов и порядку точности) он совпадает с лучшими известными; во-вторых, рассматриваемый подход дает ответ на вопросы, какие особенности должна иметь интегрируемая система и каким образом их должен использовать метод, чтобы при применении метода к дифференциальному уравнению у (ж) = /(ж, у(ж)) он не уступал методу Нюстрема?

Цель данной работы — дальнейшее исследование возможностей метода (16)—(19) интегрирования системы (15) (при т1 = 4, т2 = 3). Заметим,что в рамках трехэтапного метода (16)-(19) в [8] получены схемы четвертого порядка точности, а в [11] — для четырехэтап-ного (см. табл. 2) — пятого. С одной стороны (чисто формальной), необходимо ответить на

ТаблицаЗ т-этапный метод Нюстрема

С г Аг3 Воз В13

С1 0 В01 В11

С2 А-21 В02 В12

Ст Ат1 Ат2 • • • Атт— 1 0 В0т В1 т

Таблица4 Метод пятого порядка

Сг Аг3 В03 В13

30 10 10 1 0 11^4^6 100 13±7^6 81±34^6 250 500 4±5^6 -,-^6 12^3^6 16 + 8 16 0 0 36 36 9^6 16±^6 36 36 0 9

вопросы: каким образом увеличение числа этапов по первой компоненте повлияет на характеристики метода? Удастся ли повысить порядок точности хотя бы по первой компоненте на единицу? С другой стороны, при практическом интегрировании систем (15) довольно часто возникают ситуации неравнозначности трудозатрат при вычислении правых частей /¿, г = 1, 2. В этом случае вполне оправдана попытка улучшения характеристик метода за счет дополнительного вычисления правой части, требующей меньшего количества арифметических операций, так как затраты на вычисление правых частей (относительные) в этом случае возрастут ненамного, а выгода очевидна.

И, наконец, с точки зрения здравого смысла в схеме (22) при вычислении приближенного значения искомой функции и ее производной естественно использовать одинаковое число &2г(Л-), что означает выполнение равенства т2 = т1 — 1. А это еще один аргумент в пользу рассмотрения разноэтапных методов.

1. Разноэтапный метод четвертого порядка

Теорема 1. В рамках структурного метода (16)-(19) при М = (т1,т2) = (4, 3) численного интегрирования системы (15) не существует расчетной схемы порядка точности

Я =(?1,?2) = (5, 4).

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

т

ФДй) = уДж + й) — уДж) — ^ ЪцЬз (й), г = 1, 2,

3 = 1

метода (16)—(19) по степеням й начиналось с шестой степени при числе этапов т1 = 4, т2 = 3. Условия порядка метода (16)—(19) в этом случае с использованием упрощающих предположений

я

^а^ = Csq ,5 = 1, 2; си = 0, й1яя = 0,5 =1,..., 4, 3=1

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

]>>я = 1, (*,«) = {(1, 2), (2,1)}; (23)

я=1

^ 1 2;

с^ = -; (24)

я=1

т.

' я С'„„

1 3

я=1

т. я

а„,.А с...А =

6

!>я с2я = 3; (25)

^ у а^я3с-и3 = 6 ; (26)

Я=2 3=1

т

8

та

с39 = 4; (27)

9=1

т3 9 1

^ ] ^ ] 8; (28)

8

9=2 3=1

2 1

а«93С23 = 1п) (29)

12

9=2 3=1

т3 д 3

У^ X ^3 X/ = 24 ; (30)

9=2 3=2 ¿=1 1

т3

с^ =-; (31)

9=1

^ а«93С^3 = 10 ; (32)

10

9=2 3=1

т3 9

21

^С23 = 1- ; (33)

9=2 3=1

т3 9 3 1

^^ ] = (34)

9=2 3=2 ¿=1

т3 /9 \ 2 1

Х/^ I X/ а«93'С^3 I = 2о; (35)

9=2 \3=1 /

т3 9 1

ХХХ ^934/ = 20; (36)

9=2 3=1 20

3

9=2 3=2 ¿=1

т3 9 3

^ ] ^93 СУ3 ^ у 40 ; (37)

1

XX X = 60; (38)

9=2 3=2 ¿=1

9 3 *

1

^ ] ^93 ^ у у ^р^ир 120; (39)

9=3 3=2 ¿=2 р=1

г

У^ а^ = с«; г = 1,. • •, 4. (40)

3=1

Такая компактная запись условий порядка вынуждает при ссылке на некоторое равенство

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

через точку с запятой значением параметра 5. Например, ссылка (27; 2) подразумевает равенство (27) при значениях пары (5, V) = (2,1), а (27; 1) — равенство (27) при значениях пары (5, V) = (1, 2).

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

алгебраических уравнений (23)-(30), (31; 1)-(39; 1), (40) противоречива. Допустим, что си-

9

т

3

9

т

т

стема совместна при выполнении условия

С22С23(С22 — С2з)(с22 — С21)(С23 — С21) = 0. (41)

Преобразуя уравнения (26; 1), (29; 1) с помощью соотношений (24; 1), (40; 1), получим равенства

, . , 1 — 2(с21 + с23) + 6с21с23 ,ло,

Ъ13Й132 + Ъ14Й142 = — ^^-ГТ-Г-; (42)

12(С23 — С22)(С22 — С21)

1 — 2(С21 + С22) + 6С21С22

Ъ14Я143 = —Г(--Г". (43)

12(С23 — С21)(С23 — С22)

С другой стороны, из рассмотрения уравнений (23; 2)-(25; 2), (27; 2) имеем

Ъ22(1 — с22) = — 1 — 2(с21 + с23) + 6С21С23; (44)

12(с23 — с22 )(с22 — с21)

Ъ23(1 — С23) = 1 — 2(С21 + С22) + 6с21с22 . (45) 12(с23 — с21) (с23 — с22)

Сравнение равенств (42)-(45) с учетом (26; 1) приводит к соотношениям

Ъ^Яш + Ъ13Й131 + Ъ14Й141 = Ъ21(1 — С21),

Ъ13Й132 + Ъ14Я142 = Ъ22 (1 — С22), (46)

Ъ14Я143 = Ъ23 (1 — с23) ,

которые в теории условий порядка [4, 5] известны как упрощающие.

Дальнейшее рассмотрение системы, в частности преобразование соотношений (36; 1), (27; 2) с учетом (46)

4 г— 1 34 3 1 3 1

53 Я1г34 = 53 с33 53 Ъ1га1г3 = 53 4Ъ23 (1 — С23) = 4 — 53 Ъ23с43 =

г=2 3=1 3=1 г=3+1 3=1 3=1

приводит к равенству

Ъ21С21 + Ъ22С^2 + Ъ23С43 = 1. (47)

5

Аналогично из соотношений (37; 1), (28; 2), (46) следует ограничение

Ъ22с22 Я222С12 + Ъ23С23(Я232С12 + Я233С13) = . (48)

Уравнения (26; 2), (28; 2), (48) непротиворечивы, если параметры метода удовлетворяют соотношению

1 1 1

^С22С23 — -: (С22 + С23) + 3 = 0. (49)

3 4 5

В то же время при выполнении ограничения (49) для совместности уравнений (23; 2)-(25; 2), (27; 2), (47) требуется выполнение равенств

11

С22С23 — ^ (С22 + С23) + 3 = 0; (50)

1С22С23 — 1(с22 + С23) + 1 = 0. (51)

Но это невозможно, так как система уравнений (49)-(51) противоречива.

Это и доказывает несовместность исходной системы (23)-(30), (31; 1)-(39; 1), (40), а значит, и утверждение теоремы.

Осталось показать справедливость утверждения теоремы в случае невыполнения неравенства (41).

Предположим, что рассматриваемая система (23)-(30), (31; 1)—(39; 1), (40) непротиворечива при с23 = 0. С учетом того, что с21 = 0, 623 = 0, 614 = 0, с12 = 0, равенство (28; 2) примет вид

&22С22Й222 С12 = 1. (52)

8

Уравнения (24; 2), (25; 2), (27; 2), (26; 1), (29; 1), (36; 1) совместны при

3 3 (53)

с21 = - Т-, с22 = - ±-. (53)

21 5 10 ' 22 5 1 0 ^ '

Рассмотрение на непротиворечивость уравнений системы (30; 1), (37; 1) и (26; 2) (с учетом (28; 2)) приводит к ограничению

&14а143 =623 ¿¡Й-, (54) а уравнений (30; 1), (38; 1), (26; 2), (29; 2) — к соотношению

6-3 = 623 ^ ■ (55)

Очевидно, что последние два равенства (54), (55) справедливы при условии

2

С12 = 3 С22. (56)

С учетом того, что с12 = 0, соотношения (24; 1), (25; 1), (27; 1), (31; 1), (53), (56) могут быть непротиворечивы только при выполнении неравенства

С13С14(С13 - С14)(С12 - С14)(С12 - С13) = 0. (57)

Полагаем, что это так. Тогда одно из условий совместности уравнений (25; 1)-(28; 1), (31; 1), (32; 1) имеет вид

С12 = 2с21. (58)

Но соотношения (53), (56), (58) противоречивы. Это и доказывает утверждение теоремы при С23 = 0.

Совершенно аналогично можно показать справедливость утверждения теоремы и при выполнении любого из нерассмотренных случаев: с22 = 0, с21 = с22, с21 = с23, с22 = с23.

Таким образом, в силу доказанного утверждения теоремы увеличение количества вычислений первой правой части /1 на единицу не позволяет повысить также порядок точности по первой компоненте на единицу.

Какие все-таки преимущества можно извлечь из четвертого вычисления правой части первого уравнения исходной системы? Покажем, что дополнительное вычисление правой части позволяет построить экономичный метод четвертого порядка точности с контролем погрешности на шаге интегрирования.

2. Вложенные разноэтапные методы

Для записи вложенных методов используем уже устоявшуюся [4, 5] терминологию и их табличное представление, несколько видоизменив их (см. табл. 5).

Параметры М = (т1, т2)-этапного структурного вложенного метода (16)—(19) с точностной характеристикой Р[ф] = (рър2)[(9ъ92)] интегрирования системы (15) должны обеспечивать для величин

Ш1 т2

¿1 = У1(ж) + ^ Ьц^2 = У2(х) + 53 ^^

г=1 г=1

порядки точности р1 и р2 (приближений к решению)

У1(х + й) - = 0(йР1+1), У2(х + й) - ^2 = 0(йР2 + 1), а для ("оценщиков" погрешности)

Ш1 т2

¿1 = У1(х) + Ьн кц, ¿2 = У2(х) + У^ Ь2г^2г

г=1 г=1

— 91 и 92, соответственно

У1(х + й) - ¿1 = 0(й"1+1), У2(х + й) - ¿2 = 0(й92 + 1).

Здесь в рамках структурного подхода (16)—(19) с числом этапов М = (4, 3) для интегрирования системы (15) строим вложенный метод с точностной характеристикой Р[ф] = (4, 4)[(3, 2)], т.е. основной метод имеет четвертый порядок точности, а "оценщик" по первой компоненте — третий, а по второй — второй. Используя упрощающие предположения

53 Ьу ау = Ь2i (1 - С2г) , г = 1, . . . , Ш2, ^ а^у= 2С1д,

3=г+1 3=1

выпишем условия порядка, которым должны удовлетворять параметры вложенного мето-

Т аб л и ц а 5 М = (т1, т2)-этапный вложенный Р[ф] = (р1,р2)[(д1, )] метод (16)—(19) интегрирования системы (15)

0

С12 а121 С21 а211

с1т1 а1,т1,1 • • а1,Ш1,Ш1 — 1 с2т2 а2т21 • • а2т2т2

¿1 Ь11 • • Ъ1,т1-1 Ъ1,т1 ¿2 Ъ21 • • Ъ2т2

¿1 Ъи • • Ъ1т1-1 Ь1т1 ¿2 Ъ 21 • • Ъ2т2

да (см. табл. 5) с такими характеристиками:

4

:>>=1,

г=1 4

= -+Г, 3 = 1, 2, 3

г=2 3 +

3

Х^а1,3+1,г = С1,3+1, 3 = 1 2, 3,

г=1

3 1 Ха1,3+1,гс2г = 2с1,3+1, 3 = 1, 2, 3,

г=1

4

X Ьцащ = 623(1 - С23), 3 = 1, 2, 3,

г=3+1

3

5>г4 = 3+1, 3 = 0, 1, 2, 3,

г=1 3 +

3 г 1

а2г3с13 = 6, (59)

г=2 3=2 6

3 г 1

X Ь2гс2г ^ «2г3с13 = 8,

г=2 3=2 8

3 г 1

г=2 3=2 12

3

Ха23г = С23, 3 = 1, 2, 3,

г=1

4

£61г =1,

г=1

41 Е^г = —у, 3 = 1, 2,

г=2

31

= 7-т, 3 = 0, 1

г=1 3

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

С14 = 1, Й141 = 621, «142 = 622, «143 = 623. (60)

Решения системы (59), (60) образуют двухпараметрическое семейство. Не выписывая общее, приведем два частных ее решения. Для первого доопределим два параметра, положив с21 = 1/6, с22 = 1/2. Расчетная схема, соответствующая этому решению, приведена в

Таблицаб

М = (4, 3)-этапный вложенный Р[ф] = (4, 4)[(3, 2)] метод (16)-(19)

0

1 3 1 3

1 2 1 3 8 3 8 1 8 1 4 3 8

¿1 1 6 0 2 3 1 6

¿1 1 2 3 2 2 0

1 1

6 6

1 0 5 1 21 8

6 18 3 9

¿2 С0100 1 4 3 8

¿2 1 2 0 1 2

Таблица7

М = (4, 3)-этапный вложенный Р[ф] = (4, 4)[(3, 2)] метод (16)-(19)

0

1 ■/15 5 1 /15 5

13 8 /15 4 1 45 16 125/15 192 5 18 19 77/15 16 192 4 9 5 18

¿1 281 1308 49/15 1308 25 , 115/15 1356 + 1356 9856 36951 , 9728/15 + 184755 1 /15 2 10

¿1 569 1308 55/15 1308 415 125/15 1356 1356 9568 36951 , 4960/15 + 36951 0

1 2 /15 10 1 2 /15 10 а2у

1 2 1 32 /15 32 15 32 1 /15 + 32

1 2 1 /15 + 10 487 2180 , 33/15 + 2180 75 452 83/15 2260 6816 , 7488/15 61585 1 61585

¿2 5 18 4 9 5 18

¿2 1 2 0 1 2

табл. 6. Второе частное решение (табл. 7) базируется на трехточечной квадратурной формуле Гаусса — Лежандра, так как в качестве с2^ взяты нули сдвинутого полинома Лежандра третьей степени:

_ 1 л/15 1 1 л/15

С21 = 2 10"' С22 =2' С23 = 2 + 10" •

Характеризуя полученные схемы (см. табл. 6, 7) интегрирования, естественно использовать для сравнения известные аналоги — классическую четырехэтапную схему метода Рунге — Кутты с контрольным членом Егорова 4(2) [13] и пятиэтапную расчетную схему Зонневельда 4(3) [5]. При равном порядке точности с указанными оппонентами построенные расченные схемы более эффективны, так как по сути (в силу выполнения равенств (60)) являются трехэтапными. Подчеркнем два важных момента (не приводя здесь их доказательства): во-первых, в рамках структурного подхода (16)—(19) с числом этапов М = (4, 3) для интегрирования системы (15) нельзя построить вложенный метод с точностной характеристикой Р[ф] = (4,4)[(3, 3)]; во-вторых, не существует трехэтапной расчетной схемы вложенного метода (16)-(19) с точностной характеристикой Р[ф] = (4, 4)[(3, 2)].

Для представления вложенных методов интегрирования дифференциального уравнения (20) воспользуемся табличной формой (табл. 8), приведенной в [14]. Здесь , С1, А13- — параметры метода (21) Р = (р0,Р1)-го порядка точности вычисления приближения ZV, V = 0,1, к решению + Л) и его "оценщиков" — Я = (д0, ^1)-го порядка

1-V

Лз '

Zv = ^ -У(з+')(х) + Л1-^ Д^, V = 0,1, 3=0 -' 1=1

у(х + Л) - = О(Л90+1), у'(ж + Л) - = 0(Л91+1).

Приведенные к виду (22) для прямого интегрирования уравнения (20) расчетные схемы (см. табл. 6, 7) представлены в виде табл. 9,10 соответственно.

Существующий аналог [14] — пятиэтапный метод И,КМ4(3)4ЕМ уступает трехэтапным расчетным схемам (см. табл. 9,10), так как при одинаковом порядке точности требует четырех этапов в случае успеха. В случае же неудачи в методе ИКМ4(3)4ЕМ приходится отказываться от четырех вычислений правой части, в то время как для полученных здесь схем — от трех.

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

Таблица8 т-этапный вложенный метод типа Нюстрема

С А3 В01 В01 В11 В11

С1 0 В01 В01 В11 В11

С2 А21 В02 В02 В12 В12

Ст Ат1 Ат2 • • • Атт-1 0 В0т В0т В1 т В1т

Таблица9

Трехэтапный вложенный Р= (4, 4)[(3, 2)] метод (21)

С

"X"

6 1 2

5

6

А.

3

B0i

16 1 8 1 16

В01

4 1 4

0

В11

8 1

4 3 8

В11

2 0

1

2

Таблица 10 Трехэтапный вложенный Р= (4, 4)[(3, 2)] метод (21)

0

С1 А3 В01 В01 В11 В11

1 2 -15 10 0 5 36 + ^15 + 36 -15 18 5 18 1 2

1 2 3 -15 8 16 2 9 1 -15 2 18 4 9 0

1 2 + ^15 + 10 3 + -15 5 + 5 3 5 -15 10 5 36 -15 36 0 5 18 1 2

му и алгоритмически простому для непосредственного интегрирования дифференциального уравнения y (x) = f (x,y(x)) попадают в класс методов типа Нюстрема (21), делая

ограничение ci = 0 необязательным.

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

[1] Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Наука, 1987.

[2] Березин И.С., Жидков Н.П. Методы вычислений. Т. 2. М.: Физматгиз, 1966.

[3] Крылов В.И., Бобков В.В., МонАСТЫРНЫй П.И. Вычислительные методы высшей математики. Т. 2. Минск, 1975.

[4] Современные численные методы решения обыкновенных дифференциальных уравнений / Под ред. Дж. Холла, Дж. Уатта. М., 1979.

[5] ХАйРЕР Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М., 1990.

[6] HAIRER E. Order conditions for numerical methods for partitioned ordinary differential equations // Numer. Math. 1981. Vol. 36. P. 431-445.

[7] ОлЕМСкоЙ И.В. Численный метод интегрирования систем обыкновенных дифференциальных уравнений // Мат. методы анализа управляемых процессов. Л., 1986. C. 157-160.

[8] ОлЕМСкоЙ И.В. Структурный подход в задаче конструирования явных одношаговых методов // Журн. вычисл. математики и мат. физики. 2003. Т. 43, № 7. C. 961-974.

[9] ОлЕМСкоЙ И.В. О численном методе интегрирования систем обыкновенных дифференциальных уравнений // Оптимальное управление механ. системами. Л., 1983. С. 178-185.

[10] ОлЕМСкоЙ И.В. Экономичная расчетная схема четвертого порядка точности численного интегрирования систем специального вида // Процессы управления и устойчивость: Тр. XXX научн. конф. СПб.: НИИ химии СПбГУ, 1999. С. 134-143.

[11] ОлЕМСкоЙ И.В. Четырехэтапный метод пятого порядка точности численного интегрирования систем специального вида // Журн. вычисл. математики и мат. физики. 2002. Т. 42, № 8. С. 1179-1190.

[12] КоллАТЦ Л. Численные методы решения диференциальных уравнений. М.: Иностр. лит., 1953.

[13] Арушанян О.Б., Залеткин С.Ф. Численное решение обыкновенных диференциальных уравнений на Фортране. М.: Изд-во МГУ, 1990.

[14] DüRMAND J.R., El-Mikkawy M.E.A., PRINCE P.J. Families of Runge — Kutta — Nystrom Formulae // IMA J. Numer. Anal. 1987. Vol. 7. P. 235-250.

Поступила в редакцию 28 ноября 2003 г.