CC BY
48
Вычислительные технологии
Том 3, № 3, 1998
(m, 3)-МЕТОД ТРЕТЬЕГО ПОРЯДКА ДЛЯ ЖЕСТКИХ НЕАВТОНОМНЫХ СИСТЕМ ОДУ
Е. А. Новиков Институт вычислительного моделирования СО РАН
Красноярск, Россия
М. И. ГОЛУШКО Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: [email protected]
The L-stable (m, 3)-method of the third order for solving the Cauchy problem for stiff non-autonomous systems of ordinary differential equations has been constructed. The results of the calculations, which confirm efficiency of the algorithm of integration, have been adduced.
1. Введение
Почти все известные методы типа Розенброка [1] и их модификации предназначены для численного интегрирования задачи Коши для автономных систем вида
х' = д(х), х(Ьо) = Хо, ¿о < г < гк, (1)
где х и д — гладкие вещественные (И + 1)-мерные вектор-функции, г — независимая переменная. Такой подход оправдан тем, что неавтономную систему
у' = /(г,у), у(1о) = уо, ¿0 < г < ¿к (2)
введением дополнительной переменной можно привести к виду (1), если положить х =
(у!,...,ум ,г)Т, д =(/!,...,/м, 1)Т.
Использование автономной формы записи существенно упрощает построение схем заданного порядка точности. Однако методы типа Розенброка, предназначенные для решения автономных задач, применительно к (2) имеют вид
уп+1 = уп + Ргкг, = Е - акд/,
г=1 у
i— 1 i— 1 i— 1 rj г/, \
Dnki = hf (tn + Y,вцh,Уп + Y,виkj) + (a + ^вц)h f (^Уп), (3)
j=i j=i j=i
© Е.А. Новиков, М.И. Голушко, 1998.
в то время как методы решения неавтономных задач —
m i-i
Уп+1 = Уп + Piki, Dnki = hf (tn + Cih, yn + ^ 0jkj), (4)
i=l j=i
где a, pi, ci, 0ij, 1 < i < m, 1 < j < i — 1 — параметры. Как следует из [2], наличие частной производной df (tn,yn)/dt в расчетной формуле во многих случаях ухудшает устойчивость схемы и, в частности, может приводить к потере D-устойчивости.
В настоящей работе построен L-устойчивый (m, к)-метод [3] третьего порядка точности для численного интегрирования неавтономных задач. Предложен способ контроля точности вычислений и приведены результаты тестовых расчетов.
2. Численная схема
Для решения (2) применим (3, 3)-метод вида
Уп+1 = уп + Pl ki + Р2к2 + Plh, Dnki = hf (tn,yn), Dnk2 = hf (tn + c2h, yn + 02iki) + «21 ki,
Dnki = hf (tn + Cih, yn + p3iki + 032 k2), (5)
где Dn = E — ahf'(tn, yn), h — шаг интегрирования, E — единичная матрица, f = df/dy — матрица Якоби задачи (2), a,p1,p2,p3,021,031,032,а21 — коэффициенты, обеспечивающие заданные свойства точности и устойчивости.
Ниже потребуется разложение точного решения y(tn+1) в окрестности точки tn до членов с h4 включительно, которое имеет вид
y(tn+i) = y(tn) + hf + 0.5h2(fy f + f) + j(f2yf + fyy f2 + 2fty f+
h4 3 2 2 3
+fy ft + ftt) + 24(fy f + fy fyyf + 3fyy fyf + fyyyf +
+ 3fyytf 2 + 3fyttf + 3 fytfyf + 3fytft + 3 fyy ftf + 2fy fytf +
+ fy2ft + fy ftt + fttt)+ O(h5). (6)
Здесь нижний индекс означает соответствующую частную производную, а элементарные дифференциалы вычислены в точке tn.
Для построения метода третьего порядка точности запишем разложения ki , k2 и k3 в ряды Тейлора в окрестности точки tn до членов с h3 включительно:
ki = hf + ah2fy f + a2h3fy2f + O(h4), k2 = (1 + «2i)hf + h2[(a + 02i + 2aa2i )fy f + f] +
+h3[(a2 + 2a02i + 3a2a2i)fy2f + 0.5^1 fyy f2 + c202ifty f + ac2fy ft + 0.5cf + O(h4), k3 = hf + h2[(a + P31 + в32 + ^2ip32)fy f + C3ft] + +h3[(a2 + 2a03i + 2ap32 + 021032 + 3a«2i032)f2f + 0.5(03i + 032 + «21032)2 fyy f 2+
+С3(в31 + 032 + «21032 )fty f + (ac3 + C2 032) fy ft + 0.5c3ftt] + O(h4). (7)
Подставляя (7) в первую формулу (5) и сравнивая ряды Тейлора для точного и приближенного решений до членов с включительно, получим соотношения на параметры схемы (5), обеспечивающие третий порядок точности, т.е.
1) Р1 + (1 + «21)Р2 + Рз =1,
2) ар1 + (а + ^21 + 2аа21)р2 + (а + вз1 + вз2 + «21^з2)рз = 0.5,
3) С2Р2 + сзрз = 0.5,
4) а2Р1 + (а2 + 2ав21 + 3а2а21)р2 + (а2 + 2авз1 + 2авз2 +
+в21вз2 + 3аа21вз2)Рз = 1/6,
5) в21Р2 + (вз1 + вз2 + «21вз2)2Рз = 1/3,
6) С2^21Р2 + Сз(вз1 + вз2 + «21^з2)Рз = 1/3,
7) ас2Р2 + (асз + С2вз2)Рз = 1/6,
8) с2Р2 + сзРз = 1/3. (8)
Пользуясь произволом при выборе коэффициентов, положим с2 = в21, сз = вз1 + вз2 + а21вз2. Тогда (8) можно записать в виде:
1) Р1 + (1 + «21)Р2 + Рз = 1,
2) ар1 + (а + ^21 + 2аа21)р2 + (а + вз1 + вз2 + «21^з2)рз = 0.5,
3) в21Р2 + (вз1 + вз2 + «21вз2)Рз = 0.5,
4) а2р1 + (а2 + 2ав21 + 3а2«21)р2 + (а2 + 2авз1 + 2авз2 +
+в21вз2 + 3аа21вз2)рз = 1/6,
5) в221Р2 + (вз1 + вз2 + «21вз2)2Рз = 1/3,
6) ав21Р2 + [а(вз1 + вз2 + «21вз2) + в21вз2]рз = 1/6. (9)
Исследуем совместность нелинейной системы алгебраических уравнений (9). Из первого, второго и четвертого уравнений (9) имеем
(в21вз2 + а«21^з2)рз = 1/6 - а - а2. Из третьего и шестого уравнений (9) получим
вз2Рз = (1/6 - 0.5а)/в21, «21Р2 = -1. С учетом (10) и (11) выразим а21 и р2:
«21 = 3^21 (2а - 1)/(1 - 3а), Р2 = (1 - 3а)/(3^21 - 6а^21).
Из третьего уравнения (9) имеем
(вз1 + вз2 + «21вз2)рз = 1/(6 - 12а).
Тогда из третьего, пятого и первого уравнений (9) получим
рз = 1/[(6 - 12а)(2 - 4а - 2&1 + 6а&1)], Р1 = 1 - (1 + «21 )Р2 - Рз.
21
(10)
(11)
(12)
(13)
С учетом (11), (13) и первого соотношения (14) запишем
в32 = (1 - 3а)(1 - 2а)(2 - 4а - 2021 + ба^г)/021,
вз1 = 4(3а2 - 3а + 1)(2 - 4а - 2021 + 6а02г) - вз2- (15)
В итоге получим следующий набор параметров схемы (5):
С2 = 021, «21 = 302г(2а - 1)/(1 - 3а), Р2 = (1 - 3а)/(3в21 - 6а021), сз = 021 + (1 - 2а)(2 - 3021),
Рз = 1/(6сз - 12асз), Р1 = 1 - (1 + а21)р2 - Рз, 032 = сз(1 - 3а)(1 - 2а)/021, 0з1 = 4(3а2 - 3а + 1)сз - 0з2, (16)
при которых формула (5) имеет третий порядок точности. Здесь коэффициенты а и 021 произвольные.
Исследуем устойчивость схемы (5) на линейном скалярном уравнении
у' = Ау (17)
с комплексным А, И,е(А) < 0. Применяя (5) к (17), получим
Уп+1 = Я(х)уи, х = АН, где функция устойчивости ((х) задается формулой
((х) = {1 + х(р1 + Р2 + Рз + 0-21Р2 - 3а) + х2[3а2 - 2а(р1 + Р2 + Рз) +
+021Р2 + (0з1 + 0з2)Рз - аа21Р2 + «210з2Рз] +
+хз[-аз + а2('Р1 + Р2 + Рз) - а(021Р2 + (0з1 + 0з2)'Рз) + 0210з2Рз]}/(1 - ах)3. Учитывая значения коэффициентов (16), ее можно переписать в виде
((х) = [1 - х(3а - 1) + 0.5х2(6а2 - 6а + 1) - хз(аз - 3а2 + 1.5а - 1 /6)]/(1 - ах)з. (18)
Необходимым условием ¿-устойчивости формулы (5) является выполнение соотношения
К(х)| ^ 0 при х ^ -то.
Из (18) следует, что данное требование будет выполнено, если параметр а есть корень уравнения
аз - 3а2 + 1.5а - 1/6 = 0. (19)
Известно [4], что при 1/3 < а < 1.0685790230270 схема (5) будет обладать свойством -устойчивости. Поэтому среди корней уравнения (19) выберем а = 0.43586652150902, при котором (5) является ¿-устойчивой.
3. Алгоритм интегрирования
Локальная ошибка 8п схемы (5) с параметрами (16) имеет вид
6п = Ъ^Щ + 0(Н4), (20)
где
71 = 1/24 - аз - 3а2(021 + аа21 )Р2 - 3а(асз + аа210з2 + 0210з2)Рз. (21)
Выберем коэффициенты Ьг, 1 < г < 4 из условия выполнения соотношения
М1 + Ь2&2 + Ьзкз + Ь4й/(1п + С2й, уп + ^Л) = 72йз/2/ + 7зй/ / + 0(й4). (22)
Тогда согласно [5] для контроля точности вычислений и при выборе шага интегрирования можно использовать оценку ошибки еп, вычисленную по формуле
Т1
£п =-[Ь1^1 + Ь2&2 + Ьзкз + Ь4 й/(Ьп + С2й,Уп + ^Л)]. (23)
72
Отметим, что определение еп по (23) к существенным дополнительным вычислительным затратам не приводит, потому что в этой формуле применяются ранее вычисленные стадии.
С использованием разложений в ряды Тейлора нетрудно убедиться, что требование (22) будет выполнено, если
1) Ь1 + (1 + «21)Ь2 + Ьз + Ь4 = 0,
2) аЬ1 + (а + ^21 + 2аа21)Ь2 + (а + вз1 + вз2 + «21вз2)Ьз + ,021^ = 0,
3) С2Ь2 + СзЬз + С2Ь4 = 0,
4) в21 Ь2 + (вз1 + вз2 + «21вз2)2Ьз + = 0,
5) С2^21Ь2 + Сз(вз1 + вз2 + «21вз2)Ьз + С2^21Ь4 = 0,
6) с2Ь2 + сзЬз + с2Ь4 = 0. (24)
Учитывая, что в (16) коэффициент в21 является произвольным, положим с2 = сз = в21. В этом случае (24) можно переписать в виде
Ь1 + (1 + «21 )Ь2 + Ьз = -Ь4,
Ь2 + Ьз = -Ь4,
Ь2 = Ь4/а21. (25)
Пусть Ь4 = (4а - 2)/(а - 1). Тогда параметры (16) имеют вид
С2 = Сз = в21 = 2/3, «21 = (4а - 2)/(1 - 3а),
р1 = 5/4, р2 = (1 - 3а)/(2 - 4а), рз = 1/(4 - 8а),
вз1 = 2а2 - 3а + 5/3, вз2 = 6а2 - 5а + 1, (26)
а коэффициенты Ьг, 1 < г < 4, —
Ь1 = (2 - 4а)/(а - 1), Ь2 = (1 - 3а)/(а - 1),
Ьз = -1, Ь4 = (4а - 2)/(а - 1), (27)
при этом имеют место соотношения
71 = (48а4 - 96аз + 60а2 - 14а + 1)/(24 - 48а),
72 = (24а4 - 48аз + 38а2 - 14а + 2)/(3а - 3),
7з = (-12аз + 14а2 - 8а + 2)/(3а - 3). (28)
Теперь для контроля точности (5), (26) можно применять неравенство
1Ы|< е, (29)
где ||-1| — некоторая норма в RN, е — требуемая точность интегрирования, а еп вычисляется по формуле (23) с коэффициентами (27), (28).
Отметим одну особенность использования (29) для контроля точности вычислений и при выборе шага интегрирования. В случае решения задачи (17) имеет место соотношение
Yi + 14Х4 , ,
еп =---Тл-^Т Уп, x = Xh, (30)
Y2 (1 - ax)3
где вычисляется по второй формуле (28), а y4 — по следующей:
Y4 = (-12а4 + 14а3 - 4а2)/(3а - 3).
Численная схема (5), (26) является L-устойчивой, и, следовательно, для нее выполняется IQ(x)l ^ 0 при x ^ -то (см. (18) и (19)). Так как для точного решения y(tn+\) = exp(x)y(tn) задачи (17) выполняется аналогичное свойство exp(x) ^ 0 при x ^ -то, то естественным является требование стремления к нулю ошибки при x ^ -то. Из (30) видно, что для еп это требование не выполняется, что может приводить к неоправданному уменьшению шага и, как следствие, к понижению эффективности алгоритма интегрирования. Это связано с тем, что в [5] при получении еп все рассуждения проведены для главного члена ошибки.
С целью исправления асимптотического поведения ошибки вместо построенной выше оценки еп будем рассматривать величину еп(Зп), определенную по формуле
еп(]п) = Dl~jnеп, 1 < Зп < 3, (31)
а вместо контроля (29) станем проверять следующее неравенство:
||еп(Зп)11< е, 1 < Зп < 3. (32)
Учитывая, что
D- 1 = E + ahfy + O(h2), (33)
нетрудно видеть, что в смысле главного члена, т. е. в смысле первого члена при разложении ошибок в ряды Тейлора по степеням h, еп и еп(зп) совпадают при всех Зп, 1 < Зп < 3.
Отметим, что использование (32) вместо (29) к существенному увеличению вычислительных затрат не приводит. При x ^ 0 оценка еп(1) = еп правильно отражает поведение ошибки, и нет смысла проверять (32) при других значениях зп. При резком увеличении шага поведение еп может оказаться неудовлетворительным, что проявляется в неоправданном уменьшении шага и в повторных вычислениях решения (возвратах). Проверка (32) при Зп = 1 осуществляется редко. Однако, как показывают расчеты, использование (32) вместо (29) позволяет сократить вычислительные затраты на 10-20 %.
4. Результаты вычислений
Построенный алгоритм численно исследовался на задаче Коши [6], которая в автономной форме имеет вид
y[ = 0.2(y2 - yi), yi(0) = 0,
y2 = 10yi - (60 - 0.125уз)у2 + 0.125уз, y2(0) = 0,
y3 = 1, Уз(0) = 0,
(34)
а в неавтономной
y1 = 0.2(y2 - yi), yi(0) = 0, у2 = 10yi - (60 - 0.125t)y2 + 0.125t, y2(0) = 0.
(35)
Начальный шаг интегрирования h0 = 0.017, интервал интегрирования [0, 400], решение в конечной точке y1 = 22.2422201, y2 = 27.1107133, y3 = 400, собственные числа матрицы Якоби А1 = -60 ^ -10, Л2 = -0.17 ^ 0.
Из анализа результатов расчетов следует, что построенным алгоритмом решение (34) по сравнению с (35) при любой задаваемой точности осуществляется примерно на 15 % быстрее. Однако в конце интервала интегрирования фактическая точность соответствует задаваемой только при решении (35), в то время как для (34) она почти на порядок хуже. Следует отметить, что в конце интервала интегрирования имеет место достаточно быстрое изменение решения.
При точности расчетов е = 10-3 и грубее построенный алгоритм эффективнее метода Гира (LSODE) [7] примерно 1.5 раза, в то время как при более высокой точности расчетов эффективнее программа LSODE, причем при е = 10-6 более чем в 2 раза. Это естественно, так как в LSODE используются численные формулы до пятого порядка точности включительно, а построенный алгоритм основан на схеме третьего порядка.
Список литературы
[1] RüSENBROCK H. H. Some general implicit processes for the numerical solution of differential equations. Computer, №5, 1963, 329-330.
[2] Современные численные методы решения обыкновенных дифференциальных уравнений. Под ред. Дж. Холла и Дж. Уатта. Мир, М., 1979.
[3] Новиков Е. А., Шитов Ю. А., Шокин Ю. И. Одношаговые безытерационные методы решения жестких систем. Докл. АН СССР, 301, №6, 1988, 1310-1314.
[4] Демидов Г. В., Юматова Л. А. Исследование некоторых аппроксимаций в связи с -устойчивостью полунеявных методов. В "Численные методы механики сплошной среды". ИТПМ СО АН СССР, Новосибирск, 8, №3, 1977, 68-79.
[5] Демидов Г. В., Новиков Е. А. Оценка ошибки одношаговых методов интегрирования обыкновенных дифференциальных уравнений. Там же, 16, №1, 1985, 27-39.
[6] ENRIGHT W. H., Hüll T. E. Comparing numerical methods for the solutions of stiff systems ODE's. BIT, No. 15, 1975, 10-48.
[7] Byrne G.D., HINDMARSH A.C. Stiff ODE solvers: a review of current and comming attractions. J. Comput. Phys., No. 70, 1986, 1-62.
Поступила в редакцию 28 июля 1997 г., в переработанном виде 16 апреля 1998 г.
