Автоматизация решения обыкновенных дифференциальных уравнений Текст научной статьи по специальности «Математика»

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 10. 2009. Вып. 2

УДК 004.4’242 В. Н. Латыпов

АВТОМАТИЗАЦИЯ РЕШЕНИЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

1. Введение. Важной целью современной прикладной математики считают создание универсальных процедур решения задачи Коши для систем обыкновенных дифференциальных уравнений (ОДУ). Подобная процедура, разработанная в университете Техаса [1], основана на автоматическом дифференцировании [2]. В статье [3] рассматривается реализация метода итераций Пикара. В настоящей работе для автоматизации решения ОДУ используется другая идея, ориентированная на системы дифференциальных уравнений с полиномиальными правыми частями. Предлагается система генерации вычислительных процедур [4], автоматизирующая процесс построения пошаговых схем интегрирования, основанных на методе рядов Тейлора и методе малого параметра.

Возможность сведения нелинейных дифференциальных уравнений к полиномиальной системе ОДУ была известна еще Анри Пуанкаре [5]. В работах [6, 7] предложен метод дополнительных переменных, позволяющий сводить нелинейные системы ОДУ к полиномиальным для исключительно широкого класса уравнений. С точки зрения приближенных вычислений, полиномиальные системы ОДУ представляют особый интерес, так как алгоритмы их решения отличаются простотой и высокой эффективностью [3, 6, 7].

В отличие от подхода, основанного на итерациях Пикара [3], в котором решения приближаются последовательностью полиномов, рассматриваемый в статье метод рядов Тейлора позволяет получить простые априорные оценки погрешности (теоремы

В п. 2 и 3 будут рассмотрены методы приближенного решения для полиномиальных систем, основанные на представлении решений в форме рядов Тейлора и рядов по малому параметру. В п. 4 будет описано использование операторов Бм,е и Тм для автоматизированного построения пошаговых схем интегрирования.

2. Метод рядов Тейлора для полиномиальных систем. Для построения приближенного решения дифференциальных уравнений возможно применить метод рядов Тейлора. В отличие от традиционных численных методов, использующих вычисление правой части, в нем получают рекуррентные формулы для нахождения коэффициентов разложения решения в степенной ряд.

Рассмотрим автономную задачу Коши

Латыпов Виктор Николаевич — ассистент кафедры механики управляемого движения факультета прикладной математики-процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 4. Научные направления: численные методы решения дифференциальных уравнений, математическое моделирование. E-mail: latypov_viktor@inbox.ru.

© В. Н. Латыпов, 2009

1 и 3).

x(t) = X (x), X (x) = (X1,...,Xrl), x = (xi ,...,xn),

x(to) = (xi(to),..., xn(tg)) = (xio,.. .,xno) = xo,

(1)

(2)

где

N п іі Іт-1

Хк (х) = ак + } ' ^ ' • • • ^ ' ак;іі ,...,іт хіі X ■■■ X хіт, (3)

т=1 іі = 1 і2 = 1 іт=1

і = (І1, • • •, іт) Є Zm, ак;і ак;іі ,...,іт , ак Є С^

В формуле (3) коэффициенты ак ,ак;і - постоянные величины, і - т-компонентный мультииндекс. Так как Хк(х) голоморфны по совокупности аргументов х1, • ^х, в окрестности (хо), то по теореме Коши решение этой задачи может быть представлено степенными рядами

+ ^>

хк(Ь,Ьо,хо) = ^2 хк,д(г — го)4, к = 1, • • • , п, (4)

д=0

абсолютно сходящимися в окрестности точки і = го.

Метод рядов Тейлора приближенного решения задачи Коши (1), (2) состоит в замене точного решения отрезком ряда (4) в области {і : \і — іо| < К} С С:

м

ТМхк (г, г0, хо ) ^ ^ хк,4 (г іо)4 •

4=о

Целое положительное число М называется порядком метода. Тжх(г,го,хо) совпадает с точным решением исходной задачи.

Подставляя в выражение (4) г = го, легко убедиться, что постоянные члены ряда Тейлора выражаются через начальные условия

хк,0 хк (г0') •

Для вывода рекуррентных формул, определяющих коэффициенты хк,ч, подставим ряд (4) в уравнения (1) :

*(!Н =

N П %\ ^т-1 / \ / \

= ак + ’’У ' ^ ' • • • ^ ' ак;г1 ,.Лт ( х11,41 ^ ) Х '' Х ( ^ ' х£т ,4т г®т I •

т=1 £1 = 1 %2 = 1 1т = 1 \?1=0 / \?т = 0 /

Приравнивая коэффициенты при одинаковых степенях независимой переменной г, получаем искомые соотношения

ч +1

N П іі іт—1

ак 5 4 + ЕЕЕ-Еак;і Е <

т=1 іі = 1 і2 = 1 іт = 1 91 + ... + 9т = д

в которых 6д = 1, если ц = 1, и 6д = 0, если ц = 0.

Оценка погрешности. Рассмотрим полиномиальную систему с нулевыми свободными членами

Ь+1 П £1 £т-1

хз = Е Е '^•••^ аих£1 х ••• х хт• (5)

т=1 £1 = 1 £2 = 1 £т = 1

• • • • хі 4

°т ,Чт

іі ,41

К такому виду можно привести любую полиномиальную систему, добавив вспомогательную переменную xn+l = cn+l = const и уравнение Xn+l = О.

Пусть

гр /,ч Д (t - to)

lMx{t) = -----—

l=o

д d x Xl = to ’

M е N U {+to},

Oa(b) = {t е C||t - bl <a}.

Введем обозначения

L+l n ii im-l

j=1,...,n

т=1 £1 = 1 £2 = 1 £т = 1

11*11 = Н(Уь...,Уп)|| = тах |%|,

3 = 1,и

А 1

Р М1ИМН)'

Если х(г) - решение задачи Коши для системы (5) с начальными данными х(Ь0), то задача оценки методической погрешности 5 мх(г) = х(г) — Тмх(г) и радиуса сходимости ряда Тжх(г) решается следующим образом [7-9].

Теорема 1. В области Ор(г0) справедливо неравенство

||Лм1(()|К||1((о)||л_^у1/ь ^Г1

р ) V р )

Вращательное движение спутника в окрестности положения равновесия.

Рассмотрим уравнения движения искусственного спутника около центра масс на круговой орбите в гравитационном поле Земли. Введем системы координат и переменные, относительно которых будем решать уравнения движения.

Ох1У1Х1 - орбитальная система координат. Начало отсчета совпадает с центром масс спутника. Ось Ох1 направлена вдоль радиус-вектора, соединяющего центры масс Земли и спутника; ось Ох1 - по трансверсали орбиты. Кроме того, введем в рассмотрение систему Охуг, жестко скрепленную со спутником, оси которой направлены по главным осям инерции спутника. Взаимное расположение двух систем координат определяется таблицей направляющих косинусов:

X У Z

Xl а а' а'

У1 в в' в'

Z1 Y Y' Y''

Пусть р,ц,г - проекции угловой скорости ш на главные оси инерции спутника. Одним из положений равновесия является такое движение, при котором

р = шр, ц = шр', г = шр'', а = р' = = 1, а = ••• = 7' = 0^ (6)

Если ввести новые переменные 21 ••• хц по формулам

Х1 = р — шр, Х2 = ц — шр', хз = г — шр'', Х4 = ^, Х5 = 7', Х6 = ^'',

Х7 = а — 1, хв = а', хд = а'', хю = р, хц = р' — 1, Х12 = р'',

то они удовлетворяют системе следующих уравнений [10]:

х1 = —(N1 + 1)хз — К1ш0х12 + Ш1ш0х5 + х2х12 — хзх11+

+ 3NlШоX5X6 — ^(хзхи + х2х12 + ш0хцх12 + х2хзш0), х2 = 3^ш0х4 + хзхю — х1х12 + 3^ш0х4хб —

— ^(х1х12 + хзхю + ш0хюх12 + ш-1х1хз),

хз = —(N3 — 1)х1 — ^ш0хю + 3^ш0х4х5 — х2хю + х^ц —

— ^(х2хю + х1хц + ш-1х1х2 + ш0хюхц),

х4 = —ш- 1х2 + х5х12 — хбх11 + ш- 1(хзх5 — х2хб) +

+ ^ (х2 + хк+б — хк+з),

к=4

х5 = ш° 1х1 + хв + хю + ш° 1(х1хб — хзх4) + хбхю — х4х12, (8)

хб = х4 + хд + ш- (х2х4 — х1х5)+х4х11 — х5хю,

0-

х7 = —х4 — хд + ш0 1(хзхв — х2 хд) + х8х12 — хдхц,

хв = — х5 — х12 — ш-1хз + ш- 1 (х 1хд — хзхг) + хдхю — хгх12,

хд = ш° 1х2 + ш° 1( х2 х^ — х1хв) + хгхц — хвхю —

— 0•5 (х2+з + хк+б — хк),

=4

х10 = ш0 1х3 + ш0 1( х3 х11 — х2 х12 ),

х11 = ш-1( х1х12 — хзхю),

х12 = —ш0 1х1 + ш0 1( х2х10 — х1 х11) •

Правые части системы (8) - алгебраические полиномы второй степени. Решение задачи Коши возможно проводить описанным выше методом рядов Тейлора. У рассматриваемой системы существует несколько устойчивых положений равновесия [11] , в окрестности каждого из которых возможно применение метода малого параметра.

3. Метод малого параметра. В случае, когда коэффициенты системы линейно зависят от некоторого числового параметра е, для построения приближенного решения уравнений возможно применить метод разложения по этому параметру [8]. Рассмотрим задачу Коши

х(г) = Z (х,г,е), Z (х ,г,е) = ^1,^^и), х = ( хl,••• хП), (9)

х(г0) = (х1 (г0), •••,хи(г0)) = (х 10, •••,хи0) = х0, (10)

где

N и £1 £т — 1

Zk(х, г, е) = Ак + Ху Ху X] ••• Ху Ак;н...£„ х£1 х ''' х х£т , т=1 £1 = 1 £2 = 1 £т = 1

Ак;£1 ...£п ак;£1...£п + еЬк;£1 ...£п , Ак ак + еЬк^

Коэффициенты ак, ак;£1...£п, Ьк, Ьк;£1...£п - произвольные голоморфные функции переменной г.

Предполагается, что при е = 0 известно общее решение уравнений (9).

Так как Zj (х, г, е) голоморфны по совокупности аргументов хl,•••, хи, г,е в окрестности (х(г0),г0, 0), то по теореме Коши решение этой задачи может быть представлено степенными рядами

р=0

Метод малого параметра приближенного решения задачи Коши (9), (10) состоит в замене точного решения отрезком ряда (11) на промежутке [г0,т(е)]

Число М называется порядком метода. Бж,ех(г,г0,х0) совпадает с точным решением

щенной задачи. В каждом конкретном случае вопрос о радиусе сходимости (11) представляет собой трудную задачу, решение которой в данной работе не рассматривается. Для автономных полиномиальных систем после замены т = ге такие оценки могут быть получены с использованием теоремы 1 [12].

Для того чтобы решить задачу методом малого параметра, необходимо найти коэффициенты ряда по малому параметру е, представляющему решение задачи (9), (10):

р,9=0

абсолютно сходящимися в окрестности точки г = г0,е = 0. После перегруппировки слагаемых получаем выражения

(11)

д=0

в которых

м

д=0

исходной задачи. Из (11) следует, что Бм,0х(г,г0, х0) = х(0)(г,г0, х0) - решение невозму-

(12)

д=0

Согласно теореме Пуанкаре, функции х(ч)(г,г0, х0) удовлетворяют уравнениям

(13)

к=1

или в векторной форме

—= рЫ(г).

А

Здесь

а функции Р(9)(г) = (Р(9)(г), •••, рП9(г)) зависят только от выражений Z(х, г, е) и коэффициентов х (0)(г), • • • ,х(9 - 1)(г).

Для того чтобы рекуррентно определять коэффициенты х(9) из уравнений (13), необходимо добавить к ним начальные условия.

Из равенств (11) при г = г0 следует, что

х(г0, г0, х0) = х0 = ^2 х(9) (г0, г0, х0)е9 •

д=0

Так как х0 не зависит от е, то сравнивая коэффициенты при степенях е в правой и левой частях, получаем

х(0)(г0,г0, х0) = х0, х(ч)(г0,г0,х0) = 0, ц = 1,2,••• •

При каждом ц > 0 уравнения (13) - линейная неоднородная система, для решения которой достаточно найти общее решение однородной системы (система уравнений в вариациях)

^у-.1у = 0. (14)

Теорема 2 (Пуанкаре). Пусть С = (C1,•••,Cn) - произвольные постоянные и р(г, С) - общее решение (в области, содержащей точку г = г0,х = х0) невозмущенных уравнений, т. е. уравнений (9) при е = 0. Пусть С0 = (С10, •••, Си0) — значение С, при котором р(г, С) = х(г, г0, х0, 0). Тогда п вектор-функций

Ук = (уЬ-.У ) = (др(г,С)/дСк )с=С0 =

= (дп(г,С)/дС1,...,д<ри(г,С)/дСи)а=а0, к

1,...,п,

образуют фундаментальную систему решений уравнений в вариациях (14).

Для применения метода малого параметра кроме решения уравнений в вариациях необходимы также явные формулы для функций Р(9)(г).

Получим уравнения на коэффициенты ряда х^ч\г) подстановкой выражения (12) в уравнение (9):

Е

х(9)е9 =

N и £1 £т—1

1*2 ••• ^2 Ак; £1 ...£т

д=0 т=1 £1 = 1 £2 = 1

Введем обозначения

х(91)е91 | х ••• х

9=0

У х(Чт) е9т

£т

9=0

+ Ак •

£=£ £•••£• *9 = £

(91)

£1 = 1 £2 = 1 £т = 1 91+... + 9т = 9

(91 ) (9т)

Н9

Е

91 Н--'!)!■ 9

Чр^Ч,Р=1,т

хЦ х • • • х х£т

Е хГ х

.(9)

х

(0)

.7=1

(9т)

£1

Группируя слагаемые при степенях е, имеем

+ ^ + ^ \

]Г г^е* = ^ + £ I £ {Е К 1 + Ьк; £Б—} }еЛ . (15)

4=0 т=1 \ д=0 1 /

Уравнения (15) можно записать в форме (13). Для этого найдем —— и :

дгз

dZk N

Е EAfc

dzj 1 .

J m=1 *

Так как (Jz(q))k =Y^^=1 E* ak.Hlq, то P(q) при q > 1 выражается формулой

N

P(q) =^{ Eak; S ^ bfc; * Sj-i}. (16)

m=1 * *

Однородные уравнения Z(q) — Jz(q) = 0 могут быть решены по теореме 2, которая позволяет вычислить их фундаментальную матрицу Y(t).

Если при заданном q ^ 1 величина P(q) (t) известна, то решение (13) можно получить по формуле Коши

z(q)(t) = Y(t) f Y— 1 (r)P(q)(r) dr. (17)

Jo

В общем случае интеграл в правой части формулы (17) не выражается через элементарные функции. Алгоритмизация построения коэффициентов z(q)(t) возможна для отдельного класса систем.

Квазилинейные полиномиальные системы. Рассмотрим частный случай уравнений (9) - полиномиальную систему ОДУ

z = Az + eP (z),

где A - матрица n x n; P(z) = (Pi(z),..., Pn(z)) - полиномы по переменным zi,..., zn.

Пусть Y(t) - фундаментальная матрица линейных уравнений z = Az. Тогда общее решение невозмущенной задачи имеет вид z(t) = Y(t)zo. Оно совпадает с решением системы уравнений в вариациях (14).

Коэффициенты матрицы Y(t) - это полиномиально-тригонометрические и экспоненциальные функции:

a*tni exp (a*t) cos (pit) + bjtmj exp (y*t) sin (Sjt).

*j Поэтому интегралы в правой части (17) выражаются в элементарных функциях, и для коэффициентов z(q)(t) можно получать явные формулы, применяя какую-либо систему компьютерной алгебры.

Алгоритм приближенного решения задачи Коши (9), (10) состоит в вычислении приближенного решения Sm,eX(t0 + h) = M=0 z(q)(t0 + h), используя формулу (17) для выражения z(q)(t).

Оценка погрешности. Рассмотрим вспомогательную функцию

s(y) = Y-1 .m&^n max j ( E Ym E aJ';4 A E Ym E j

\ \m=1 i / \m=1 i У

Пусть

L = N — 1, a > 0,

ro = ((1 + a)Ls) 1,

Ko=[l-t^

ro

Koaro

-1/L

1C :

aro — \e(t — to) |

Для оценки методической погрешности можно воспользоваться следующей теоремой [8, 9].

Теорема 3. Справедлива следующая оценка:

М ~М+1 / ^ - *0|4 м+1

\\z(t,to,zo,e) — J2 z(l)(t)£l\\ < K||zo||

l=o

ro

Пример. Рассмотрим подробно способ применения полученного алгоритма к построению приближенного решения системы уравнений Лотки-Вольтерры [13]

У1 = У1(а - вУ2),

У2 = У2 — + 5у1).

Анализ устойчивости решений подобных уравнений проведен в книге [13]. Данная система имеет два положения равновесия - точки (0,0) и (^/д, а/в).

Произведем линейную замену переменных

Х1 = -у 1 - 1,

7

х2 = ~У2 ~ 1-

а

Для системы уравнений относительно переменных (21,22) точка (0, 0) является устойчивым положением равновесия.

Рассмотрим задачу Коши

21 = -а.22 - а.Х1Х2, (18)

22 = 721 + 121X2, ( )

2 (0) = (21(0),22(0)). (19)

Чтобы применить метод разложения по параметру, домножим квадратичные члены

в правых частях системы на величину е:

21 = -а22 + е ((-аг^)), (20)

22 = 721 + е ((72122)) . ( )

При е = 1 уравнения (20) и (18) совпадают. Общее решение невозмущенных (е = 0) линейных уравнений следующее:

еЛ1* + еЛ2* • (еЛ1*-еЛ2*

— | / 1 , V 2 ' | Г г1^ ] — УШ ( г1^

(} 1 I V -2(0) ) [ ) V -2(0)

В обозначениях, введенных выше, п = 2 и N = 2. Все величины ак;г1,г2 и Ък-г равны нулю. Выпишем коэффициенты ак-г и Ък-г1,г2:

®1 = а2 = а1 ; 1 = а2 ; 2 = 0, а1 ; 2 = -а, а2; 1 = ^,

Ъ1 ; 2,1 = -а, Ъ2; 2,1 = 1.

(?)

Выразим величины Pf,q) по формуле (16):

Pl(q) = -аБ^, P2q) = 1S(2’1), S(2,1) = £ z(l)z

(l)z (q-l)

2 z1 .

l=o

Используя формулу (17), получим выражение

+

A2 - а Y - Al

j Y9lit) + ^i92it)'\eX2t

где

с, = Ь(о)22(о), C2 = \ (az^trl(0))>

mmAi = —i^Ja 7, Л2 = 7.

/iW = зіеЗЛ2‘ + ieAlt - зі - XT.

= ieA2t + зІ7еЗЛі4 - ЗІ -

Ш = ik^-b^-ik + b

ш = і^-зі^ + зі-і-

Приближения z(q)(t) более высоких порядков (q > l) можно получить аналогичным образом при помощи формулы (17).

4. Автоматическое построение вычислительных процедур. Пошаговые методы численного интегрирования системы ОДУ позволяют построить последовательность приближенных значений решения в точках tl = to + hl,. ..,tm = tm-l + hm.

Действие введенных ранее операторов Sm,e и Тм на точных решениях задачи Коши задается явными формулами, и поэтому таблица значений формируется следующим образом:

x(t1) = Am x(t1,to,x(to)),

x(tm) = AMx(tm,tm—1,x(tm—1) ) •

Здесь Am - один из операторов SM}E и Тм.

Рассмотренные в предыдущих пунктах алгоритмы являются настолько общими, что их можно использовать для создания пакета прикладных программ, позволяющего автоматизировать составление универсальных процедур на языках высокого уровня приближенного решения задачи Коши для ОДУ. Один из вариантов такой системы был создан на факультете прикладной математики-процессов управления СПбГУ [4].

Для реализации алгоритмов интегрирования, основанных на разложении решения в ряд, необходима организация хранения массивов, содержащих коэффициенты полиномов в правых частях, обеспечение ввода этих коэффициентов, а также вычисление операторов SM}E и Тм. Предлагаемый генератор программ принимает в качестве исходных данных полиномиальную систему уравнений, а на выходе выдает набор процедур на заданном алгоритмическом языке программирования (C или Fortran), вычисляющих приближенные решения.

Пользователь освобождается от необходимости написания однообразных фрагментов программы, имеет возможность изменить используемые типы данных для представления числовых величин или дополнить генерируемые процедуры в соответствии со спецификой решаемой задачи. Например, рассматривавшийся выше пример с интегрированием уравнений вращательного движения спутника может быть составной частью системы построения программного управления для гашения колебаний искусственного спутника [14].

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

Литература

1. Jorba A., Zou M. A Software Package for the Numerical Integration of ODEs by Means of High-Order Taylor Methods // Experimental Mathematics. 2005. Vol. 14, N 1. P. 99-117.

2. Rail L. B. Automatic Differentiation: Techniques and Applications. Lecture Notes in Computer Science 120. Berlin: Springer Verlag, 1981. 171 p.

3. Parker G. E., Sochacki J. S. Implementing the Picard Iteration // Neural Parallel Sci. Comput. 1996. Vol. 4, N 1. P. 97-112.

4. Бабаджанянц Л. К., Латыпов В. Н. Разработка программно-методического обеспечения по решению обыкновенных дифференциальных уравнений методом малого параметра и рядов Тэйлора. Пилотный проект № 22: Разработка и внедрение инновационной образовательной программы «Прикладные математика и физика». СПб.: Изд-во С.-Петерб. ун-та, 2006. 50 с.

5. Пуанкаре А. О кривых, определяемых дифференциальными уравнениями / пер. с фр. Е. Леон-тович, А. Мейер; под ред. А. А. Андронова. М.: Гостехиздат, 1947. 392 с.

6. Бабаджанянц Л. К. Продолжаемость и представимость решений в задачах небесной механики // Труды Ин-та теор. астрономии АН СССР. 1978. Вып. 17. С. 3-45.

7. Бабаджанянц Л. К., Саркисян Д. Р. Методы рядов Тейлора для динамических систем с управлением: сходимость и оценка погрешности // Современная математика и ее приложения. Т. 24: Динамические системы и оптимизация. 2005. С. 14-34.

8. Алферов Г. В., Бабаджанянц Л. К., Ковригин Д. А., Сенатова С. В. Лабораторный практикум по механике управляемого движения с использованием мини-ЭВМ. Л.: Изд-во Ленингр. ун-та, 1989. 84 с.

9. Бабаджанянц Л. К., Мгоян П. Б. Оценка голоморфных решений обыкновенных дифференциальных уравнений // Вестн. Ленингр. ун-та. Сер. 1: Математика, механика, астрономия. 1984. № 7. С. 9-16.

10. Белецкий В. В. Движение спутника относительно центра масс в гравитационном поле. М.: Изд-во Моск. ун-та, 1975. 308 с.

11. Сарычев В. А., Сазонов В. В. Оценка границы колебательных движений спутника // Препринты Ин-та прикл. математики АН СССР. 1974. № 130. 38 с.

12. Бабаджанянц Л. К., Мгоян П. Б. Оценка голоморфных решений обыкновенных дифференциальных уравнений // Изв. АН Арм. ССР. Сер. математич. 1982. Т. 17, № 2. С. 83-91.

13. Александров А. Ю, Платонов А. В., Старков В. Н., Степенко Н. А. Математическое моде-

лирование и исследование устойчивости биологических сообществ: учебн. пособие. СПб.: Соло, 2006. 186 с.

14. Бабаджанянц Л. К., Потоцкая И. Ю. Управление по критерию расхода в механических системах. СПб.: Изд-во С.-Петерб. ун-та, 2003. 137 с.

Статья рекомендована к печати проф. А. П. Жабко.

Статья принята к печати 25 декабря 2008 г.