CC BY
58
Сер. 10. 2010. Вып. 3
ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА
УДК 519.62:517.93:521.1 Л. К. Бабаджанянц МЕТОД РЯДОВ ТЕЙЛОРА
Введение. Предлагаемая работа продолжает статью [1], в которой для полиномиальной задачи Коши были рассмотрены метод рядов Тейлора, составной метод рядов Тейлора и метод Рунге-Кутта. Важность изучения полиномиальных систем обыкновенных дифференциальных уравнений (ОДУ) еще раз подтверждена результатами работы [2], в которой введен в рассмотрение исключительно широкий класс дифференциальных уравнений, сводимых к полиномиальной форме в автоматизируемом режиме. В настоящей статье рассматривается только метод интегрирования рядами Тейлора вдоль вещественной прямой или на комплексной плоскости. Вводится понятие схемы, обеспечивающей экономное представление задачи Коши и быстрое вычисление коэффициентов разложения решений по точным рекуррентным формулам. Всесторонне исследуется вопрос об оценке погрешностей (в автоматическом режиме). В качестве примера приводится задача многих тел. Численные эксперименты не проводятся, они в достаточном объеме представлены в [1]. Основная цель, на которую ориентирована предлагаемая статья (как и [1, 2]), - построение необходимых базовых средств, обеспечивающих возможность создания программ автоматизированного решения дифференциальных уравнений методами рядов по различным системам функций.
1. Задача Коши. Рассмотрим задачу Коши
где / = (/1,...,/п); х = (х1,...,хп); хо = (х1го,...,хп,о) € Сп; а = (аи...,аш) € Сш; Мо е С.
Полиномиальной называют систему (1) (и задачу (1), (2)), в которой все функции /г - алгебраические полиномы по х1 ,...,хп с постоянными коэффициентами, зависящими от параметра а. Естественно, что для таких систем метод рядов оказывается сравнительно простым и эффективным средством нахождения приближенных решений. Реализация того или иного метода рядов для полиномиальных уравнений зависит, помимо прочего, и от формы представления этих уравнений. В настоящем разделе предложены три варианта записи полиномиальной системы, которые далее будут использованы. Полиномиальная система в форме (3) рассмотрена в п. 3 для получения явных априорных оценок погрешности в методе рядов Тейлора, остальные две
Бабаджанянц Левой Константинович — доктор физико-математических наук, профессор кафедры механики управляемого движения факультета прикладной математики—процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 84. Научные направления: дифференциальные уравнения, классическая и небесная механика. E-mail: [email protected].
© Л. К. Бабаджанянц, 2010
dx/dt = f (x, t, a) x(to) = xo (a),
(1) (2)
формы - в п. 2 для вывода различных рекуррентных формул расчета коэффициентов рядов, представляющих решения полиномиальной системы. Ради простоты, возможную зависимость уравнений и начальных данных от параметров явно указывать не будем, а решение задачи Коши будем обозначать х(г,го,хо) или х(г).
Первая форма представления полиномиальной задачи Коши:
Ь+1
¿х/йг = а + ^^ а[%]хг
т=1 геТ (т)
(3)
х(го) = хо, (4)
здесь % = (%1,... ,%п), %1,...%п € Z; Ь € [0 : те); I (т) = {г € Zn\il,...,in > 0, |г| = т},
|г| = %1 + ... + %п; х = (х1, ...,хп) € Сп; хг =
'; хо = (х1,о,..., хп,о) =
(х1 (го), . . .,хп(го)) € Сп; а = (аь ...,ап) € Сп; а% = (а1%], . . .,ап%) € Сп; € С. Вторая форма представления системы (3) - наиболее простая и экономная:
и
йхз/¿г = ^ а^}кхг(к), 3 € [1 : п],
к=0
(5)
хп, а хг(п+1),..., хг(и) - все различные нелинейные
где хг(0) = 1, хг(1) = х1, ..., хг(п) мономы в правых частях уравнений (3).
Третья форма - иерархическая, она представляется формулами
(1)
йхз/¿г = ^ а'8,ку1(1'к, 3 € [1 : п],
к=о
(6)
и(т + 1)
г(г,з) _ г(г+1,к)
уг =2-^1 аг,3,к уг + 1 ,
к=о
3 € [1 : т(г)], г € [1 : в — 1],
(7)
в которых использованы обозначения: уг = (уг,1,..., Уг,т(г)), г € [1 : в]; у8 = х, т. е.
( \ г(т,0) -, г(т,1) г(т,т(т))
уе,1 = х1,...,у8т(в) = хп, т(в) = п; уГ = 1, уГ = уг,1,...,ут = уг,т(г),
г-| и __ Г1 1 г(т,т(т) + 1) г(т,и(т))
г € [1 : в]; при каждом г € [1 : в] величины ут ,...,ут различные нели-
нейные мономы по уГ}1,..., ут,т(т) (в частности, при г = в это мономы по х1,..., хп).
В формулах (6), (7) величина в - любое фиксированное натуральное число. Будем называть его числом уровней третьей формы. При в = 1 третья форма совпадает со второй. При в € [2 : те) третья форма является обобщением второй. Вместе с тем простым перемножением полиномов третью форму всегда можно свести к форме с меньшим числом уровней и, в конце концов, - ко второй. Необходимость в третьей форме может возникнуть по разным причинам, но важно, что к такой форме полиномиальной системы сводятся системы ОДУ в результате применения к ним метода дополнительных переменных в самом общем случае в автоматизированном режиме [2].
2. Схемы и коэффициенты Тейлора. Рассмотрим упорядоченный набор Т = (хг(1),..., хг(п), хг(п+1) ,...,хг(м)) при условии, что 2 < 1%(п + 1)| < 1%(п + 2)1 < ... < 1%(М)| ^ Ь +1, и обсудим задачу вычисления мономов хг(п+1) ,...,хг(К) набора Т в предположении, что известны первые п его мономов, т. е. мономы хг(1) = х1,..., хг(п) = хп.
Вычислить мономы хг(п+1),... в Т возможно всегда, так как каждый из мономов можно получить, перемножив соответствующее число раз исходные величины
х
х1,...,хп, например х1х3хз = Ж1Ж1Ж2Ж2Ж2Ж3, но ясно, что это не самый эффективный способ определения Т. Вычисляя мономы последовательно слева направо в Т, можно для получения очередного монома использовать в качестве сомножителей уже рассчитанные ранее мономы, если таковые есть. Например, х1х3хз можно вычислить как (х!х3)хз или как х\(х3хз), если моном х\х3 или мономы ж3жз и ж2 ранее были получены.
Если очередной моном хг(т) в наборе можно найти как произведение ранее вычисленных мономов хг(р),хг(ч), то это означает, что г(г) = г(р) + ъ(ч), причем 1 ^ р < г, 1 ^ ц < г. Если при |г(г)| ^ 3 все мономы х1(т) можно рассчитать таким образом последовательно (мономы в Т упорядочены), то можно ввести в рассмотрение схему -упорядоченный набор Б = ((р(п +1), ц(п +1)),..., (р(М), д(М))) из N — п пар натуральных чисел р(г), ц(г), удовлетворяющих условию г > р(г), ц(г) для любого г € [п +1: N].
Набор Т может не иметь схемы, а может иметь одну или несколько. Если у набора Т есть схема, ее можно использовать при расчете Т: х1(п+1'),..., х1(к) вычисляются слева направо в этом порядке, причем очередной моном х1(т) определяется как произведение уже полученных х1(р(г)) и х1(ч(г)). Если у набора Т нет схемы, то его всегда можно дополнить до какого-то набора Т', имеющего схему. Это очевидно, так как любой набор Т мономов степени не выше Ь +1 можно дополнить до набора Т, содержащего все мономы степени не выше Ь.
Пусть Т, Т' - упорядоченные наборы мономов. Будем писать Т С Т' или Т' Э Т в том случае, когда все мономы набора Т содержатся ив Т'. Если Т' Э Т и Т' имеет схему, то будем называть Т' оболочкой для Т. Если мы хотим найти Т на основе некоторой схемы, а он ее не имеет, то следует рассчитывать не Т, а его оболочку Т', причем чем меньше дополнительных мономов в Т', тем быстрее, вообще говоря, будет процедура вычисления набора Т.
С одной стороны, для каждого данного Т нахождение оболочки Т' с небольшим, по возможности, количеством дополнительных мономов - вполне решаемая задача. С другой стороны, построение универсального алгоритма нахождения подходящей оболочки Т' для произвольного Т требует отдельного рассмотрения, и здесь данным вопросом заниматься не будем.
Дополним, если необходимо, набор мономов Т до оболочки Т' со схемой Б = ((р(п + 1), ц(п + 1)),..., (р(и), ц(и))) и обратимся к уравнениям (5) с тем, чтобы получить искомые формулы для коэффициентов Тейлора. Подставив
то
х'(к) =]Т хиА* — *о)р, к € [0: и], (8)
р=0
в равенство
х*(к) = х1(р(к))х1(<1(к)), к € [п +1:и], (9)
и в уравнения (5), находим
то р
^2[хк,р — хр(к),1хч(к),р-1](* — *0)Р = ° (Ю)
р=0 1=0
то и
У^[(р + 1)х^^р+1 — аз,кхк,р](* — *о)р = 0. (11)
р=0 к=0
Приравнивая нулю коэффициенты при (г — ¿о)р в равенствах (10), (11) и вспоминая начальные условия (4), получаем расчетные формулы
Хк,0 = хк(¿0), к е [1 : п], (12)
р
хк,р = 53 хр(к),1хч(к),р-1, к е [п +1 : и]
1=0
}, Р = 0,1,.... (13)
хк
,р+1 = (р + 1) ак,1Х1,р, к е [1 : п]
1=0
Формулы (12), (13) позволяют вычислять последовательно все коэффициенты хк<р разложений (8) величин Х1,..., хп, х1(п+1),..., х1(и) в ряд Тейлора. Для того чтобы использовать эти формулы, достаточно иметь в своем распоряжении начальные данные хк (¿о) при к е [1 : п], коэффициенты ак,1 при к е [1 : п], I е [0 : и] и соответствующую схему $ = ((р(п + 1), д(п +1)),..., (р(и), д(и))). Важно, что теми же величинами задается и сама полиномиальная задача Коши, и такой способ задания может быть предпочтительным не только при ее пошаговом интегрировании (методом рядов Тейлора, методом Рунге-Кутта и т. п.), но и при решении других задач, связанных с дифференциальными уравнениями.
Аналогичным способом выводим теперь формулы для коэффициентов Тейлора на основе представления уравнений (3) в форме (6), (7). При любом г е [1 : в] схему для набора мономов (уГ(г'т(г')+1,..., уГ(г'и(г^) обозначим символом
= ((р(г, т(г) + 1), ц(г, т(г) + 1)), ..., (р(г, и(г)), ц(г, и(г)))).
Подставляя
то
уЧг,к) =£угЛр(г — 1о)р, к е [0: и(г)], г е [1 : в], (14)
р=0
в равенства
у?г'к) = угг(г'р(г'к))угг(г'^(г'к)), к е [т(г) + 1:и(г)], г е [1 : в], и в равенства (7), (9), находим
то р
"^2\уг,к,р — ^ уг,р(г,к),1утл(т,к),р-1}(г — г0 )р = 0, (15)
р=0 1=0 то и(т+1)
У^\ут,],р — а,т,],кут+1,к,р](г — ¿0)р = 0, (16)
р=0 к=0 ТО и(1)
У^[(р + 1)хз,р+1 — у1,к,р](г — ¿0)р = 0. (17)
р=0 к=0
Приравнивая нулю коэффициенты при (г — ¿0)р в равенствах (15)-(17) и вспоминая начальные условия (4), получаем расчетные формулы
ув,к,0 = хк,0 = хк(¿0), к е [1 : п], (18)
Ут,к,\
Р
Е
1=0
Ут,р(т,к),1Ут,д(т,к),р-1, к е [т(г) + 1 : и(г)],
г(т)
}, р = 0,1,... . (19)
Ут-1,к,р = ^2 аг,к,1 Уг,1,р, к е [1 : т(г — 1)],
1=0
г = в, в — 1,. ..,2,
и(1)
Уз,к,р+1 = Хк,р+1 = (р + 1)-1 ^ У1,1,Р, к е [1 : п]
1=0
Формулы (18), (19) позволяют вычислять последовательно все коэффициенты Ут,к,Р разложений (14) величин У1т(т'к в ряд Тейлора по степеням г — ¿о, в том числе и искомые коэффициенты хк,р = УВ}к}Р при к е [1 : п], р е [1 : те), т. е. коэффициенты разложения решения х = (х1,..., хп) задачи (6), (7), (4) в ряд по степеням г — г0.
Для того чтобы использовать формулы (18), (19), достаточно иметь в своем распоряжении начальные данные хк(¿0) при к е [1 : п], коэффициенты ат,к,1 при г е [1 : в], к е [1 : п], I е [0 : и(г)] и схемы Б1,...,Б(см. выше, перед (14)). Теми же величинами задается и сама полиномиальная задача Коши.
3. Оценки. Рассмотрим полиномиальную задачу Коши в форме (3), (4). Ее решение обозначим х(Ь,Ь0,х0) или х(г). Кроме того, будем использовать обозначения
Лк)
Зкх/дьк, х0к) = х(к) (г0), \х\ = тах \хн\, ор(г0) = [г е С\\г — г0\ <р},
ге[1:п]
м
(то) — ¿оУ
Тмх(г,г0,х0) = жд ■—5Тмх(г,го,х0) = х(г,г0,х0) -Тмх(г,г0,х0), (20) ^—' т!
т=0
где Тм и ЗТм - операторы, которые решению х(г, г0,х0) сопоставляют полином Тейлора Тмх(г, г0, х0) и остаточный член ЗТмх(г, г0, х0) соответственно. Радиус сходимости ряда Тейлора Тжх(г,г0,х0) обозначим Е(г0,х0).
Метод рядов Тейлора решения задачи Коши (3), (4) заключается в построении таблицы приближенных значений хы = х(гы) по формулам
х1 = Тмгх(г1,г0,х0), ..., хт = Тм„х(гт,гш-1,хш-{),...,
(21)
где М1, М2,... - натуральные числа, г1 = г0 + , г2 = г1 + Н2,..., а Н1, Н2,... удовлетворяют неравенствам \Нт\ < К(гт-1,хт-1).
Вычисление каждого значения х(гы) называют шагом метода, а число - величиной этого шага. В общем случае интегрирования вдоль кривой на комплексной плоскости - комплексные числа, и точки г^ лежат на данной кривой. При интегрировании вдоль вещественной прямой величины гы, вещественны.
Для вычисления х(г,г0,х0) при некотором заданном г е С с высокой точностью по формулам (21), даже для г из круга сходимости Од((о,Хо)(г0) = [г е С\\г — г0\ < Е(г0,х0)}, число шагов г может оказаться большим, что может стать причиной ускоренного накопления ошибок округления и увеличения времени счета. Потому на каждом шаге целесообразно использовать по возможности больший шаг. Этого можно добиться, если иметь в своем распоряжении априорные гарантированные оценки величин Е(г0,х0) и ЗТмх(г,г0,х0) при г е Ощ(оХо)(г0). Получению таких оценок посвящен п. 3.
3.1. Оценки для линейных систем. Если в задаче (3), (4) Ь = 0 (т. е. уравнения линейны), то формально можно перейти к нелинейному случаю Ь = 1 (т. е. к квадратичным уравнениям), если ввести дополнительную переменную хп+1 и задачу Коши для нее: хп+1 = ех^п+1, хп+1(го) = 0. И все же мы рассматриваем оценки для линейных уравнений отдельно, так как они лучше, чем в общем случае, выводятся проще и, кроме того, в п. 3.3 оценки для нелинейных уравнений будут получены, как их естественное обобщение.
3.1.1. Основные оценки для линейных систем. Рассмотрим задачу Коши
йх/йг = а + Ах, (22)
х(г0) = хо, (23)
где х = (хь . .. ,хп) £ Сп; хо = (х^о, ..., хГ1,о) £ Сп; а = (ах,. . .,ат) £ Сп; А = (а^); I, ^, а^^ £ С.
Кроме (20), будем использовать также обозначения
(Ах)г =У2 а*>ЗхЗ, Р =1/s, в = \\А\\ю = ^ ^ = ЕКу ^ (24) 1=1 1=1
м Тт
Тмет = У] —, 5Тмет = ет - Тмет. (25)
т!
т=0
Предложение 1. При любом г € С решение х(г,г0,х0) задачи (22), (23) удовлетворяет неравенству
\5Тых(Мо,хо)| < (|хо| + |а| р)5Тме11-101/р. (26)
В дополнение к (26) отметим, что из (25) следует неравенство
5Тыеь-1о11р < eJt-t0|/p(\t - £о| /р)м+1/(М + 1)! (27)
Доказательство. Так как х(1) = А1х + А1-1 а, |(А1 х)^ | ^ х /р1 при г £ [1 : п], I £ [0 : ж), то
^Тм х(Ь,Ьо,хо)
1=М +1
1 и
<
^(|х0| + Нр) £ ^ ,Р = (|х0\ + \а\р)5Тме^/р. (28)
1=М+1
Это и требовалось доказать.
Чем меньше в, тем лучше оценки (26), (27). В задаче (22), (23) введем замену:
ху = а уу, ] £ [1 : п], а1,...,ап > 0, (29)
чтобы получить возможность выбором параметров а = (а1,..., ап) уменьшить в. Получаем задачу Коши
йу/йг = Ь + Ву, (30)
y(to)= yo, (31)
где y = (yi, ■ ■ ■ ,yn); yo = (yio, ■ ■ ■,yn,o); b = (bi,. .. , bn); B = (bij); yi = а-1хн; yi,o =
а- Xi,o; bi = а- ai; bi,j = а- a.jaiyj.
Применяя предложение 1 к задаче (30), (31), приходим к следующему результату.
Следствие 3.1. Если используются обозначение (25) и (вместо (24)) обозначения
n
р = р(а) = 1/s, s = s(a) = max s-, si = Si (а) = а-1 a.j \aij \ ,
j=1
\yo\ = шах(а-1 \xiio\), \b\ = max а-1\а\ (32)
i^[1:n] i^[1:n]
то при любых t G C, i G [1 : n] компоненты xi решения задачи (22), (23) удовлетворяют неравенствам (см. (28)-(31))
\STmxH(t,to,xo)\ < аМ + \b\р)6ТмeJt-tol/p. (33)
Возможность выбора масштабирующих множителей аj для уменьшения величины s = s(a) делает следствие 3.1 реальным средством автоматического назначения шага интегрирования с априорной гарантированной оценкой локальной погрешности.
Как выбрать масштабирующие множители а.1,...,а.п в (33), чтобы уменьшить величину s^)? Ввиду важности этого вопроса, рассмотрим два различных подхода к нему в п. 3.1.2, 3.1.3.
Замечание 3.1. При любом подходе к выбору масштабирующих множителей стоит иметь в виду, что в практических вычислениях достаточно использовать их грубые приближения с относительной погрешностью порядка примерно 10%, так как следствие 3.1 остается истинным при любых положительных масштабирующих множителях, с одной стороны, и, с другой - небольшое их изменение (во второй значащей цифре) не приведет к заметному ухудшению оценки (33).
3.1.2. Масштабирующие множители на основе теоремы Перрона. Существуют ли величины а* = (а*,..., аП), s* = s^*) = min s^) и как их найти?
ai ,...,ап >o
Чтобы ответить на этот вопрос, рассмотрим в удобной для нас форме следующий известный результат Перрона [3].
Теорема Перрона. Пусть n х п-матрица P = (pi}j) положительна, т. е. pi}j > 0 при всех i,j G [1 : п]. Тогда истинны следующие утверждения:
существует единственное собственное число X(P) этой матрицы с наибольшей абсолютной величиной;
это собственное число положительное и простое, а соответствующий ему собственный вектор может быть выбран положительным;
имеет место равенство A(P)= min maW Y^n=-\ Pi jxj/xA .
xi,...,Xn>o ie[1:n] V j J
Замечание 3.2. Более общую теорему Фробениуса и другие результаты о собственных числах и собственных векторах неотрицательных матриц можно найти в [4]. Здесь не рассматривается эта теорема, так как за счет незначительного ухудшения оценки (33) можно считать матрицу A+ = (аj\) положительной.
Таким образом, если \ai j \ > 0 при всех i,j G [1 : п], то в качестве масштабирующих множителей а.1,...,а.п в следствии 3.1 естественно использовать компоненты положительного собственного вектора а* = (а*,..., аПП) матрицы A+ = (а j\), соответствующего ее максимальному по абсолютной величине собственному значению A(A+).
3.1.3. Выбор масштабирующих множителей по коэффициентам Тейлора. Начнем с простого, но важного замечания: для системы (22) линейных дифференциальных уравнений с постоянными коэффициентами вычислять масштабирующие множители потребуется всего один раз, так как их выбор не зависит от начальных данных. Тем не менее вычисление, даже с небольшой точностью, собственного значения А(А+) и соответствующего ему собственного вектора а* = (а*, . ..,а*п) при больших п может оказаться слишком трудоемким. Поэтому мы предлагаем здесь еще один способ построения масштабирующих множителей, основанный на предварительном нахождении коэффициентов = х(т)(г0)/т! решения х(г) = (х1(г),...,хп(г)) задачи (22), (23). Стоит отметить, что необходимость предварительного нахождения этих величин не увеличивает трудоемкости метода рядов Тейлора, так как они все равно должны определяться на каждом шаге.
Предлагаемый эвристический способ проводится по следующей схеме:
1) начальное приближение для р рассчитывается по формуле (24);
2) вычисление а = (а1,..., ап). Обратимся к неравенствам (33) при
М = 0, ¿-¿о = рг, г = те^Т1р,
где р £ [0,2п), а т - некоторое фиксированное положительное число (например, т = 0, 5). Положим
/ N N
£хмХ1(} (ет - 1)-1 ,
1=1 )
где N - количество известных (вычисленных) коэффициентов Тейлора. Из неравенств (33) получаем
^ < аМ + ( Щ). (34)
Величины а = (а1,..., ап) будем искать из системы равенств
= аI(|уо| + ( |Ь|) (35)
ai = max
в предположении, что все ai отличны от нуля. Тогда можно положить
ai = ai ai/ai, i £ [1 : n], i = I, (36)
где I - любое фиксированное значение индекса из [1 : n]. Правые части в (35) (как и в (32)-(34)) не изменятся, если все ai,.. . ,an умножить на одно и то же положительное число, поэтому в (36) за ai можно принять любое положительное число. Так как вариантов формул (36) всего n (индекс I может принимать значения 1,...,n), то следует выбрать тот из них, который дает наибольшее значение p(a) (см. (32)).
Замечание 3.3. Величины ai можно вычислять достаточно грубо, однако на практике (в программе) формулу (36) следует видоизменить подходящим образом, чтобы избежать деления на нуль или чего-либо подобного. Мы не будем обсуждать здесь этот вопрос - каждый программист обычно решает такие проблемы по-своему. Отметим только, что в критических случаях всегда остается возможность использовать теорему Перрона (см. п. 3.1.2).
3.2. Введение в метод бесконечных систем. Применим метод бесконечных систем для доказательства неравенств (см. предложение 3), при помощи которых в п. 3.3 получим оценки для нелинейной задачи (3), (4), обобщающие оценки п. 3.1
для линейной задачи. Опуская детали, идею метода бесконечных систем можно описать так:
Шаг 1. Задачу Коши (3), (4) сводят к задаче Коши для бесконечной системы линейных обыкновенных дифференциальных уравнений первого порядка, разрешенных относительно производных.
Шаг 2. Используя линейность, получают результаты для последней задачи.
Шаг 3. Эти результаты интерпретируют в терминах задачи Коши (3), (4).
Рассмотрим задачу Коши (3), (4) при а = 0, т. е. рассмотрим задачу
Ь+1
<Х/& = ^ Е аЩх*, (37)
т=1 геТ (т)
x(to) = Хо, (38)
где x, x0, a[i] G Cn; t,t0 G C.
Для любого i G UmmZi I(m) введем новые переменные:
x[i] = xi. (39)
Множества Хт = {x[i] | i G I(m)} упорядочим по возрастанию m, а элементы в Хт -так, чтобы из ii = ki,..., ij = kj, ij+i > kj+i следовало, что x[i] предшествует x[k]. Упорядоченные переменные x[i] обозначим символами zi,z2,..., а вектор из них - z = (zi,z2,...). Кроме того, введем в рассмотрение начальный вектор Z0 = (zi,0,z2,0,...,zn,0,,zn+i,0,...) = (xi,o,..., xn,o, xl,0, xi,0x2,0,..., x'n,o,...), соответствующий начальному вектору x0 задачи (37), (38).
Номера p компонент zp вектора z связаны с мультииндексами i G |Jmm==i I(m) переменных x[i] некоторым взаимно-однозначным соответствием
p = w(i), i = Q(p)(Q = ^-i), (40)
которое обладает очевидным свойством:
^(Xq) = (aq-i(n) : aq(n)], (41)
где aq(n) - количество элементов множества |Jqm=i Хт, причем
aq(n) = (q + n)!/(q!n!) - 1. (42)
Предложение 2. Пусть a[i] доопределены так, что a[i] = 0, если |i| > L + 1 или хотя бы одна из компонент i отрицательна. Тогда функции x[k] удовлетворяют уравнениям
L
dx[k]/dt = ^ £ a[k, i] x[k + i], (43)
m=0 ieJ (m)
здесь
n
J (m) = {i G Zn Iii > -1,..., in > -1; Iii = m}, a[k,i] = ^ kj aj [ii, ...,ij + 1, .. .,in].
j=i
Доказательство. Действительно:
n L+i n L+i
dx[k]/dt = E kj (xk/xj) £ £ aj [i]xi = E £ £ kjaj [i]xk+i/xj =
j=i m=i iel(m) j=i m=i ieI(m)
= Е кз аз [г1,...,гз + 1,...,гп]х[к + ¡] =
у = 1 т=0 i£J(т)
= ^ кз аз [г1,...,гу + 1,...,гп}\ х\к + ¡}.
т=0 ieJ(т) \3=1 /
Что и требовалось.
Запишем систему (43) и начальные условия как задачу Коши для функции г:
дг/дг = Ог, (44)
г(1о) = го, (45)
где г = (¿1, ¿2,...), а О - комплексная бесконечная матрица. Каждая строка матрицы О и каждый ее столбец содержат лишь конечное число ненулевых элементов. Поэтому, в частности, определены все целые неотрицательные степени матрицы О. Дифференцируя равенство (44), получаем
¿¿/Л = Ог, д2 ¿/¿г2 = Одг/дг = ООг = О2г,...,дг г/Аг = Ог г. (46)
Замечание 3.4. Переход от нелинейной задачи (37), (38) к линейной задаче (44), (45) Р. Беллман [5] называет линеаризацией Карлемана. Общие формулы этой линеаризации (43) предложены в работе [6] (см. также [1, 7, 8]). Они применимы и в случае, когда ау [¡] не постоянные, а комплекснозначные функции времени г и каких-то параметров, причем для Ь допускается и значение Правда, в последнем случае матрица О не будет матрицей с конечными строками, и тогда правые части в (44) будут не конечными суммами, а бесконечными рядами.
Теперь используем формулы (43) для получения оценок величин (46) (см. [6]). Предложение 3. Пусть а[г] = (а1[г],..., ап[г}), х0 = (х1,0,..., хп,0) - коэффициенты и начальные данные задачи (37), (38), О - матрица уравнения (44), (Ог- ¡-тая координата вектора Ог и используются обозначения
Ь+1
ву, = ^21т-1 ^2 а [г]1
зе[Ъп] т=1
= Т. Чт Е а [¡1,...,гу + 1,...,гп]\), 7 = Ы . (47)
т=0 ieJ(m)
Тогда для любых натуральных р, д, связанных отношением
р е (а— (п): ач (п)] (48)
(см. (40)-(42)), и для любого натурального ] выполняются неравенства
з-1
|(О3г)р| < в37' Ц (д + тЬ). (49)
т=0
Доказательство. Из определения (39) и 7 (см. (47)) выводим неравенство
\х[к + ¡]| = \х[к]\\хЩ\ < 77№. (50)
в
Если p, q связаны отношением (48), то это означает, что zp = x[k] при |k| = q, а тогда, используя (43), (44), (50), имеем следующие неравенства:
Ln
|(Gz)p = ^ = |x[k]| <Е Е J2kj a[ii,...,ij + 1,...,in]||x[k + i]| <
m=0 ieJ (m) j=i
nL
<yqY.kEYm Y, a[n,...,ij + 1,...,in]|<«7qq
j=i m=0 ieJ (m)
и, значит, неравенство (49) доказано при j = 1.
Введем обозначения yp = drzp/dtr, y[k] = drx[k]/dtr, y = (yi,y2,...), из которых следует, что yp = y[k], y = drz/dtr и y = Grz (см. (44)). Дифференцируя (43) и (44) r раз, получаем аналогичные уравнения относительно y[k] и y:
L
dy[k]/dt =£ ]Г a[k,i]y[k + i], (51)
т=0 ieJ (т)
dy/dt = Gy. (52)
Предположим теперь, что неравенство (49) истинно при j = r. Вспоминая определение функции ш (см. (40)), выводим неравенства
r-i
b[k + i]| = \y^(k+i)\ = |(Gr z)u{k+i) \ < sr7lk+il П (|k + ^ + mL). (53)
m=0
Используя (51)-(53), получаем
\(Gr+iz)p\ = |(Gy)p| = ^ = |y[k]| <
Ln
^E E J2kj a [ib . ..,ij + 1,...,in]||y[k + i]| <
m=0 ieJ (m) j=i
L n r-i
E Y.kja[n,...,ij + 1,...,in]|«r7'k+i'(|k+i| + vL)<
m=0 ieJ (m) j=i v=0
r-i r
< sr+i7qq П(q + L + vL) = sr+y П(q + vL),
v=0 v=0
т. е. неравенство (49) выполнено при j = r +1. А это и требовалось доказать. 3.3. Оценки для нелинейных систем.
3.3.1. Основные оценки для нелинейных систем. Рассмотрим задачу (37), (38) (см. также задачу (3), (4)) и будем использовать обозначения (20) и
m- i
Ь(т) = (1 - т)-i/L, Ь^ = (dmЬ(т)/dTm)т=0 = П (1/L + l),
l=0
M
TmЬ(т) = E Ь^тm/m!, STmЬ(т) = Ь(т) - TmЬ(т). (54)
0
Предложение 4. Истинны следующие утверждения:
а) решение х задачи (37), (38) регулярно в круге Ор(г0), где
Ь+1
Р = 1/(Ьв) в = тах , = V \х0\т-1 У] \ау[г]\;
1 ] т=1 (т)
б) при г е Ор(г0) имеет место неравенство
\5Тмх(Ь,Ь0,х0)\ < \х0\ 5ТмЬ(\г — ^0\ /р). Замечание 3.5. При \г — г0\ /р < 1 из (54) следует неравенство
бТмЬ(\г — ¿0\ /Р) < Ь(\г — ^0\ /Р) \(г — и)/р\м+1.
Доказательство. Набор из первых п компонент решения задачи (44), (45) совпадает с решением задачи (37), (38), поэтому из (46) и (49) (при д = 1, р е [1 : п]) получаем
1-1
|x(l)i < sl 7 П (1+mL)■
m=0
Используя это неравенство, находим
ISTm x(t,t0,x0)l
Е x(0l)(t - t0)1 /l!
x{l)(t - U)l /l!
l=M+1
= |xo| STmb(lt - tol /p)■
l-1
\i
< ixol (it - tol /Р)1Ц (1/L + m)/l!
l=M+1 m=0
Что и требовалось доказать.
В (3), (4) введем дополнительную переменную xn+i = 1, заменим свободный член a на axn+i и используем замену (см. (29) и далее) xj = a.jyj, j G [1 : n], где a.j - произвольные положительные параметры. Применяя к полученной задаче предложение 4, получаем следующий результат.
Следствие 3.2. Пусть начальные данные xj,0 в (4) и положительные параметры а = (ai, ■ ■ ■ ,an) (масштабирующие множители) удовлетворяют неравенствам
0 < lxj,ol < a.j,j G [1 : n]■ (55)
Тогда истинны следующие утверждения:
а) решение x(t,t0,x0) задачи (3), (4) регулярно в круге Op(t0), где
p = p(a) = 1/(Ls(a)), s = s(a) = max Sj(a),
je[1:n]
L+1
sj = sj(a) =a-1 • (laj1 + E E a lajw|);
m=1 \i\=m
б) при t G Op(t0), j G [1 : n] имеют место неравенства
ISTm xj (t,to,xo)i < a.j STm b(\t - to | /p) ■ (56)
Замечание 3.6. Условие 0 < \х^,о\, 3 € [1 : п], в (55) можно считать выполненным за счет незначительного уменьшения р и ухудшения оценок (56).
Возможность выбора масштабирующих множителей а^ для уменьшения величины в = в (а) делает следствие 3.2 (как и следствие 3.1 в линейном случае - см. п. 3.1.1) реальным средством автоматического назначения шага интегрирования с априорной гарантированной оценкой локальной погрешности. В линейном случае (см. п. 3.1.2) для выбора а^ мы воспользовались теоремой Перрона, что позволило найти наилучшие возможные значения этих множителей. Кроме того, в п. 3.1.3 был предложен практический прием выбора а^ на основе известных коэффициентов Тейлора. В нелинейном случае вопрос выбора масштабирующих множителей осложняется в основном следующими обстоятельствами:
1) так как масштабирующие множители зависят от начальных данных (см. (55)), то их придется выбирать на каждом шаге метода рядов Тейлора;
2) в нелинейном случае мы не располагаем аналогом теоремы Перрона, а если бы его имели, то использование такой теоремы на каждом шаге слишком трудоемко.
Поэтому воспользуемся практическими приемами: в п. 3.3.2, как и в линейном случае, рассмотрим выбор а^ на основе коэффициентов Тейлора, а в п. 3.3.3 применим линеаризацию задачи (3), (4) в окрестности начальных данных хо.
Замечание 3.7. Как и в замечании 3.1, отметим, что при любом подходе к выбору масштабирующих множителей стоит иметь в виду, что в практических вычислениях достаточно использовать их грубые приближения с относительной погрешностью порядка 10%, так как следствие 3.2, с одной стороны, остается истинным при любых положительных а^ и, с другой стороны, небольшое их изменение (во второй значащей цифре) не приведет к заметному ухудшению оценок (56).
3.3.2. Выбор масштабирующих множителей по коэффициентам Тейлора. Предложим здесь способ построения масштабирующих множителей, основанный на предварительном нахождении коэффициентов Тейлора хгт = х(т(Ьо)/т! решения х(£) = (х1(£),...,хп(£)) задачи (3), (4). Он повторяет, с необходимыми изменениями, аналогичный способ, рассмотренный в п. 3.1.3 для линейных уравнений. Необходимость предварительного нахождения коэффициентов Тейлора не увеличивает трудоемкости метода рядов Тейлора, так как они все равно должны определяться на каждом шаге.
Предлагаемый эвристический способ проводится по следующей схеме:
1. Начальное приближение для р. На первом шаге начальное приближение для р вычисляется по линейному приближению (см. п. 3.3.3). На следующих шагах в качестве начального значения р принимается значение этой величины на предыдущем шаге.
2. Вычисление а = (а1,..., ап). Обратимся к неравенствам (56) при М = 0 и £ — = рг, г = те^^-^, где ср € [0, 2тг), а г - некоторое фиксированное положительное число (например, т = 0.5). Положим
где N - количество известных (вычисленных) коэффициентов Тейлора. Из неравенств (56) получаем а^ ^ а^ и полагаем а^ = шах{<г^, \х^,о\} в предположении, что ху,о = 0 (см. (55) и замечание 3.6 после следствия 3.2).
3.3.3. Выбор масштабирующих множителей по линейному приближению. Введем в задаче (3), (4) замену £ = (£1,...,£п) = х — х0, и правые части
полученных из (3) уравнений разложим по степеням £ь...,£п. Удерживая только линейные слагаемые, выписываем линейную задачу Коши
(1 п
-¿^з = сз + Е сзА»1 €з(*о) = °> 3 е [1 : п],
^=1
где
L+1 L+1
Cj = a3 +Е Е aj C3V = E E aj /xv,0-
m=1 ieI(m) m=1 iel (m)
К этой задаче применим способы, рассмотренные в п. 3.1.2 и 3.1.3, с тем, чтобы найденные таким образом масштабирующие множители использовать для задачи (3), (4).
3-4- Автоматический выбор шага. Для того чтобы обеспечить в методе рядов Тейлора автоматический априорный выбор величины шага с гарантированной заданной точностью решения, мы предложили следующие средства для линейных и нелинейных уравнений:
а) для линейных уравнений:
• следствие 3.1 (см. п. 3.1.1 после предложения 1),
• алгоритм нахождения масштабирующих множителей при помощи теоремы Перрона (см. п. 3.1.2),
• алгоритм нахождения масштабирующих множителей при помощи коэффициентов Тейлора (см. п. 3.1.3);
б) для нелинейных уравнений:
• следствие 3.2 (см. п. 3.3.1 после предложения 4),
• алгоритм нахождения масштабирующих множителей при помощи коэффициентов Тейлора (см. п. 3.3.2),
• алгоритм нахождения масштабирующих множителей на основе линейного приближения (см. п. 3.1.3).
Вместо следствий 3.1 и 3.2 удобнее использовать вытекающие из них следствия, которые мы сейчас сформулируем. При любом фиксированном M e [1 : введем
в рассмотрение функцию и аргумента т > 0 такую, что
и(т) = STm eT, (57)
а символом и-1 обозначим обратную ей функцию.
Следствие 3.1.1. Пусть используются обозначения из следствия 3.1 и формула (57). Пусть Х1,...,Хп, e1,...,en - положительные числа и e = (|у0| + р Щ)-1 min (eiXi/ai). Тогда
ie[1:n]
(|t - toI < р ■ u-1(e)) ^ (ISTmXi(t,t0,X0)| < eiXi, г e [1 : n]). (58)
При любом фиксированном M e [1 : введем в рассмотрение функцию v аргумента т e [0,1) такую, что
v(T) = STm Ь(т), (59)
а символом v-1 обозначим обратную ей функцию.
Следствие 3.2.1. Пусть используются обозначения из следствия 3.2 и формула (59). Пусть х1,...,Хп,е1,...,еп - положительные числа и е = min (eixi/ai). Тогда
i€[1-.n]
(|t - toI < р ■ v-1(e)) ^ (ISTmXi(t,to,xo)| < eiXi, г e [1 : n]). (60)
Как можно видеть из неравенств (58), (60), для определения максимально возможного шага р • и-1(е) (в линейном случае) и р • v-1(e) (в нелинейном случае) необходимо уметь вычислять функции u-1, v-1 при заданных значениях е. Как сделать эти вычисления быстрыми - дело программиста. По нашему мнению, проще всего организовать быстрые вычисления, если задать такие функции в виде таблиц.
На этом мы заканчиваем здесь обсуждение вопроса об автоматическом выборе шага - приведенных выше средств вполне достаточно для создания соответствующих алгоритмов и программ.
4. Задача многих тел. Пусть l + 1 материальных точек Mo,...,Mi с массами то,.. .,mi соответственно движутся под действием взаимного притяжения по закону Ньютона. Пусть Oxyz - некоторая инерциальная система координат, а Xi,yi, zi - координаты точки Mi в ней. Введем в рассмотрение относительные координаты с началом в точке Мо с осями Mogi,Mog2,Mogo, которые одинаково направлены с Ox,Oy,Oz. Пусть gi,i,gi,2,gi,3 - координаты точки Mi в этой системе. Тогда gi,i = Xi - X0, gi,2 = yi - yo, gi,3 = Zi - Z0 при любом i e [1 : l]. Как известно, функции времени gj удовлетворяют уравнениям
d2gi,j/dt2 = -k2(mo + mi)gi,j/ro, i + k2 E ms[(gs,j - gi,j)/r°,i - gs,j/гоЛ (61)
sE[1:l] s=i
где rs¡i = ^Jj23j=i(gi,j —9a,j)2, i £ [1 : l], 3 = 1, 2, 3, s G [0 : /], а к - постоянная Гаусса.
Вводя при s < i, s e [0 : l], i e [1 : l], j e [1 : 3] переменные pitj = dgi}j/dt, ds¿ = r-1, получаем полиномиальную систему относительно gi,j, pi,j, di,j
dpi,j/dt = -k2(mo + mi)gi,jd%,i + k2 ^ ms[(gs ,j - gi,j)d3s¡i - gs,jd%,s],
s£[1:l], s=i
dgi,j/dt = pi,j, (62)
з
j3
Ме,1/& = (91,3 — да,3 )(рг,3 — Ра,]), г € [1 : /], 3 = 1, 2, 3, в € [0 : /], в <г.
3=1
Отметим, что порядки систем (61) и (62) равны 6/ и 6/ + /(/ + 1)/2, а количества различных нелинейных слагаемых в их правых частях равны 3/2 и 3/2 + 6/(/ — 1). Правые части уравнений (62) содержат следующие различные нелинейные мономы:
дг,3 ¿0,ъ, дг,3 ¿э ,г, дг,3 ¿э ,грг,3 , дг,3 ¿э ,гра,3 (63)
при г, в € [1 : /], г = в, 3 € [1:3]. Задачу эффективного вычисления этих мономов можно поставить в связи с решением уравнений (62) тем или иным методом, например методом Рунге-Кутта, методом рядов Тейлора и т. д. Исходный набор мономов Т здесь образуют переменные дг,3, рг,3 при г € [1 : /], 3 € [1:3] и мономы (63). Набор Т не имеет схемы, но, добавив к нему мономы ^^ г, г при в < г, в € [0 : /], г € [1 : /], получаем оболочку для Т, т. е. искомый набор мономов Т', имеющий схему. В наборе Т' группы мономов дг,3,рг,3,с!а,г; ¿2],г; г; дг,3г; дг,3гРг,3, дг,3гРа,г должны располагаться в этом порядке, а внутри групп расположение мономов может быть любым.
Все коэффициенты х^,р разложений величин х1,...,хп, хг ,...,хг в ряд Тейлора можно вычислить последовательно по формулам (12), (13). Для того чтобы можно
было получить ориентировочное представление об эффективности таких вычислений, приведем величины n, N и количество K ненулевых aj, k для полиномиальных уравнений (62), описывающих относительное движение l материальных точек:
n = 6l + l(l + 1)/2, N = n + 5l(2l - 1), K = 3l(5l +1)
Заключение. Прежде всего отметим, что как представление и использование дифференциальных уравнений реальных задач прикладной математики в полиномиальной форме, так и метод рядов Тейлора имеют давнюю и богатую историю, особенно в небесной механике (см., в частности, [1, 2, 6-25]). Основные результаты настоящей работы представлены новыми алгоритмами вычисления коэффициентов Тейлора (12), (13) и (18), (19) в п. 2 и новыми алгоритмами априорного выбора шага интегрирования с гарантированной оценкой погрешности в п. 3. Для задачи многих тел в п. 4 построена схема, позволяющая воспользоваться алгоритмом вычисления коэффициентов Тейлора по формулам (12), (13). Что касается выбора шага с гарантированной априорной оценкой погрешности в этой задаче, то ограничимся ссылкой на работы [19, 25], в которых предложены соответствующие результаты для планетной и спутниковой задач и на реальных данных показана их практическая эффективность.
Литература
1. Babadzanjanz L., Sarkissian D. Taylor series method for dynamical systems with control // J. of M. Sciencies. New York: Springer, 2006. Vol. 139, N 6. P. 7025-7046.
2. Бабаджанянц Л. К. Метод дополнительных переменных // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2009. Вып. 4. С. 3-11.
3. Беллман Р. Введение в теорию матриц / пер. с англ.; под ред. В. Б. Лидского. М.: Наука. 1969. 368 с.
4. Гантмахер Ф. Р. Теория матриц. М.: Наука, 1988. 550 с.
5. Беллман Р. Процессы регулирования с адаптацией / пер. с англ. Ю. П. Леонова; под ред. А. М. Летова. М.: Наука, 1964. 360 с.
6. Бабаджанянц Л. К. Продолжаемость и представление решений в задачах небесной механики // Труды Ин-та теор. астрономии АН СССР. 1978. Вып. XII. С. 3-45.
7. Babadzanjanz L. K. Existence of the continuations of the W-body problem // Celestial Mechanics. 1979. N 20. P. 43-57.
8. Babadzanjanz L. K. On the global solution of the W-body problem // Celestial Mechanics and Dynamical Astronomy. 1993. N 56. P. 427-449.
9. Пуанкаре А. О кривых, определяемых дифференциальными уравнениями / пер. с франц. Е. Леонтович, А. Майер. М.: ОГИЗ, 1947. 392 с.
10. Steffensen J. F. On the restricted problem of three bodies // Mat.-Fys. Medd. Danske, Videnskab. Selskab. 1956. Vol. 30, N 18. P. 75-83.
11. Steffensen J. F. On the problem of three bodies in the plane // Mat.-Fys. Medd. Danske, Videnskab. Selskab. 1957. Vol. 31, N 3. P. 98-123.
12. Мерман Г. А. О представлении общего решения задачи трех тел сходящимися рядами // Бюл. Ин-та теор. астрономии АН СССР. 1958. Т. 6, № 10. С. 713-769.
13. Rauch L. M. Iterative solution of the W-body problem for real time // ARS Journal. 1960. N 30. P. 284-286.
14. Rauch L. M., Riddel W. S. The iterative solution of the analytical W-body problem //J. Soc. Industr. and appl. math. 1960. N 8. P. 568-581.
15. Брумберг В. А. Ряды полиномов в задаче трех тел // Бюл. Ин-та теор. астрономии АН СССР. 1963. Т. 9, № 4. С. 234-256.
16. Estes R. H., Lancaster E. R. An algorithm for integrating stepwise the restricted problem in Thiele's coordinates // Celestial Mechanics. 1970. N 1. P. 297-300.
17. Griffith J. S. On a generalized Taylor scheme for numerical integration // Astronomy and Astrophysics. 1970. Vol. 8, N 2. P. 267-272.
18. Broucke R. Solution of the W-body problem with recurrent power series // Celestial Mechanics. 1971. N 4. P. 110-115.
19. Чернышева, Н. А. Метод вычисления возмущений в поле вращающегося тела // Вестн. Ленингр. ун-та. Сер. 1: Математика, механика, астрономия. 1987. Вып. 4. С. 83—89.
20. Бабаджанянц Л. К., Чекашкин Ю. В. Аналитический метод вычисления возмущений в координатах планет // Вестн. Ленингр. ун-та. Сер. 1: Математика, механика, астрономия. 1990. Вып. 3. С. 101-107.
21. Parker G. E., Sochacki J. S. Implementing the Picard Iteration // Neural, Parallel and Scientific Computation. 1996. Vol. 4. P. 97-112.
22. Parker G. E., Sochacki J. S. A Picard-McLaurin Theorem for Initial Value PDE's // Abstract and Applied Analysis. 2000. Vol. 5. P. 47-63.
23. Pruett C. D, Rudmin J. W., Lacy J. M. An Adaptive N-Body Algorithm of Optimal Order // J. of Comput. Physics. 2003. Vol. 187. P. 298-317.
24. Carothers D. C., Parker G. E., Sochacki J. S., Warne P. G. Some properties of solutions to polynomial systems of differential equations // Electron. J. Diff. Eqns. 2005. Vol. 2005, N 40. P. 1-17.
25. Babadzanjanz L. K. Error estimation for numerical solution of N-body problem // Letters to AG. 1981. Vol. 7, N 12. P. 752-755.
Статья рекомендована к печати проф. Л. А. Петросяном. Статья принята к печати 1 апреля 2010 г.
