Вложенные методы пятого порядка Текст научной статьи по специальности «Математика»

УДК 519.62/.642 Вестник СПбГУ. Сер. 10, 2004, вып. 2

И. В. Олемской

ВЛОЖЕННЫЕ МЕТОДЫ ПЯТОГО ПОРЯДКА

1. Введение. В работах [1-3] для численного интегрирования систем обыкновенных дифференциальных уравнений

(1)

(у,-, /,- - функции размерности г,-, г = 1,2) в рамках структурного подхода [4] был рассмотрен одношаговый явный метод типа Рунге-Кутты. Приближенное решение ищется в виде

771 {

у,-(ж + Л) и г,' = + ^ » = 1, 2, (2)

3=1

где т,- - число этапов по г'-й компоненте (т1 > ггаг); функции = вычисляются

в строгой последовательности &21, &12, ^22, •• • по'схеме

_ Г ЛЛ(х,у2(г)); ¿ = 1,

у" I ЛЛ(® + суА(№(®) + Е;:}ау|,Л2,), ;>1,

]

Ь] - Л/2(® + с2]к,ух{х) + ^аг^щ), С21 Ф 0. (4)

77=1

В [3] построены четырехэтапные расчетные схемы пятого порядка, которые на треть экономичнее схем метода Рунге-Кутты того же порядка.

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

У" = У), У{х о) = уо, у'(х о) = Уо> (5)

его вычислительную схему запишем так:

7711—1

у(х + Л) » - у(я) + /1у'(ж) + Л Аи*.-,

1=1

т2 ¿=1

¿-1

= Ь = hf(x + СхЬ,,у{х) + С,-Ау'(г) + НУ^А^к^),

}=1

ТП\ I

-Во» = ^ ВЦ = ¿>2г") -А,^- = ^ а,2ааи], С,-= сг;.

¿=» + 1 + 1

И. В. Олемской, 2004

Все расчетные схемы этого типа, полученные в [5], не уступают по своим характеристикам (соотношение порядка точности & и числа этапов ггц) известному методу Нюстрёма [6] прямого интегрирования дифференциального уравнения (5).

Практическая реализация преддаженных. в ЭД щштшх. жж татадатш трёбует для их использования эффективного алгоритма автоматического выбора шага. Поэтому здесь в рамках структурного подхода (2)-(4) на базе расчетных схем [3] строятся вложенные методы типа Дормана-Принса [6,7]. Для них характерно, что члены погрешности для результата старшего порядка минимизированы, сам результат принимается в качестве искомого приближения, а результат младшего порядка используется в алгоритме управления величиной шага.

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

(табл.1). _ , , ,п\ /л\

Таблица 1. то-Этапныи вложенный р^) метод (2)-(4)

Сц а щ

* С12 ^121

С17П1 «1,7111,1 ••• &1, »711,7711 — 1

Ьц Ьц Ьх^тпх-Х

йи ¿11

С2г О-Из

С21 0211

С2 ,т2 02,т2Д • • • 7712,7712

&2г &21 &2,т2

¿2* (¿21 ■ ^2,7712

Параметры тп-этапного вложенного р^) метода (2)-(4) интегрирования системы (1) должны обеспечивать для величин

т

2] - Уз(х) + ^Ьцкц, у3{х + К) - zi = 0{ЪР+1), 3 - 1,2,

г=1

порядок точности р (приближений к решению), а для (оценщиков «погрешности»)

т

Ь ~ Уз (х) + X) (х + К)- Щ = 0(№+1), = 1,2,

¿=1

порядок д.

Наиболее популярный из известных вложенных методов пятого порядка - семиэтап-ный метод Дормана-Принса [7] с оценщиком четвертого порядка - ООРШ5(4)7Г. Его последний этап на текущем шаге является первым на следующем (идея [7,8] - ЕБАЬ). Аббревиатура метода сообщает о типе метода (БОРИ - Дормана-Принса), порядке приближения («5») и оценщика («4»), числе этапов («7»), а также об использовании идеи РБАЬ (Р).

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

Вложенные методы типа Нюстрёма интегрирования дифференциального уравнения (5) будем так же, как в [8], представлять в табличной форме (табл. 2).

Таблица 2. т- Этапный вложенный метод типа Нюстрёма

Сг Ац

Сх 0

с2 ¿21

Сщ ¿тД Ат> 2 Ат, т—1 0

Вы Во1 #02 Во,т-1 Во>т

Ан Dol ^>02 А),т-1 ,т

Ви Вп В\2 В\,т-\

#12

Здесь Виг, С{, Ац - параметры метода (6) р-го порядка точности вычисления приближения и = 0,1, к решению у^\х + К) и его оценщиков - д-го порядка:

" из т

= £ + Л1-" ]Г V. = 0,1,

¿=0 3' 1=1

у(х + Л) - 2о = 0(Л9+1), у'(х + к)-2! = 0(М+1).

В [8] получены вложенные методы типа Нюстрёма с четвертого порядка по шестой, причем метод четвертого порядка (ВР4(3)4Г) требует четырех этапов, а шестого порядка (БР6(4)Е) - шести.

В результате приведения полученной РС5(3)5Г к виду, экономичному и алгоритмически простому (6) для интегрирования дифференциального уравнения второго порядка (5), получены две пятиэтапные расчетные схемы вложенного метода (РС1Н5(3)5Г, РС2Н5(3)5Е) типа Нюстрёма пятого порядка с оценщиком третьего порядка. Последний этап этих схем на текущем шаге является первым на следующем. 2. Расчетные схемы пятого порядка.

Теорема 1. В рамках пятиэтапного структурного метода (2)-(4) численного интегрирования системы (1) существует вложеннный метод пятого порядка точности с оценщиком четвертого порядка.

Доказательство. Условия порядка для параметров четырехэтапного метода (2)-(4) пятого порядка с использованием упрощающих предположений:

ч

]Рав9у = сзд, в = 1,2, сц = 0, а1ЯЯ = 0, д = 1,...,4,

¿=1

образуют систему 31 нелинейного алгебраического уравнения с 25 неизвестными Ьу, Су, а^а вида

г' = 0'1'---'4' М = {(1,2), (2,1)},

3=1

4 • 9 1

Е Х^ а*91ст>] = 2Ц + зч' ^ = 2>

3=2 7 = 1 1 ^ 1

4 3

* 9 3 2

9=2 ¿ = 1 *=1 ^ ^ ;

4 /9 V 1

ЕМЕа'«сИ =20' П)

д=2 \, = 1 /

4 9 х

ЕЬ'<?Еа*«с^' = 20'

9=2

4 9 3 ^

У^ Е а1>43с-"3 Е а»7<с4* = ^п ' 9=2 ¿=2 «=1

4 9 .7 ^

XIЕXI = бо>

9=2 7 = 2

4 9 7« ^

Е ^ ЕЕ Е аз*рс*>р ~ ^20'

9=3 7 = 2 *=2 р=1

^ ^ — с»«1 г* = 1) ■ • •, 4.

7=1

Параметры пятиэтапного оценщика , с4 ^, четвертого порядка должны удовлетворять 18 равенствам

5 1

= ^О'1'2'3' М = {(1,2), (2,1)},

9= 1 5 3

Е = 2(1+3)' г ~

5 <7 х

3=2 ¿=1

5 9 3 ^

Е Е Е = 94'

д=2 7=2 ¿=1

5

Е =: с,5<

7 = 1

Для доказательства теоремы достаточно найти любое частное решение системы а; еб* раических уравнений (7), (8). В [3] достаточно подробно исследована система уравнениЦ (7) и выписано ее решение. Воспользуемся им. Это позволит свести исходную задачу к поиску решения системы (8). Таким образом, все параметры метода, кроме (1ад) с^, аг*> определены.

В результате исследования на непротиворечивость системы (8) (с учетом выполнения равенств (7)) ограничения на параметры у, с$о,принимают вид

5 ,

£ а1дс\д = —, г = 0,1,2,3, (5, г;) = {(1, 2), (2,1)},

9=1

4 5 1

ЕС2 9 X! = 12'

9=2 7=9+1

5

Е ¿1дСЧд1 = О,

9=2

4 5 х

Е X) =2' ^

9=2 7=^+1

5 • 1

= — г = 0,1,

7 = 1

5 1

Еа257Сг7 = дС?51 ¿21 = 0,

¿=1

достаточно простой для того, чтобы представить решение системы (7), (9), образующее шестипараметрическое семейство. Не выписывая все семейство, приведем в табл. 3 одно принадлежащее ему частное решение. Его существование и доказывает утверждение теоремы.

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

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

а 157 = Ьгу, а257 = Ьу, ; = 1,..., 4, (10)

<2255 = С21, С25 = 1 + С21, С15 = 1.

Таблица 3. РС5(4)5 метода (2)-(4)

Си

0 4 ч/б 15 15 4 ^/б 15 15

1 у/6 2 8 9 9л/б 32 128 7 7л</6 32 128

10 ^ 20 4977 4419^/6 9400 18800 2213 , 9809\/б 9400 112800 61 , 4469\/5 940 22560

1 63 , 36\/ё 47 47 35 1б7л/б 94 564 185 265л/6 94 564 0

Ьц 82 , 77ч/б 285 1140 297 351-\/б 1337 764 2432 , 64^6 2415 345 18184 250401 , 51676ч/б 250401 0

дц 41 13\/б 570 285 0 32 _ц64\/6 345 345 964 1311 32\/б 437 10 15

С2г' «2 г;

2 •/& 15 30 • 2 ¿/б 15 30

5 10 1 ^ 10 40 3 10 3^6 10

2 1 ^ 5 10 1337 , 1947\/б 1250 5000 4551 1750" 1083%/б 1000 8448 ', 496\/б 4375 625

1 103 83\/б 38 76 2901 , 382 11721^6 5348 72 272 л/б 23 161 62874 , 49236\/б 83467 1 83467

1 1439 , 869\/б 500 "т" 750 4857 700 1149л/б 350 4428 , 5576\/б 875 "т" 2625 1 3 1 3

0 4 9 36 9 ' 36 1 9 0

С?2г 0 4 9 л/б " 36 4 , ч/б 9 36 0 1 9

Теорема 2. В рамках пятиэтапного структурного метода (2)-(4) численного интегрирования системы (1) существует вложенный метод пятого порядка точности с оценщиком третьего порядка, удовлетворяющий равенствам (10).

Доказательство. Считаем, что параметры bsj,csj,asji удовлетворяют ограничениям (7) и всем упрощающим предположениям в [3], обеспечивающим пятый порядок точности четырехэтапного метода, и, как следует из [3], определены единственным образом. Параметры же пятиэтапного оценщика с^-, а^, помимо упомянутых выше равенств, должны еще удовлетворять системе

1

£

= —> г = 0,1, 2, (з, и) = {(1,2), (2,1)},

¿ + 1 9-1

-V] - б,

9=2 ]=1

'У ^ д,$д У ] азд]С

1

Е

= ¿-^ус'^1» * = О»1"» (И)

■

«15> = Ь2], а2ы = ] = !,..., 4,

а255 =С21, С25 = 1 + С21, С15 = 1.

Решения системы (11) образуют трехпараметрическое семейство

d25 = <, d21 = -t, t^O, d22 = ^ - ~ \/б + t - ^ ,

1131 591 /1620 г- 10354 /351 198 ,-Д

= + (12)

3034 2051 г- /160 190 >-\ /22468 9252 /Л 14 = ~4393 + 4393 Р V19I ~ 191 J + " V^393 4393 J ' ,136 4 г- /16 144 А 8 г

Формально полученный метод использует пять этапов по каждой компоненте, но фактически в случае удовлетворительного соотношения между оценкой /«погрешности:» и максимально допустимой заданной ее величиной последние вычисленные на текущем шаге коэффициенты kss(h), применяемые в оценщике, являются первыми на следую-

1 1 14 л/б

щем шаге. Полагая в (12) п = -, р = t = — — + ——, получим расчетную схему

3 3 69. 23

PC5(3)5F (табл. 4), которая по своим характеристикам превосходит известный [6] метод Дормана-Принса (DOPRI5(4)7F).

Таблица 4. PC5(3)5F метода (2)-(4)

си auj

0 15 15 JL 15 15

2 8 9 9\/б 32 128 7 32 7л/б 128

10 ~ 20 4977 4419л/б 9400 18800 2213 , 9400 9809V6 112800 61 | 4469V6 940 22560

1 0 4 9 ¿6 36 9 36 1 9

Ьц ' 82 , 77\/б 285 "г 1140 297 1337" 351V6 764 2432 , 64V6 2415 ' 345 ■18184 250401 ■ 51676л/б 250401 0

du 1 3 . 2103 1337" 117V6 1337 296 | 316V5 483 483 5682 4393 7469л/6 13179 1 3

C2i a2ij

2 yß[ 15 30 2 15 30

2 ¿6 5 10 1 & 10 40 3 3V6 10 10

5~ 10 1337 , 1947^6 1250"г 5000 4551 1083>/е5 1750 1000 8448 | 496\/б 4375 625

1 103 83л/б 38 76 2901 , 11721\/б 382 "г 5348 72 272 т/б 23 161 62874 , 49236\/б 83467 1 83467

17 V6 15 30 82 , 77\/б 285 1140 297 351\/б 1337 764 2432 , 64л/6 2415 345 18184 , 51676\/6 250401 1 250401 2 ^/6 15 30

¿>21 0 4 л/б 9 36 4 , Уб 9 36 1 9 0

d2i 1 i/6 46 23 55 , 5л/6 138 138 55 1 5\/б 138 "Г" 138 14 5\/б 69 69 46 23

Приведение РС5(3)5Г к виду (б) для прямого интегрирования Дифференциального уравнения второго порядка (5) в силу отсутствия алгоритмической симметрии метода (2)-(4) допускает два представления.

Первая - РС1Н5(3)5Г (табл. 5) вида (б) получена при замене переменных у\ = у', У2 = У ъ уравнении (5) и применении к системе первого порядка метода РС5(3)5Е. Вторая - РС2Н5(3)5Р (табл. 6) - при замене уг = у, у2 — у1.

Таблица 5. РС1Н5(3)5Р метода (б)

.С,- Ац

0 0

±_ у/5 15 15 1 & 2 8 11 4 У5 225 225 209 19У? 2560 640 231 21\/б 2560 640

7 , у/5 •10+201 . А41 82 , 77\/б 285 1140 А42 927 1881-А 2674 5348 Ащ 1552 , 17Вт/б 2415 "г 805 6986 , 16412л/в 83467 1 250401

Во} Ак 82 . 77\/6 285 1140 3403 , 96л/6 6900 575 927 1881 у/5 2674 5248 427863 399867т/<Г 615020 615020 1552 , 17в\/б 2415 805 12896 . 1248 л/б" 12075 4025 6986 , 1в412л/б 83467 ' 250401 4664 | 732У? 13179 1 4393 0 4 ■ Зу/5 345 460

Ви £>и 82 1 77у5 285 "г" 1140 1 3 297 351-\/б 1337 764 2103 117\/б 1337 1337 2432 , 64^6 2415 ■* 345 296 ■ 316у/5" 483 ' 483 18184 | 51676Л/6* 250401 1 250401 5682 74в9\/б 4393 13179 0 1 3

л _ 4902303 | 1320067У% л _ 14304369 5743641У% д _ 840941 ■ 340324У5 ^41 — 9400000 9400000 ' 42 ~ 13160000 13160000 ' 43 ~~ 1028125 1028125 "

Таблица 6. РС2Н5(3)5Р метода (6)

Сг А*]

2 у/5 15 30 2 у/5 5 10 2 , ^ 5 т 10 1 17 у/5 15 30 0 и ¿/в 100 25 13 7-у/б 81 I 17-4/5" 250 250 500 ' 250 1 . 5у/5 \/5 з Зу/5 4"1" 16 8 4 16 о _2__¿5 54 ' 108 135 180 135 270

£>о; 0 0 0 749 у/5 775 , 403ч/б 12373 , 6137у/5 1 п 752 30 10152 1 5076 20304 1 20304 27 и

Ви Ои о ±+У£ I 0 " 9 36 9 ~ 36 9 и 1 у/5 55 ,5у/б 55 .Бт/б 14 5 ч/б 1 , у/5 46 23 138 "Г 138 138 "" 138 69 69 46 23

Обе схемы интегрирования предназначены для прямого интегрирования уравнения (5) и могут быть классифицированы [7] как методы типа Нюстрёма. При одинаковых точностных и затратных характеристиках (в построенных расчетных схемах) между ними есть и отличие. РС2Н5(3)5Г содержит два момента, расширяющих традиционные постановочные требования к методам типа Нюстрёма. Для параметров метода (6) равенства С\ — 0 и Ст = 1 являются необязательными.

3. Численное тестирование. Численное сравнение проводилось на решении задачи Коши

у" = -у + a cos bx, у(х0) = уо, у'(хо) = y'Q. При 6 ф ±1 общее решение имеет вид

. . а , у{х) = ——jp cos Cl sm ж + c2 cos xi

' / \ ab . ,

у (x) ~ — -—sm bx + cj cos x — C2 sin x.

В серии расчетов изменялась допустимая погрешность toi (контрольный член [9]) от 10~2 до КГ11'8, причем каждое следующее значение было в л/ÏO раз меньше предыдущего. Интегрирование производилось на интервале х 6 [0, 5, 5-тг] при параметрах системы а = 5 и b = 1/2. В качестве оппонента полученному методу - PC5(3)5F -был выбран [7] метод Дормана-Принса - DOPRI5(4)7F как лучший из класса вложенных с подобными точностными характеристиками. Результаты тестирования представлены на рис. 1-4 и в табл. 7. Так, на рис. 1 и 2 приведены зависимости точности (<5,- = maxrg[o, 5,5тг] \yi(x)—yi\ - покомпонентная максимальная глобальная погрешность) и трудоемкости (/с - количество вычислений правой части) от величины ограничения «погрешности» (toi - предельно допустимое значение контрольного члена).

-lg(j№)

одм = max |у,(ж;) - у,|; xj € [0;5;5/?]);

от допустимого значения контрольного члена toi.

Для как характерно то, что по оси абсцисс (со знаком «минус») отложены десятичные логарифмы значений допустимой погрешности, а по оси ординат - логарифмы максимальной глобальной погрешности (рис. 1) и логарифмы количества вычислений правой части (рис. 2) соответственно. Эволюцию глобальной погрешности обеих компонент решения иллюстрирует рис. 3, а принятого шага интегрирования - рис. 4.

-Igfc

Рис. 3. Покомпонентное развитие глобальной погрешности

fi = \Vi(х)- pil методов DOPRI5(4)7F (itoi = HT**,/с = 1070), PC5(3)5F (toi = 1<Г"*,Л = 1090).

Рис. 4■ Длина принятого шага для методов DOPRI5(4)7F (;toi = 1(Г5*,/С = 1070), PC5(3)5F (toi = 10~**,fc = 1090).

Таблица 7. Зависимость максимальной глобальной погрешности и трудоемкости от величины допустимого значения контрольного члена

-lg (toi) DOPRI5(4)7F PC5(3)5F

-Wll /с -Ы6\\ /с

3,00 2,510 350 0,705 154

3,50 3,018 434 1,253 181

4,00 3,494 542 1,203 208

4,50 3,926 686 2,379 253

5,00 4,264 854 2,397 289

5,50 4,449 1070 2,869 334

6,00 4,522 1346 3,251 388

6,50 4,547 1682 3,522 451

7,00 4,555 2114 4,003 523

7,50 . 7,017 2 726 4,227 604

8,00 7,503 3434 4,411 694

8,50 7,957 4310 4,474 802

9,00 8,338 5414 4,496 919

9,50 8,591 6806 5,948 1090

10,00 8,716 8570 6,368 1261

10,50 8,764 10778 6,726 1450

11,00 10,525 13634, 7,054 1675

г

Анализ результатов численного тестирования (рис. 1-4, табл. 7) показал, что: - 1) метод DOPRI5(4)7F требует меньшего (чем PC5(3)5F) допуска на погрешность, чтобы получаемые в действительности погрешности были сравнимы;

2) алгоритм формирования шага интегрирования PC5(3)5F чутко реагирует на изменения глобальной погрешности;

3) при равных вычислительных затратах глобальная погрешность PC5(3)5F на порядок меньше, чем DOPRI5(4)7F;

4) при одном и том же порядке глобальной погрешности метод DOPRI5(4)7F использует на треть больше обращений к процедуре вычисления правых частей.

Следовательно, и численное сравнение подтверждает эффективность полученного в статье метода.

В заключение отметим, что рассмотренный структурный подход указывает способ конструирования методов интегрирования систем обыкновенных дифференциальных уравнений первого порядка специального вида, которые после их приведения к форме, экономичной и алгоритмически простой для прямого интегрирования дифференциального уравнения у"(х) = f(x,y(x)), попадают в класс методов типа Нюстрёма (6), делая ограничение С\ = 0 необязательным.

Summary

Olemskoy I. V. The fifth order embedded method.

The obvious one-step-by-step method of numerical integration of the system of the ordinary differential equations of a special kind is considered. The statements of correspondence of accuracy order and number of stages are proved. Within the framework of the suggested method economical numerical schemes up to the fifth order with the mechanism of an automatic choice of a step of integration are constructed.

Литература

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

2. Олемской И. В. Экономичная расчетная схема четвертого порядка точности численного интегрирования систем специального вида //Труды XXX науч. конференции «Процессы управления и устойчивость». СПб., 1999. С. 134-143.

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

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

5Г Олемской И. В. Методы типа Рунге-Кутты интегрирования систем и дифференциальных уравнений второго порядка специального вида // Вычисл. технологии. 2004. Т. 9, № 2. С. 67-81. *

6. Современные численные методы решения обыкновенных дифференциальных уравнений /Под ред. Дж. Холла, Дж. Уатта; Пер. с англ.; Под ред. А. Д. Горбунова. М., 1979. 312 с.

7. Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи/Пер. с англ.; Под ред. С. С. Филлипова. М., 1990. 512с.

8. Bormand 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.

9. Арушанян О. В., Залеткин С. Ф. Численное решение обыкновенных дифференциальных уравнений на Фортране. М., 1990. 336 с.

Статья поступила в редакцию 10 мая 2004 г.