ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2007, том 47, < 2, с. 234-244
УДК 519.622
ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧИ КОШИ ДЛЯ ЖЕСТКИХ СИСТЕМ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ НА ОСНОВЕ МНОГОЗВЕННЫХ ИНТЕРПОЛЯЦИОННЫХ
ПОЛИНОМОВ ЭРМИТА1}
© 2007 г. А. Ф. Латыпов, Ю. В. Никуличев
(630090 Новосибирск, ул. Институтская, 4/1, Ин-т теор. и прикл. механ. СО РАН)
e-mail: latypov@itam.nsc.ru.
Поступила в редакцию 01.03.2005 г. Переработанный вариант 26.01.2006 г.
Построены семейства A-, L- и L(5)-устойчивыx методов решения задачи Коши для системы обыкновенных дифференциальных уравнений (ОДУ). Дано определение L(5)-устойчивости метода с параметром 5, 5 е (0, 1). Методы основаны на представлении правых частей системы ОДУ на шаге h в виде, соответственно, двух- и трехточечных интерполяционных полиномов Эрмита. Приведены сравнительные результаты решений тестовых задач. На основе многозвенных интерполяционных полиномов Эрмита получены формулы для вычисления определенных интегралов. Даны оценки точности. Библ. 10. Табл. 1.
Ключевые слова: системы обыкновенных дифференциальных уравнений первого порядка, задача Коши, устойчивость, интерполяция полиномами Эрмита, оценка точности.
1. ПОСТАНОВКА ЗАДАЧИ
Рассматривается задача Коши для системы ОДУ
у'(х) = f (у), х е [0, X], Уо = у(0). (1.1)
Предполагается, что выполнены условия существования и единственности решения задачи. Пусть для х е [0, хг] решение получено. На отрезке [х, + h < X] представим систему (1.1) в следующем виде:
dy= hf (у) = Ф(%), у(хг) = уг, % = (х - хг)Ш, % е [0, 1 ]. (1.2)
Описываемые ниже методы основаны на представлении правых частей системы ОДУ на шаге h в виде двух- и трехточечных интерполяционных полиномов Эрмита. В Приложениях приводятся необходимые соотношения.
2. СЕМЕЙСТВО ЬЯ-МЕТОДОВ
Зададимся целыми положительными величинами L и R. Предполагается, что функции правых частей на решении имеют порядок гладкости по аргументу х не ниже Q = L + R + 2. Представим их на шаге h интерполирующими полиномами Эрмита по двум точкам (см. Приложение 1, (П1.3)):
L Я
Ф (%) = £ Ф0\о, + X Ф1"Ч г(%) . (2.1)
1 = 0 г=0
Верхним индексом в круглых скобках обозначены соответствующие производные по независимой переменной. Коэффициенты во второй сумме (2.1) зависят от значений фазовых переменных в правой точке, подлежащих определению. Проинтегрировав уравнения (1.2) на отрезке [0, 1] при
Работа выполнена при частичной финансовой поддержке РФФИ (код проекта 02-07-90359) и Междисциплинарного интеграционного проекта СО РАН < 119 2003 г.
Ф© ~ Ф ©, получим следующую систему алгебраических уравнений для определения приближенного решения на шаге И (у (И) = у (х, + И)):
У (h) = уг + X So, i( 1) Ф0г) + X Si. r(1 )ф1Г) •
(2.2)
i = o
r = 0
Методы решения задачи Коши для систем ОДУ, основанные на пошаговом представлении функций правых частей в форме (2.1) и решении систем алгебраических уравнений (2.2), названы ЬЯ-методами (см. [1]). Метод, определенный заданными значениями параметров Ь и Я будем называть ЬЯ(Ь, Я)-методом. Погрешность ЬЯ-методов по шагу интегрирования И согласно (П1.8) удовлетворяет оценке
|у(И) - У(И)|~ Ив + 1, т.е. порядок погрешности ЬЯ-метода равен Q.
Для системы линейных ОДУ с постоянными коэффициентами
у' = Ау + Ь, у(0) = Уо, х е [0, X] (2.3)
производные вектор-функции у(х) равны
y(k) = Aky + Ak -1 b •
(2.4)
Для получения приближенного решения системы (2.3) на шаге И, используя (2.4) в системе (2.2), получим систему линейных алгебраических уравнений для определения решения у (И):
E- XSi,r( 1)A
—.r + 1
У (h) =
E + X S o, i( 1) A
-¡i +1
i = o
У.
+
X So, i( 1) Al + X S1, r (1) A
Li = o
hb.
(2.5)
(Здесь А0 = Е, А = ИкАк, Е - единичная матрица.)
Теорема 1. ЬЯ-Метод А-устойчив при Ь < Я < Ь + 2 и Ь-устойчив при Ь < Я < Ь + 2. Доказательство. Для исследования устойчивости метода используется линейное уравнение (см. [2])
г'(х) = Хг(х), г(0) = г0, х > 0, Яе(Х) < 0. (2.6)
Применяя схему (2.5), получаем
^(ц) = ¿oq(ц), q(ц) =
R
-1
1 + XSo,i(1 )Ц +1 1- XS1,r(1 )Ц +1
V i = o
/V r = o
, ц = Xh.
(2.7)
Выражение (2.7) с учетом (П1.6) представляет собой аппроксимацию Падэ для функции ехр(ц) (в обозначениях, принятых в [3, с. 24, [(Ь + 1)/(Я + 1)]]). Согласно теореме об А-устойчиво-сти аппроксимации Падэ (см. [4]), ЬЯ-методы А-устойчивы при значениях параметров, удовлетворяющих неравенствам Ь < Я < Ь + 2. Это доказывает первую часть теоремы.
Представим комплексное число Х в форме (в прямые скобки берется модуль комплексного числа) Х = р[со8(ф) + / 8т(ф)], со8(ф) < 0, р = Ир. Тогда имеем
4 (ц)| = \д (р)| = 7[ Яе (д )]2 + [ 1т( д )]2 = Т4+45/(43+4Г),
Ь Ь
41 = 1 + X50,1( 1 )Р1 +1С05 [ф(I +1 )], 42 = X 50,1( 1 )р' +151П [ф(I +1 )],
(2.8)
i=o
R
i=o
q3 = 1- X S1, m( 1 )pm + 1cos [ф( m +1)], q4= X S1, m( 1 )pm +1sin [ф( m +1)] •
m = o
m=o
Поскольку при R > L степень полинома в знаменателе в (2.8) выше, чем в числителе, то lim |q(р)| = 0, что доказывает вторую часть теоремы.
L
R
R
L
- L
R
r=o
r=o
R
Приведем вычислительные схемы двух методов из описанного семейства.
ЬЩ1,1)-Метод. Положив L = 1, Я = 1 в (2.2), получим Л-устойчивый одношаговый метод 4-го порядка погрешности:
у(h) = уг + (Ф0 + Ф1)/2 + (Ф01) - Ф11))/12.
ЬЩ0,1)-Метод. Положив L = 0, Я = 1 в (2.2), получим L-устойчивый одношаговый метод 3-го порядка погрешности:
у(h) = уг + (Ф0 + 2Ф1)/3- Ф11)/6.
3. СЕМЕЙСТВО ЬМЯ-МЕТОДОВ
Зададимся целыми положительными величинами L, М и Я. Предположив, что функции правых частей на решении имеют порядок гладкости Q = L + М + Я + 3, представим их на шаге h интерполирующими полиномами Эрмита по трем точкам (Приложение 2, (П2.3)):
L М Я
Ф (%) = X Ф01)Х0,1(%) + X Ф^ - (%) + X Ф1г)Х1, г(%). (3.1)
1=0 т=0 г=0
Проинтегрировав полиномы (3.1) на отрезках [0, 5], и [0, 1], 5 е (0, 1), получим следующую систему алгебраических уравнений:
L М Я
у (= уг + X /0,1 (5 ) Ф01) + X ^ т(5 ) Ф5т) + X 71, г(5 ) Ф^,
1 = 0 т = 0 г = 0
L М Я V • /
у^) = уг + X/0,1( 1)Ф01) + X т( 1)Ф5т) + X 71. г( 1)Ф1г),
1 = 0 т = 0 г = 0
где
% % %
/0,1(%) = |%0, l(%)d%, 15, т(%) = т(%)d%, Л, г(%) = г(%)d%, (3.3)
0 0 0
а подынтегральные функции даются формулами (П2.3). Методы решения задачи Коши для систем ОДУ, основанные на пошаговом представлении функций правых частей в форме (3.1) и решении систем алгебраических уравнений (3.2), назовем ЬМЯ-методами. Метод, определенный заданными значениями параметров L, М и Я, будем называть ЬМЯ(Ь, М, Я)-методом. Погрешность ЬМЯ-методов по шагу интегрирования h согласно (П2.4) удовлетворяет оценке
|у(h) - у(h)|~ ^ +1.
Для линейных систем ОДУ производные от функций правых частей определяются формулами (2.4), что позволяет применять для таких систем ЬЯ- и ЬМЯ-методы с произвольными значениями параметров L и Я, которые определяются требуемой точностью решения. При этом решение на шаге h определяется из системы линейных алгебраических уравнений. Решение линеаризованной по у системы (1.2) можно использовать для получения первого приближения в итерационных методах решения алгебраических систем (2.2) и (3.2).
Определение (см. [2, с. 122]). Одношаговый метод назовем L(5)-устойчивым с параметром 5, 5 е (0, 1), если метод Л-устойчив и если для задачи (2.6) будет Нт ^(ц)| = 5.
Ниже на примерах двух алгоритмов из данного семейства будет показано, что для определенных численных методов величина 5 может быть выбрана сколько угодно малой и в этом случае L(5)-устойчивость, как и L-устойчивость является полезным свойством метода для эффективного численного решения "сильно жестких" (см. [5, с. 133]) систем.
Рассмотрим два метода, порождаемых семейством ЬМЯ-методов при значениях параметров, соответственно, L = 1, М = 0, Я = 1 и L = 1, М = 1, Я = 1.
ЬМЩ1,0,1)-Метод. Положим Ь = 1, М = 0, Я = 1 в (3.1), получим одношаговый метод 5-го порядка погрешности
где
У (sh) = уг + Fo «1 (s) + Fo1) U2 (s) + Ф« (s) + Ф1 u4( s) + Ф11) U5 (s), У (h) = уг + Fo^1 (s) + Fo1) W2 (s) + Ф^( s) + Ф^4 (s) + ф11) W5 (s),
/4 s(2o-5s-6s2 + 3s3) s2 (1o - 1o s + 3s2)
«1(s) =-3o-, «2(s) =-6o-,
/4 s( 1o - 15s + 6 s2) s4 (- 1o + 12 s -3s2) s4 (5-3 s)
«3( s ) = "-2--, «4( s ) = —*-2-«5 ( s ) = --—-
3o( 1-s)2 3o( 1-s)2 6o( 1-s)
15s2 - 2s - 1 , ч 5s -2
(3.4)
(3.5)
W1 (s ) = -2-, W2 (s ) = ——, W3 (s ) =
3os2 6os ' 3os2(1- s)2'
12-28s +15s2 5s -3
W4 (s) = --—, w5( s) =
30( 1-* )2 ' 60 (1-* )•
Подставив в (3.5) значение * = 1/2, получим метод из [6].
Теорема 2. ЬМЯ(1,0,1)-Метод А-устойчив при * > 1/2 и Ь(5)-устойчив с параметром 5 = (1 - *)/* при * е (0.5, 1).
Доказательство. Обратимся к уравнению (2.6). Как и ранее, представим комплексное число ц = ХИ в тригонометрической форме, полагаем с = ео8(ф), с < 0, ё = 2* - 1 и из (3.4), (3.5) получаем (р - модуль комплексного числа Х, р = Ир)
г( И) = ¿0 4(Ц), \д (Ц)1 = \4 (р)1 = 7[ Яе (д)] 2 + [ 1т( 4 )]2 = 7ОД / Н (р), (3.6)
где
G(p) = X ' Я(Р) = X Р'Р''
go = -144oo, po = go, g1 = 288oc( 5-d), px = gx-288ooc,
g2 = 144 (4o С + d2-2oc2d + 5), p2 = g2 + 576o c2 d,
g3 = 96oc3(1 - d) + 144c(5 - 2d + d2), p3 = g3 - 192oc3 - 144oc - 288cd2,
g4 = 12(2oc3 + 4d2c2 - 24dc2 + d2 + 2), p4 = g4 + 576dc2, .
g5 = 12 c (2 - 3 d + d2), p5 = g5-48c - 24 d2c, g6 = (1-d)2, p6 = (1+ d)2. Разность между знаменателем и числителем в (3.6) имеет вид
H- G = - 24cр[ 12oo + 4р2(15 + 2oc2 + 3d2) + р4(2 + d2)] + 4dр2[ 144c2(p2 + 1o) + р4]. Разность положительна при достаточном условии d > 0 (c < 0 по условию задачи (2.6)), и так как G > 0, H > 0, то |q( р )| < 1 при s > 1/2 для любого р , откуда следует справедливость первого утверждения теоремы.
Из приведенных соотношений следует lim |q(р)| = /— = ^—d = 11—s-, что доказывает вто-
р^^ л/ p6 1+ d s
рое утверждение теоремы.
ЬМЯ(1,1,1)-Метод. Положим L = 1, M = 1, R = 1 в (3.1), получим одношаговый метод 6-го порядка погрешности с вычислительной схемой
У ( sh) = У,. + Фo «1 ( s ) + Ф01) U2( s ) + ФА ( s ) + Ф^ «4( s ) + Ф1«5 ( s ) + Ф11) «6( s ), У ( h ) = У, + Фo W1 (s) + Ф01) W2( s) + Ф^3 (s) + Ф^ W4 (s) + Ф1 W5( s) + ф11) W6( s),
(3.7)
6
6
i
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.