CC BY
50
Вычислительные технологии
Том 16, № 2, 2011
Численный метод решения задачи Коши для жестких обыкновенных дифференциальных уравнений на основе многозвенных интерполяционных
полиномов Эрмита*
А. Ф. Латыпов, О. В. Попик Институт теоретической и прикладной механики СО РАН, Новосибирск, Россия
e-mail: latypov@itam.nsc.ru
Рассмотрен одношаговый метод решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений. Метод основан на представлении функций правых частей системы на шаге интерполяционными полиномами Эрмита, заданными значениями функций в трех точках. Доказана А- и £(£)-устойчи-вость метода. Дана оценка точности, приводится алгоритм расчета глобальной ошибки. Рассмотрены результаты решений тестовых задач.
Ключевые слова: системы обыкновенных дифференциальных уравнений первого порядка, задача Коши, устойчивость, интерполяция полиномами Эрмита, оценка точности.
В [1] разработано семейство 1.НМ — методов численного интегрирования систем обыкновенных дифференциальных уравнений (ОДУ). Однако метод 1ЛМ(0, 0, 0, в) отдельно не представлен и соответственно не изучен. Он прост в реализации, А- и Ь(6)~ устойчпв, имеет четвертый порядок точности по шагу интегрирования. Эти свойства делают его привлекательным для использования.
1. Постановка задачи, описание метода
Рассматривается задача Коши для системы обыкновенных дифференциальных уравнений, число уравнений и:
у'(х) = {(у,х), х € [0,X], у(0) = Уо. (1)
Условия существования и единственности решения предполагаются выполненными. Пусть решение на [0,х^] получено, тогда на отрезке [х^х^ + к < X] представим систему (1) в виде
^ = М(у,х) = Ф(у,0, у(хг) = уг, ^ [0,1]. (2)
Обозначим у (х) = у (х^ + к£) = у (£) = у^, Ф (у (£), £) = Ф (£) = Ф^. Приближенно представим функцию Ф(у,£) полиномом Эрмита вида
* Работа выполнена при финансовой поддержке Междисциплинарного интеграционного проекта СО РАН № 119 и проекта № 21.26 программы РАН.
при условиях
Получим
Ф(С) = a(£ - 1)(£ - s) + b£(£ - 1) + е£(С - s)
Ф(0) = Фо, Ф(s) = Фs, Ф(1) = Ф1, s G [0.5,1.0). (3)
m =-;-+ ф-i) + (S-i) • (4)
Подставляя (4) в (2) и интегрируя на отрезках [0,s], [0,1], получим систему из 2n алгебраических уравнений
Ps = bi(ys - Уо) + «оФо + «1АФ5 + АФ1 = 0,
Ч>1 = &2(yi - Уо) - Ь2Фо + АФs + «2АФ1 = 0, (5)
где
АФs = Фs - Фо, АФ1 = Ф1 - Фо.
Ü0 = -- -~> al = а2 = —s(3s — 2),
bi = - 6(gg71}, h = 6ф-1). Погрешность метода по шагу интегрирования h удовлетворяет оценке [1] |y (h)-у (h) | ~
h4.
Система уравнений (5) решается квазиньютоновским методом [2]. Зададим функционал I = 0.5(р,р), р = {ps,p1}. Итерационный процесс поиска решения строится по соотношению dp = -pdr, которое определяет сходимость процесса из любой точки области притяжения искомого решения, так как
dl = (р, dp) = -(р, p)dr = -2Idr, lim I = 0.
т ^го
Итерационный процесс для системы (5) имеет вид
(а)
. Л -г
öyi = у1~ у1~\ = у\ - yi~\ J =
где
9(р* - h F 4- л Г-1 9(р* - Г"1 —— — bitj + aiJs , -— — ,
dys dy1
dP1 Ti-1 dP1 un , „ Ti-1
_ j;-1, 7P = 62E + a2Jri. (7)
dys dy1
Здесь E — единичная матрица. Подставляя (7) в (6), находим г-е приближение решения
D = J1-1 - (bJ-1)-1 + ацЕ)^-1 + fcE),
iyj = d-1[(61(j;-1)-1 + aEm-1 -
iyj = №-1)-1[М-1 + b2E)iyi - ^i-1r ]. (8)
Начальное приближение задается решением системы ОДУ (2) методом Рунге — Кутты третьего порядка точности на отрезке [0, s] и методом Эйлера на отрезке [s, 1], Из анализа устойчивости LRM(0, 0, 0^)-метода (см, ниже) применительно к системам ОДУ для параметра s целесообразно значение s « 0.9. Параметр т G (0,1] определяется на каждом шаге из условия Iг < Iг-1. Итерационный процесс закапчивается при выполнении условий Iг < е1; |$yi | < е2, где е1 и е2 — заданные величины.
s
би системы уравнений (5) неограниченно возрастает. Поэтому при обращении матриц необходимо применять способы повышения точности вычислений. Для приведения матрицы к верхней треугольной форме используются операции умножения и алгебраического сложения. При сложении двух чисел, заданных с абсолютной ошибкой е (ошибка округления), абсолютная погрешность результата может максимально только удвоиться. Однако при умножении c = (а + е)(Ь + е) ошибка результата может возрасти многократно < (|а| + |6|)е. Для уменьшения ошибки в последнем случае целесообразно использовать следующий прием. Пусть необходимо найти решение системы n линейных алгебраических уравнений Ax = b.
Опишем первый шаг: исключение элемента ani, Определим максимальные по модулю элементы в (n — 1)-й и n-й строках, которые обозначим соответственно qn-i и qn. Вычисляем
kn-i = entier (log2 qn-i), kn = entier (kg qn ) и производим нормирование посредством операций
a 1 • = а т .2-(fcn-1+i) а • = а .2-(кп+1)
h Л ■ — h л .O-i.kn-1+l) T. . — h o-(fcn+l) A — TT7 un-i,J un-l,]^ 1 Un,J un,JL' 1 J
Нормирование проводится в двоичном коде, изменяются лишь порядки чисел, при этом модули элементов строк матрицы будут меньше единицы. Далее элементы (n — 1)-й строки умножаются на âni, а n-й строки — на an-i,i. Сложением или вычитанием строк (n, 1)-й элемент обращается в ноль. Точно так же обращаются в ноль все поддиаго-
A
A b E
2. Устойчивость метода
Для исследования устойчивости метода используется линейное уравнение [3]
z'(x) = Az(x), z(0) = zo, x > 0, Re(A) < 0, ^ = Ah. (9)
Решение для одношагового метода записывается в виде
z(u) = zoq(u). (10)
Приведем определение ¿(^-устойчивости, данное в [1].
Определение (следуя [3], с, 122), Одношаговый метод назовем ¿(0)-устойчивым с параметром ö £ (0,1), если метод A устойчив и lim |q(p) = ö,
ö
случае £(0)-устойчивость, как и L-уетойчивоеть, является полезным свойством метода для эффективного численного решения "сильно жестких" систем.
Теорема 1. LRM(0, 0, 0^)-метод A-устойчив при s £ [0.5,1.0) и Ь(0)-устойчив с параметром ö = (1 — s)/s при s £ (0.5,1.0),
Доказательство. Обратимся к уравнению (9), Представим комплексное число Л в тригонометрической форме и обозначим
Л = p[cos(0) + i sin(0)j,
c = cos(0) < 0, d = 2s — 1, p = hp. (11)
Из (5), (10) получим
, s a + ¿в
QW = , -г, Y + го
где
4 4
g(p) = e 9ipi, H (p) = E pipi, i=0 i=0
go = 144, po = 144,
gi = c(144 - 48d), pi = c(-144 - 48d),
g2 = 48c2 + d2 - 48c2d + 12, p2 = 48c2 + d2 + 48c2d + 12,
g3 = 4d2c - 16dc + 12c, p3 = -4d2c - 16dc - 12c,
g4 = 1 - 2d + d2, p4 = 1 + 2d + d2.
Разность между знаменателем и числителем в (12)
H - G = - cp( 288 + 8p2d2 + 24р2) + dp2 (96c2 + 4р2)
положительна при достаточном условии d > 0 (c < 0 по условию задачи (9)), Так как G > 0, H > 0, то |q(p)| < 1 при s > 1/2 для любо го р, откуда следует справедливость первого утверждения теоремы. Из приведенных соотношений следует
[gi 1 - d 1 - s
lim q(p) = J—= —— =-,
р^ж у p4 1 + d s
что доказывает второе утверждение теоремы, 3. Расчет локальной и глобальной ошибки
Вычисление ошибки = y - у на шаге h производится путем следующей процедуры. Исходная система
^ = Ф(у,е), ее [0,1], у(о) = уо- (13)
Приближенная система
)
Ф(е), е е [0,1], у(0) = Уо. (14)
^е
При |$у| << 1 с погрешноетью о(|$у|) справедлива оценка
ф(у,е)« ф(у,е) + ЗД^у, (15)
где 1(е) = ЛЛ(£). Из (13), (14), (15) получим уравнение для вычисления ошибок ^у
¿е
С(е) + 5(е)*у = ,е), С(е) = ф(у,е) - ф(е). (16)
Ошибка аппроксимации функции ф(г) интерполяционным полиномом Эрмита Нт(г) равна [4]
вд = (г-гоГ(г-^Г ••• (г-ьп)ап,
т +1 = ао + а + ... + а„. (17)
Для приближенного значения СС(е) будем использовать представление (17), В случае ЬЕМ(0, 0, 0, з)- метода а0 = а 3 = а1 = 1, т +1 = 3и
С(е)« С(е) = ад(е), вд = е(е - *хе -1). (18)
Коэффициент С1 определим из условия
СС (е*) = С(е*),
где е* соответствует максимуму П(е), П(е*) = тах П(е), е е (0,1). Получим
(з + 1) -у (з + 1)2 - 3з е* = з ,
Аппроксимация Якобиана. Якобиан аппроксимируем квадратичной функцией
е
5(е) = (1 - е)а+ев+е(1 - ер (20)
при условиях 5(0) = 50, 5(1) = 51, 5(з) = 5 Получим
А = 30, В = 3Ъ Б = -[(5!-
з 1 - з
Все функции в (16) определены, и интегрирование системы квазилинейных диффе-
(1, 1)
случай неавтономной системы, пятого порядка погрешности при начальном значении $у0, равном глобальной ошибке интегрирования на предыдущем шаге.
Модификация Ы1(1,1)-метода для интегрирования системы квазилинейных ОДУ. В данном случае интерполяционный полином Эрмпта записывается следующим образом:
К(£у,£) - g(е) = Ко^о,о(£) + к'^о,1(е)+ВД1,о(£) + к/^мФ,
^ме) = 1 - зе2 + 2£3, ^0,1(0 = е - 2£2+е3,
) = зе2 - 2е3, ^,1(£) = -е2+е3, 0 = ^ + ^ + 3(0(0(0 + ЭД^У)- (21)
Интегрируя (16) с учетом (21), получаем систему линейных алгебраических уравнений для определения глобальной ошибки на шаге к
¿У1 = ^Ио + ^К-'о + 7^1 - К'11*0 = Зо^Уо, = З^Уь
К'о = С^ + [(-А + В + Б) + ^уо, БЛ = С1(1 - в) + [(-А + В - Б) + 5?]£у1. (22)
Для регулирования шага интегрирования к используется локальная ошибка 6у1,
ос
— ^Уо- Задается какая-либо монотонно убывающая функция К (и), и = ^ I
К(1) = 1, е — заданная локальная точность. Если на полученном решении К(м) > 1,
к
то расчет повторяется с уменьшенным шагом к» = ¡{(и-)' ^ ПР0ТИВН0М случае вычнс-
к.
ляется величина следующего шага //... | — —-
К{щУ
4. Тестовые расчеты
Для тестовых расчетов взяты следующие системы дифференциальных уравнений (А^ собственные числа матрицы Якоби при £ = ¿о).
1, Жесткая система с пограничным слоем на левом конце
У1 = -10 У2 + 3000(1 - у2),
У2 = 0.04(1 - У2) - (1 - ш)У2 + 10-4(10-1 - у2),
0 < £ < 2.6, У1(0) = 1, У2(0) = 0,
А1 = -6000, А2 = -0.04.
2, Жесткая система с периодически повторяющимся пограничным слоем
У1 = -2000/1 + 1000/2 + 1 + яп(10£),
У/2 = У1 - У2,
0 < £ < 4, /1(0) = /2(0) = 0, А1 = -2000, А2 = -0.5.
3, Жесткая система с пограничным слоем на левом конце для первых двух компонент и на правом конце для третьей компоненты вектора решения
у/1 = -(55 + уз)У1 + 65^2, /2 = 0.0785(/1 - /2), У/з = 0Т/ъ
0 < г < 500, /1(0) = /2(0) = 1, /з(0) = 0, Л1 = -55, Л2,з = 0.00625 ± 0.01г.
4, Слабо жесткая задача с пограничным слоем на правом конце. Сильная чувствительность к возмущениям начальных условий и очень быстрое возрастание на правом конце (тест Троеша [5]):
2/1 = /2,
у2 =
0 < г < 10,
/1(0) = 0, /2(0) = 3.585 ■ 10-4.
Расчеты выполнены 1ЛМ(0, 0, 0, з)-методом, а также неявным методом Рунге — Кут-ты пятого порядка точности (РК5), использующим вторую производную от решения [6] (программная версия 1995 г.), и стандартным методом Рунге —Кутты четвертого порядка точности (РКСТ), Результаты решения сведены в таблицу, где представлены шаги интегрирования: минимальный Л,т;п, макспмальный Л.тах, Эти данные характеризуют эффективность используемого способа автоматического выбора шага интегрирования по локальной точности. Затрачиваемые ресурсы характеризуются числом вычислений правых частей системы к/, £ — задаваемая локальная точность интегрирования, $/ё10ь _ получаемая глобальная точность. Глобальная точность решения задач методами РКСТ и РК5 оценивалась сравнением с решениями, полученными методом 1ЛМ(0, 0, 0, 0.9) (столбец #/),
Из приведенных в таблице данных видно, что эффективность рассматриваемого метода по критерию быстродействия, которому отвечает количество вычислений правых
Сравнительная таблица решений тестовых задач различными методами
Задача Метод ъ ■ 1 <-111111 Ъ 1 'тах е 6у
1 РКСТ 8.0Е-05 1.9Е-03 51352 1.0Е-07 _ 3.9Е-06
РК5 1.7Е-04 1.3 53 1.0Е-07 _ 1.0Е-09
Ы1М(000,0.9) 1.0Е-04 1.0 70 1.0Е-07 1.0Е-06 _
2 РКСТ 2.7Е-05 2.0Е-03 88864 1.0Е-07 _ 5.7Е-07
РК5 1.7Е-04 1.0Е-01 586 1.0Е-07 _ 1.0Е-07
Ы1М(000,0.9) 1.0Е-05 0.4 553 1.0Е-07 1.0Е-06 _
3 РКСТ 1.7Е-03 1.9Е-01 14261 1.0Е-07 _ 4.7Е-07
РК5 1.0Е-03 12 2405 1.0Е-07 _ 1.0Е-07
Ы1М(000,0.9) 0.12 49 1107 1.0Е-07 1.0Е-08 _
4 РКСТ 3.8Е-04 1.8Е-01 3588 1.0Е-07 _ 7.3Е-05
РК5 1.4Е-04 2.2Е-01 2405 1.0Е-06 _ 1.0Е-04
Ы1М(000, 0.9) 1.0Е-04 1.0 1330 1.0Е-07 1.0Е-03 _
частей системы, для разных задач различна. Для задачи 1 в методе РК5 количество вычислений правых частей минимально, однако для других задач видны преимущества нового метода. Представленный метод позволяет получать решения с высокой точностью, Это свойство особенно важно для систем с наличием неустойчивости решений к возмущениям начальных значений фазовых переменных, к которым, например, относятся системы уравнений движения небесной механики.
Для нежестких систем и систем с умеренной жесткостью предпочтительно значение s = 0.5.
Список литературы
[1] Латыпов А.Ф., Никуличев Ю.В. Численные методы решения задачи Коши для жестких систем обыкновенных дифференциальных уравнений на основе многозвенных интерполяционных полиномов Эрмита // Журн. вычисл. математики и матем. физики. 2007. Т. 47, № 2. С. 234-244.
[2] Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Наука, 1987.
[3] Современные численные методы решения обыкновенных дифференциальных уравнений / Ред. Дж. Холл и Дж. Уатт. М.: Мир, 1979.
[4] Березин И.е., Жидков Н.П. Методы вычислений. М.: Наука, 1970.
[5] Roberts S.M., Shipmen J.S. Solution of Troesch's two-point boundary value problem by a combination of techniques //J. Comput. Phvs. 1972. Vol. 10, No. 2.
[6] Хариер Э., Bahhep Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М.: Мир, 1999.
Поступила в редакцию Ц октября 2009 г., с доработки — 29 сентября 2010 г.
