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

УДК 519.622 Е.А. Новиков, Ю.А. Шитов

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

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

Введение. Почти все известные методы типа Розенброка [1] и их модификации предназначены для численного интегрирования задачи Коши для автономных систем обыкновенных дифференциальных уравнений вида

х = £(х), х(г0) = х0, г0 < г < гк, (1)

где х и д - вещественные (М +1)-мерные вектор-функции; I - независимая переменная. Такой подход оправдан тем, что неавтономную систему

У= I(г,У)> У(го) = Уo, го <г < гк (2)

введением дополнительной переменной можно привести к виду (1), если положить х = (Ур---, УМ, г )Т,

д = (/1>•••,/м,1) . Использование автономной формы записи существенно упрощает построение схем

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

У.+1 = У, + ЁРк, , О, = Е-аН 0/{(-'У-},

(=1 (У

I-1 I-1 I-1 0/ (г У )

О,к,= Н/ (г, + (Ё в, )Н, У, + Ё в, к,) + (а + Ё Р, )Н 0 I ’У , (3)

,=1 ,=1 ,=1 01

в то время как методы решения неавтономных задач следующий:

т ,-1

Уп+1 = У, + Ё Рiki ’ °пк, = Н/ (гп + С,Н’ У, + Ё вк, ) , (4)

1=1 ,=1

где а, р{, , в,, 1 <, < т, 1 < , <,-1, - коэффициенты. Как следует из [2], наличие частной произво-

дной д/ (г,, у, )/ дг в расчетной формуле во многих случаях ухудшает устойчивость схемы и может приводить к потере О -устойчивости.

Здесь построен L-устойчивый (т, к)-метод [3] третьего порядка точности для численного интегрирования неавтономных задач вида (2). Получены оценка ошибки и неравенство для контроля точности вычислений и автоматического выбора величины шага интегрирования.

Численная схема. Для решения (2) применим (3, 3)-метод вида

У,+1 = У, + Р1к1 + Р2к2 + Р2к3 ,

°ик1 = Н/ (г, , У, ) , ОА = Н/ (г, + С2Н У, +в21к1) + а21к1 , (5)

О,к3 = Н/ (г, + С3Н, У, + в31к1 + Р31к2 ) ,

где Б, = Е - аН/ '(г,, у,); Н - шаг интегрирования; Е - единичная матрица;

/'(г,,У,) = д/(г,,у,)/ду - матрица Якоби задачи (2), а, р^, рг, р3, р1Х, Р31, Р32, а21 - коэффици-

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

к1 = Н/ + аН2 /у/ + а2 Н3 /У / + 0(Н4),

к2 = (1 + 0,21 )Н/ + Н2[(а + в21 + 2аа1Х) /у/ + с 2 / ] + Н^[(а2 + 2ав21

+3а2а21)/У/ + 0.5в221/уу/2 + св/у/ + ас2/у/ + 0.5с1 / + 0(НА),

к3 = Н/ + Н2[(а + Р31 + Р32 +а21в32) /у/ + С3/г] + + 2аРз1 (6)

+2аР32 + в2р32 + 3аа2р32) /У / + 05(в31 + Р32 + а21Р32)1 /уу/2 +С3(Р31 + Р32 + а21Р32 ) /гу/ + (аС3 + С2р) /у/ + 0-5с3 /гг] + 0(Н^ .

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

разложения к1, к2 и к3 в ряды Тейлора в окрестности точки г, до членов с Н3 включительно

к1 = Н/ + аН2 /у/ + а2 Н3 /у / + 0(Н4), к2 = (1 + 0,21 )Н/ + Н2[(а + Р21 + 2а&21) /у / + с 2 / ] + Н^[(а2 + 2аР21 +3а2а21)/У/ + 0.5Р21/уу/2 + с2Р21/у/ + ас2/у/( + 0.5с\/ + 0(НА),

к3 = Н/ + Н2[(а + Р31 + Р32 +а21Р32) /у/ + с3/г] + + 2аР31 (7)

+2аР32 + Р21Р32 + 3аа21Р32) /У / + 05(Р31 + Р32 + а21Р32)1 /уу/2 +с3(Р31 + Р32 + а21Р32 ) /у/ + (ас3 + с2Р32 ) /у/г + 0-5с3 Л] + 0(Н4) .

Подставляя (7) в первую формулу (5) и при условии у, = у( гп), сравнивая ряды Тейлора для точного у(г,+1) и приближенного уп+1 решений до членов с Н3 включительно, получим соотношения на коэффициенты схемы (5), обеспечивающие третий порядок точности, то есть

Р1 + (1 + СХ21) Р2 + Р3 = 1, ар1 + (а + Р21 + 2а&21) Р2 + (а + Р31 + Р32 + ОС21Р32) Р3 = 0.5,

с2 р2 + с3 р3 = 0'5 , с2 р2 + с3р3 = 1/3 , Р21Р2 + (Р31 + Р32 + а21Р32) Р3 = 1/3 ,

а2 Р1 + (а2 + 2аР21 + 3а 2а>21) Р2 + (а2 + 2а ^1 + 2аР32 + Р21Р32 + 3а&21р32) Р3 = 1/6,

с2 Р21Р 2 + с3(Р31 + Р32 +а21Р32) Р3 = 1/3 , ас2 Р2 + (ас3 + с2Р32) Р3 = 1/6 .

Пользуясь произволом при выборе коэффициентов, положим с2 = Р21 и

с3 = Р31 + Р32 + а21Р32. Тогда данную систему можно записать в виде

1) р1 +(1 + ^21) р2 + р3 =1,

2) ар1 + (а + Р21 + 2а&21) Р2 + (а + Р31 + Р32 + ОС21Р32) Р3 = 0.5,

3) Р21 р2 + (Р31 + Р32 + а21Р32)р3 = 0-5 ,

4) а2р1 + (а2 + 2аР21 + 3а2а21)р2 + (а2 + 2аР32 +Р21Р32 + 3аа21Р32) = 1/6, (8)

5) Р21 р2 + (Р31 + Р32 +а21Р32) р3 = 1/3 ,

6) аР21 р2 + \-а(Р31 + Р32 +а21Р32) + Р21Р32]Р3 = 1/6 .

Исследование совместности (8) приводит к следующему набору коэффициентов: р1 = 1 -(1 + а21)р2 - р3, р2 = (1 -3а)/(3Р21 -6аР21), р3 = 1/(6с3 - 12ас3),

(0,21 = 3Р2Х(2а -1)/(1 -3а), Р31 = 4(3а -3а + 1)с3 -Р32, Р32 = с3(1 -3а)1 - 2а)/Р21, (9)

с2 = Р21 , с3 = Р21 + (1 - 2а)(2 - 3Р21 ) ,

при которых численная формула (5) имеет третий порядок точности. Здесь коэффициенты а и Р21 произвольные.

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

У' = Лу, у(0) = У0, г > 0 (10)

с комплексным Л, Яе(А) < 0 . Применяя (5) к (10), получим уп+1 = Q(х)уп, х = НЛ, где функция устойчивости Q(x) задается формулой

Q(x) = {1 + х(р1 + р2 + р3 +а21 р2 -3а) + х2[3а2 -2а(р1 + р2 + р3) + Р21 р2 + (Р31 + Р32)р3

-аа21 Р2 +@21032 Р3] + х [-а + а (р^ + Р2 + Р3) - а(^2\Р2 + (Р31 + Р'32) Р3 + $21032 Р3]}/(1 - ах) . Учитывая значения коэффициентов (9), ее можно переписать в виде

Q(х) = [1 - х(3а -1) + 0.5 х2 (6а2 - 6а +1) - х3 (а3 - 3а2 + 1.5а -1/ 6~)]/(\ - ах)3. (11)

Необходимым условием L-устойчивости численной формулы (5) является выполнение соотношения I Q(х) 0 при х ^ -». Из (11) следует, что данное требование будет выполнено, если коэффициент

а есть корень уравнения

а3-3а2 + 1.5а -1/6 = 0. (12)

Известно [4], что при 1/3 < а < 1.0685790230270 схема (5) будет обладать свойством А-устойчивости. Поэтому среди корней данного уравнения (12) выберем а = 0.43586652150902, при котором (5) является L-устойчивой.

Контроль точности вычислений. Локальная ошибка 8, схемы (5) с коэффициентами (9) имеет вид

8, = гН /у/+0(НА),

где У1 = 1/24 - а3 - 3а2 (Р21 + а&2\)Р2 -3а(ас3 + а&2\Р32 + Р21Р32)Р3. (13)

Выберем коэффициенты Ъ1, I < I < 4, из условия выполнения соотношения

Ь1к1 + Ъ2к2 + Ъ3к3 + Ъ4Н/(гп + с2Н Уп +Р21Ю = 72Н/У/ + 73Н'/у/, + 0(Н4) . (14)

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

£п = (71 / 72 )[Ъ1к1 + Ъ2к2 + Ъ3к3 + Ъ4Н/ (гп + с2Н У, + Р21к1 )] . (15)

Отметим, что задание е, по формуле (15) к существенным дополнительным вычислительным затратам не приводит, потому что в (15) применяются ранее вычисленные стадии. С использованием разложений в ряды Тейлора нетрудно убедиться, что требование (14) будет выполнено, если Ъ + (1 + СХ21 )Ъ2 + Ъъ + Ъ4 = 0,

аЪ1 + (а + Р21 + 2аа21 ')Ъ2 + (а + Р31 + Р32 + а21Р32 )Ъ3 + Р21Ъ4 = 0 ,

с2Ъ2 + с3Ъ3 + с2Ъ4 = 0 , Р21Ъ2 + (Р31 + Р32 + а21Р32 ) Ъ3 + Р21 = 0 , (16)

с2Р21Ъ2 + с3 (Р31 + Р32 + а21Р32 ')Ъ3 + с2Р21Ъ4 = 0 , с2Ъ2 + с3Ъ3 + с2Ъ4 = 0 .

Учитывая, что в (9) коэффициент Р21 является произвольным, положим с2 = с3 = Р21. В этом случае (16) можно переписать в виде

Ъ + (1 + (Х21 )Ъ2 + Ъ3 = -Ъ4 , Ъ2 + Ъ3 = -Ъ4 , Ъ2 = Ъ4 / (Х21 .

Пусть Ъ4 = (4а - 2)/(а -1). Тогда коэффициенты (5) имеют вид р1 = 5/4, р2 = (1 -3а)/(2-4а), р3 = 1/(4 -8а), с2 = с3 = Р21 = 2/3, аа21 = (4а - 2) /(1 - 3а), Р31 = 2а2 - 3а + 5/3, Р32 = 6а - 5а +1, (17)

а коэффициенты Ъ, 1 < I < 4,

Ъ1 = (2 -4а)/(а -1), Ъ2 = (1 -3а)/(а -1), Ъ3 =-1, ЪА = (4а - 2)/(а -1). (18)

При этом имеют место соотношения

Y = (48a4 - 96a3 + 60a2 - 14a +1) /(24 - 48a),

Y2 = (24a4 — 48a3 + 38a2 — 14a + 2) /(3a - 3), (19)

Y3 = (-12a3 + 14a2 - 8a + 2)/(3a - 3).

Теперь для контроля точности (5), (17) можно применить неравенство || £п ||< £, где || ■ || - некоторая

норма в Rn , £ - требуемая точность интегрирования, а £ п вычисляется по формуле (15) с коэффициен-

тами (17) и (18).

Отметим одну важную особенность использования неравенства || £п ||< £ для контроля точности вычислений и при выборе шага интегрирования. В случае решения задачи (10) имеет место соотношение

£ п = (Y2 X + Y4 xA) У n /(1 - ax)3, x = HA, (20)

где у1 вычисляется по второй формуле (19), а уА имеет вид

Ya = (-12a4 + 14a3 -4a2)/(3a -3).

Численная схема (5), (17) является L-устойчивой и, следовательно, для нее выполняется | Q(x) ^ 0 при x ^ -^ (см. (11) и (12)). Так как для точности решения y(tn+1) = exp(x)y)tn) задачи (10) выполняется аналогичное свойство exp(x) ^ 0 при x ^ -те, то естественным является требование стремления к нулю ошибки при x ^ -те. Из вида (20) следует, что для £п это требование не выполняется, что может приводить к неоправданному уменьшению шага и, как следствие, к понижению эффективности алгоритма интегрирования. Это связанно с тем, что в [5] при получении £п все рассуждения проведены для главного члена ошибки.

С целью исправления асимптотического поведения ошибки (20) вместо построенной выше оценки £п будем рассматривать £п (jn), определенную по формуле

£ (j ) = D1 j£ ,1 < j< 3, (21)

а вместо контроля || £п ||< £ будем проверять неравенство || £п (jn) ||< £ . Учитывая, что

D- = E + ahfy + O(H2), нетрудно видеть, что в смысле главного члена, то есть в смысле первого члена при разложении ошибок в ряды Тейлора по степеням H, £п и £п (jn) совпадают при всех jn, 1 < j п < 3. Отметим, что применение

(21) вместо £п к существенному увеличению вычислительных затрат не приводит. При x ^ 0 оценка £ п (1) = £ п правильно отражает поведение ошибки и нет смысла проверять (21) при других значениях j п. При резком увеличении шага поведение £ может оказаться неудовлетворительным, что проявляется в неоправданном уменьшении шага и повторных вычислениях решения (возвратов). Проверка (21) при j п Ф 1 осуществляется редко. Однако, как показывают расчеты, использование (21) вместо || £п И<£ сокращает вычислительные затраты на 10-20% .

Литература

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

2. Современные численные методы решения обыкновенных дифференциальных уравнений / под ред. Дж. Холла и Дж. Уатта. - М.: Мир, 1979. - 312 с.

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

4. Демидов, Г.В. Исследование некоторых аппроксимаций в связи с А-устойчивостью полуявных методов / Г.В. Демидов, Л.А. Юматова // Числен. методы механ. сплош. среды. - 1977. - Т. 8. - №3. - С. 68-79.

5. Демидов, Г.В. Оценка ошибки одношаговых методов интегрирования обыкновенных дифференциальных уравнений / Г.В. Демидов, Е.А. Новиков // Числен. методы механ. сплош. среды. - 1985. - Т. 16. -№ 1. - С. 27-39.