Численный метод решения задачи Коши для жестких обыкновенных дифференциальных уравнений на основе многозвенных интерполяционных полиномов Эрмита.

Вычислительные технологии
Том 16, ќ 2, 2011
Численный метод решения задачи Коши
для жестких обыкновенных диеренциальных
уравнений на основе многозвенных интерполяционных
полиномов Эрмита?
А. Ф. Латыпов, О. В. Попик
Институт теоретической и прикладной механики СО АН, Новосибирск, оссия
e-mail:
latypovitam.ns.ru
ассмотрен одношаговый метод решения задачи Коши для жесткой системы
обыкновенных диеренциальных уравнений. Метод основан на представлении
ункций правых частей системы на шаге интерполяционными полиномами Эрмита, заданными значениями ункций в трех точках. Доказана A- и L(?)-устойчивость метода. Дана оценка точности, приводится алгоритм расчета глобальной
ошибки. ассмотрены результаты решений тестовых задач.
Ключевые слова : системы обыкновенных диеренциальных уравнений первого порядка, задача Коши, устойчивость, интерполяция полиномами Эрмита, оценка точности.
В [1? разработано семейство LRM методов численного интегрирования систем
обыкновенных диеренциальных уравнений (ОДУ). Однако метод LRM(0, 0, 0, s) отдельно не представлен и соответственно не изучен. Он прост в реализации,
A-
и
L(?)-
устойчив, имеет четвертый порядок точности по шагу интегрирования. Эти свойства
делают его привлекательным для использования.
1. Постановка задачи, описание метода
ассматривается задача Коши для системы обыкновенных диеренциальных уравнений, число уравнений
n:
y? (x) = f(y, x),
x ? [0, X],
y(0) = y0 .
(1)
Условия существования и единственности решения предполагаются выполненными.
Пусть решение на
[0, xi ]
получено, тогда на отрезке
[xi , xi + h ? X]
представим систе-
му (1) в виде
dy
= hf(y, x) ? ?(y, ?),
d?
Обозначим
y(xi ) = yi ,
?=
(x ? xi )
,
h
? ? [0, 1].
y (x) = y (xi + h?) = y (?) = y? , ? (y (?) , ?) = ? (?) = ?? .
?(y, ?) полиномом Эрмита вида
(2)
Приближенно
представим ункцию
?
абота выполнена при инансовой поддержке Междисциплинарного интеграционного проекта СО
АН ќ 119 и проекта ќ 21.26 программы АН.
78
79
Численный метод решения задачи Коши для жестких...
??(?) = a(? ? 1)(? ? s) + b?(? ? 1) + c?(? ? s)
при условиях
??(0) = ?0 ,
??(s) = ?s ,
??(1) = ?1 ,
s ? [0.5, 1.0).
(3)
Получим
??(?) =
(1 ? ?)(s ? ?)?0 ?(1 ? ?)?s ?(s ? ?)?1
+
+
.
s
s(s ? 1)
(s ? 1)
Подставляя (4) в (2) и интегрируя на отрезках
[0, s], [0, 1],
получим систему из
(4)
2n
алгебраических уравнений
?s = b1 (ys ? y0 ) + a0 ?0 + a1 ??s + ??1 = 0,
?1 = b2 (y1 ? y0 ) ? b2 ?0 + ??s + a2 ??1 = 0,
(5)
где
??s = ?s ? ?0 , ??1 = ?1 ? ?0 .
6(s ? 1)
3 , a = ?s(3s ? 2),
a0 =
, a1 = 2s ?
2
2
2
s
s
6(s ? 1)
b1 = ?
, b2 = 6s(s ? 1).
s3
4
Погрешность метода по шагу интегрирования h удовлетворяет оценке [1? |y(h)?y?(h)| ? h .
Система уравнений (5) решается квазиньютоновским методом [2?. Зададим унк-
I = 0.5(?, ?), ? = {?s , ?1 }. Итерационный процесс поиска решения строится
соотношению d? = ??d? , которое определяет сходимость процесса из любой точки
ционал
по
области притяжения искомого решения, так как
dI = (?, d?) = ?(?, ?)d? = ?2Id?, lim I = 0.
? ??
Итерационный процесс для системы (5) имеет вид
??s i ??s i
?y +
?y = ??i?1
s ?,
?ys s ?y1 1
??1 i ??1 i
?y +
?y = ??i?1
1 ?,
?ys s ?y1 1
(6)
где
?ysi = ysi ? ysi?1 ,
?y1i = y1i ? y1i?1 ,
??s
= b1 E + a1 Ji?1
s ,
?ys
??1
= Ji?1
s ,
?ys
Здесь
J=
??
,
?y
??s
= Ji?1
1 ,
?y1
??1
= b2 E + a2 Ji?1
1 .
?y1
(7)
E единичная матрица. Подставляя (7) в (6), находим i-е приближение решения
?1
D = Ji?1
? (b1 (Ji?1
+ a1 E)(a2 Ji?1
+ b2 E),
1
s )
1
80
А. Ф. Лапытов, О. В. Попик
?1
?y1i = D?1 [(b1 (Ji?1
+ a1 E)?i?1
? ?i?1
s )
1
s ]?,
?1
i?1
?ysi = (Ji?1
+ b2 E)?y1i ? ?i?1
s ) [(a2 J1
1 ? ].
(8)
Начальное приближение задается решением системы ОДУ (2) методом унге Кутты
третьего порядка точности на отрезке
[0, s] и методом Эйлера на отрезке [s, 1]. Из анали-
за устойчивости LRM(0, 0, 0, s)-метода (см. ниже) применительно к системам ОДУ для
параметра
s
целесообразно значение s ? 0.9. Параметр ? ? (0, 1] определяется на кажI i < I i?1 . Итерационный процесс заканчивается при выполнении
дом шаге из условия
i
i
условий I < ?1 , |?y1 |
< ?2 ,
где
?1
?2 заданные величины.
s к единице число обусловленности
и
Замечание. При стремлении
матрицы Яко-
би системы уравнений (5) неограниченно возрастает. Поэтому при обращении матриц
необходимо применять способы повышения точности вычислений. Для приведения матрицы к верхней треугольной орме используются операции умножения и алгебраического сложения. При сложении двух чисел, заданных с абсолютной ошибкой
? (ошибка
округления), абсолютная погрешность результата может максимально только удвоиться. Однако при умножении
гократно
?c ? (|a| + |b|)?.
c = (a + ?)(b + ?)
ошибка результата может возрасти мно-
Для уменьшения ошибки в последнем случае целесообразно
использовать следующий прием. Пусть необходимо найти решение системы
алгебраических уравнений
Ax = b.
Опишем первый шаг: исключение элемента
лю элементы в
(n ? 1)-й
n линейных
и
n-й
an1 .
Определим максимальные по моду-
строках, которые обозначим соответственно
qn?1
и
qn .
Вычисляем
kn?1 = entier(log2 qn?1 ), kn = entier(log2 qn )
и производим нормирование посредством операций
a?n?1,j = an?1,j 2?(kn?1 +1) ,
a?n,j = an,j 2?(kn +1) ,
b?n?1,j = bn?1,j 2?(kn?1 +1) ,
b?n,j = bn,j 2?(kn +1) ,
j = 1, n.
Нормирование проводится в двоичном коде, изменяются лишь порядки чисел, при этом
модули элементов строк матрицы будут меньше единицы. Далее элементы
строки умножаются на
(n, 1)-й
a?n1 , а n-й строки на a?n?1,1 . Сложением
(n ? 1)-й
или вычитанием строк
элемент обращается в ноль. Точно так же обращаются в ноль все поддиаго-
нальные элементы матрицы
Для обращения матрицы
A
A.
Затем решение находится путем обратной прогонки.
необходимо вектор
b
заменить единичной матрицей
E.
2. Устойчивость метода
Для исследования устойчивости метода используется линейное уравнение [3?
z ? (x) = ?z(x),
z(0) = z0 ,
x ? 0,
Re(?)
< 0,
µ = ?h.
(9)
ешение для одношагового метода записывается в виде
z(µ) = z0 q(µ).
Приведем определение
L(?)-устойчивости,
данное в [1?.
(10)
81
Численный метод решения задачи Коши для жестких...
Определение (следуя [3?, с. 122). Одношаговый метод назовем
с параметром
? ? (0, 1),
если метод
A
устойчив и
Ниже будет показано, что величина
случае
L(?)-устойчивость,
как и
L(?)-устойчивым
lim |q(µ)| = ? .
|µ| ??
? может быть выбрана достаточно малой. В этом
L-устойчивость,
является полезным свойством метода
для эективного численного решения сильно жестких систем.
(0, 0, 0, s)-метод A-устойчив
? = (1 ? s)/s при s ? (0.5, 1.0).
Теорема 1. LRM
раметром
при
s ? [0.5, 1.0)
и
L(?)-устойчив
с па-
Доказательство. Обратимся к уравнению (9). Представим комплексное число
?
в тригонометрической орме и обозначим
? = ?[cos(?) + i sin(?)],
c = cos(?) < 0,
d = 2s ? 1,
?? = h?.
(11)
Из (5), (10) получим
? + i?
,
? + i?
s
s
2
2
? +?
G(??)
|q(µ)| =
,
2
2 =
H(??)
? +?
q(µ) =
(12)
где
G(??) =
4
P
gi ??i ,
H(??) =
i=0
g0
g1
g2
g3
g4
4
P
pi ??i ,
i=0
= 144, p0 = 144,
= c(144 ? 48d), p1 = c(?144 ? 48d),
= 48c2 + d2 ? 48c2 d + 12, p2 = 48c2 + d2 + 48c2 d + 12,
= 4d2 c ? 16dc + 12c, p3 = ?4d2 c ? 16dc ? 12c,
= 1 ? 2d + d2 , p4 = 1 + 2d + d2 .
азность между знаменателем и числителем в (12)
H ? G = ? c??( 288 + 8??2 d2 + 24??2 ) + d??2 (96c2 + 4??2 )
d ? 0 (c < 0 по условию задачи (9)). Так как
s ? 1/2 для любого ??, откуда следует справедливость
положительна при достаточном условии
G > 0, H > 0,
то
|q(??)| < 1
при
первого утверждения теоремы. Из приведенных соотношений следует
lim |q(??)| =
????
r
g4
1?d
1?s
=
=
,
p4
1+d
s
что доказывает второе утверждение теоремы.
3. асчет локальной и глобальной ошибки
Вычисление ошибки
?y = y ? y?
на шаге
h
производится путем следующей процедуры.
Исходная система
dy(?)
= ?(y, ?),
d?
? ? [0, 1],
y(0) = y0 .
(13)
82
А. Ф. Лапытов, О. В. Попик
Приближенная система
dy?(?)
= ??(?),
d?
При
|?y| << 1
с погрешностью
o(|?y|)
? ? [0, 1],
y?(0) = y0 .
(14)
справедлива оценка
?(y, ?) ? ?(y?, ?) + J?(?)?y,
где
J?(?) = hJ(?).
(15)
Из (13), (14), (15) получим уравнение для вычисления ошибок
d?y
= Q(?) + J?(?)?y = R(?y, ?),
d?
Ошибка аппроксимации ункции
?(t)
Q(?) = ?(y?, ?) ? ??(?).
интерполяционным полиномом Эрмита
(16)
Hm (t)
равна [4?
?(m+1) (?)
?(t) ? Hm (t) =
?(t),
(m + 1)!
?(t) = (t ? t0 )?0 (t ? t1 )?1 · · · (t ? tn )?n ,
m + 1 = ?0 + ?1 + ... + ?n .
(17)
Q(?) будем использовать
?0 = ?s = ?1 = 1, m + 1 = 3 и
представление (17). В случае
Для приближенного значения
LRM(0, 0, 0, s)- метода
Q(?) ? Q?(?) = C1 ?(?),
Коэициент
C1
?(?) = ?(? ? s)(? ? 1).
(18)
определим из условия
Q?(?? ) = Q(?? ),
где
??
?(?), ?(?? ) = max ?(?), ? ? (0, 1).
q
(s + 1) ? (s + 1)2 ? 3s
?? =
,
3
соответствует максимуму
C1 =
?(y?? ) ? ??(?? )
,
?? (?? ? s)(?? ? 1)
Получим
y?? = y?(?? ).
(19)
Аппроксимация Якобиана. Якобиан аппроксимируем квадратичной ункцией
от
?
J?(?) = (1 ? ?)A + ?B + ?(1 ? ?)D
при условиях
J?(0) = J?0 , J?(1) = J?1 , J?(s) = J?s .
A = J?0 ,
B = J?1 ,
(20)
Получим
J?1 ? J?s
1
D = [(J?1 ? J?0 ) ?
].
s
1?s
Все ункции в (16) определены, и интегрирование системы квазилинейных диеренциальных уравнений (16) производится методом LR(1, 1) [1?, модиицированным на
случай неавтономной системы, пятого порядка погрешности при начальном значении
?y0 ,
равном глобальной ошибке интегрирования на предыдущем шаге.
83
Численный метод решения задачи Коши для жестких...
Модиикация LR
(1, 1)-метода
для интегрирования системы квазилиней-
ных ОДУ. В данном случае интерполяционный полином Эрмита записывается следу-
ющим образом:
R(?y, ?) ? g(?) = R0 ?0,0 (?) + R? ?0 ?0,1 (?)+R1 ?1,0 (?) + R??1 ?1,1 (?),
?0,0 (?) = 1 ? 3? 2 + 2? 3 ,
?1,0 (?) = 3? 2 ? 2? 3 ,
R? ? (?y, ?) =
?0,1 (?) = ? ? 2? 2 + ? 3,
?1,1 (?) = ?? 2 + ? 3 ,
dQ(?) dJ?(?)
+
?y + J?(?)(Q(?) + J?(?)?y).
d?
d?
(21)
Интегрируя (16) с учетом (21), получаем систему линейных алгебраических уравнений для определения глобальной ошибки на шаге
h
1
1
1
1
?y1 = R0 + R? 0 + R1 ? R?1 R0 = J?0 ?y0 ,
2
12
2
12
R1 = J?1 ?y1 ,
R? 0 = C1 s + [(?A + B + D) + J?20 ]?y0 ,
R? 1 = C1 (1 ? s) + [(?A + B ? D) + J?21 ]?y1 .
?yloc =
|?yloc |
?y1 ? ?y0 . Задается какая-либо монотонно убывающая ункция K(u), u =
? ,
K(1) = 1, ? заданная локальная точность. Если на полученном решении K(ui ) > 1,
hi . В противном случае вычисто расчет повторяется с уменьшенным шагом hi =
K(ui )
hi .
ляется величина следующего шага hi+1 =
K(ui )
Для регулирования шага интегрирования
h
(22)
используется локальная ошибка
4. Тестовые расчеты
Для тестовых расчетов взяты следующие системы диеренциальных уравнений (?i собственные числа матрицы Якоби при
t = t0 ).
1. Жесткая система с пограничным слоем на левом конце
y?1 = ?10 y?2 + 3000(1 ? y12),
y?2 = 0.04(1 ? y2 ) ? (1 ? y1 )y2 + 10?4 (10?1 ? y12),
0 ? t ? 2.6, y1 (0) = 1, y2 (0) = 0,
?1 = ?6000, ?2 = ?0.04.
2. Жесткая система с периодически повторяющимся пограничным слоем
y?1 = ?2000y1 + 1000y2 + 1 + sin(10t),
y?2 = y1 ? y2 ,
0 ? t ? 4, y1 (0) = y2 (0) = 0,
?1 = ?2000, ?2 = ?0.5.
84
А. Ф. Лапытов, О. В. Попик
3. Жесткая система с пограничным слоем на левом конце для первых двух компонент и на правом конце для третьей компоненты вектора решения
y?1 = ?(55 + y3 )y1 + 65y2,
y?2 = 0.0785(y1 ? y2 ),
y?3 = 0.1y1,
0 ? t ? 500, y1 (0) = y2 (0) = 1, y3 (0) = 0,
?1 = ?55, ?2,3 = 0.00625 ± 0.01i.
4. Слабо жесткая задача с пограничным слоем на правом конце. Сильная чувствительность к возмущениям начальных условий и очень быстрое возрастание на правом
конце (тест Троеша [5?):
y?1 = y2 ,
y?2 = sh(y1 ),
0 ? t ? 10,
y1 (0) = 0, y2 (0) = 3.585 · 10?4 .
асчеты выполнены LRM(0, 0, 0, s)-методом, а также неявным методом унге Кутты пятого порядка точности (К5), использующим вторую производную от решения [6?
(программная версия 1995 г.), и стандартным методом унге Кутты четвертого порядка точности (КСТ). езультаты решения сведены в таблицу, где представлены шаги
интегрирования: минимальный
hmin ,
максимальный
hmax .
Эти данные характеризуют
эективность используемого способа автоматического выбора шага интегрирования
по локальной точности. Затрачиваемые ресурсы характеризуются числом вычислений
правых частей системы
kf , ? задаваемая локальная точность интегрирования, ?yglob получаемая глобальная точность. лобальная точность решения задач методами КСТ
и К5 оценивалась сравнением с решениями, полученными методом LRM(0, 0, 0, 0.9)
(столбец
?y ).
Из приведенных в таблице данных видно, что эективность рассматриваемого метода по критерию быстродействия, которому отвечает количество вычислений правых
Сравнительная таблица решений тестовых задач различными методами
Задача Метод
1
КСТ
К5
LRM(000, 0.9)
2
КСТ
К5
LRM(000, 0.9)
3
КСТ
К5
LRM(000, 0.9)
4
КСТ
К5
LRM(000, 0.9)
hmin
8.0Е-05
1.7Е-04
1.0E-04
2.7Е-05
1.7Е-04
1.0E-05
1.7Е-03
1.0Е-03
0.12
3.8Е-04
1.4Е-04
1.0E-04
hmax
1.9Е-03
1.3
1.0
2.0Е-03
1.0Е-01
0.4
1.9Е-01
12
49
1.8Е-01
2.2Е-01
1.0
kf
51352
53
70
88864
586
553
14261
2405
1107
3588
2405
1330
?
1.0Е-07
1.0Е-07
1.0Е-07
1.0Е-07
1.0Е-07
1.0Е-07
1.0E-07
1.0E-07
1.0Е-07
1.0E-07
1.0E-06
1.0Е-07
?yglob
1.0Е-06
1.0Е-06
1.0Е-08
1.0Е-03
?y
3.9Е-06
1.0Е-09
5.7Е-07
1.0Е-07
4.7E-07
1.0E-07
7.3E-05
1.0E-04
Численный метод решения задачи Коши для жестких...
85
частей системы, для разных задач различна. Для задачи 1 в методе К5 количество
вычислений правых частей минимально, однако для других задач видны преимущества
нового метода. Представленный метод позволяет получать решения с высокой точностью. Это свойство особенно важно для систем с наличием неустойчивости решений к
возмущениям начальных значений азовых переменных, к которым, например, относятся системы уравнений движения небесной механики.
Для нежестких систем и систем с умеренной жесткостью предпочтительно значение
s = 0.5.
Список литературы
[1?
[2?
[3?
Численные методы решения задачи Коши для жестких систем обыкновенных диеренциальных уравнений на основе многозвенных интерполяционных полиномов Эрмита // Журн. вычисл. математики и матем. изики. 2007.
Т. 47, ќ 2. С. 234244.
Латыпов А.Ф., Никуличев Ю.В.
Бахвалов Н.С., Жидков Н.П., Кобельков .М.
Современные численные методы решения обыкновенных диеренциальных уравнений / ед. Дж. Холл и Дж. Уатт. М.: Мир, 1979.
[4?
Березин И.С., Жидков Н.П.
[5?
Roberts S.M., Shipmen J.S.
[6?
Численные методы. М.: Наука, 1987.
Методы вычислений. М.: Наука, 1970.
Solution of Troesh's two-point boundary value problem by a
ombination of tehniques // J. Comput. Phys. 1972. Vol. 10, No. 2.
ешение обыкновенных диеренциальных уравнений. Жесткие
и диеренциально-алгебраические задачи. М.: Мир, 1999.
Хариер Э., Ваннер .
Поступила в редакцию 14 октября 2009 г.,
с доработки 29 сентября 2010 г.