Разработка экономичных схем интегрирования структурно разделённых систем обыкновенных дифференциальных уравнений тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат наук Коврижных Николай Александрович

• Коврижных Николай Александрович
• кандидат наук
• 2018

Оглавление диссертации кандидат наук Коврижных Николай Александрович

Введение

Глава 1. Структурные методы интегрирования

1.1 Классы структурных особенностей систем ОДУ

1.1.1 Класс А

1.1.2 Класс В

1.1.3 КлассС

1.2 Метод интегрирования

Глава 2. Структурные особенности математических моделей

небесной механики

2.1 Линейные модели движения ИСЗ

2.2 Модели ограниченной задачи трёх тел

2.3 Задача об электродинамическом тросе на орбите Земли

2.4 Другие модели, имеющие структурные особенности

2.4.1 Сведение уравнений высших порядков к классу В

2.5 Выделение структурных особенностей

Глава 3. Условия порядка для структурных методов

3.1 Общий алгоритм нахождения условий порядка с учётом структуры СОДУ

3.2 Численное решение системы условий порядка

3.3 Использование упрощающих предположений

3.3.1 Описание программы построения системы-следствия

Глава 4. Схемы шестого порядка

4.1 Семейство схем класса В

4.1.1 Алгоритм построения явного метода класса В

4.1.2 Интегрирование с переменным шагом. Оценка локальной погрешности

4.1.3 Схема 11КВ6(4){7Е}

Стр.

4.2 Семейство схем класса C

4.2.1 Схема RKC6(4){8F,7F,7F}

4.3 Численное исследование моделей

4.3.1 Описание программы ode46b

4.3.2 Описание использованных моделей

4.3.3 Обсуждение результатов моделирования

Глава 5. Исследование устойчивости

5.1 Понятие устойчивости в случае структурных СОДУ

5.1.1 Области устойчивости методов класса A

5.1.2 Области устойчивости методов класса C

Глава 6. Непрерывные расширения структурных методов

6.1 Непрерывные методы и задачи с запаздыванием

6.2 Непрерывные структурные методы

6.2.1 Непрерывные методы для задач классов A и B

6.2.2 Непрерывные методы для задач с общей группой уравнений

Заключение

Список литературы

Список рисунков

Список таблиц

Приложение А. Текст программы построения системы условий

порядка

Приложение Б. Текст программы ode46b

Приложение В. Система условий порядка для методов класса C

Приложение Г. Семейство методов шестого порядка класса B

Приложение Д. Семейство методов шестого порядка класса C

Введение диссертации (часть автореферата) на тему «Разработка экономичных схем интегрирования структурно разделённых систем обыкновенных дифференциальных уравнений»

Введение

Системы дифференциальных уравнений являются одним из основных способов математического моделирования. Математическое описание сложных динамических процессов, протекающих в окружающем мире, зачастую приводит к системам, не имеющим аналитического решения. Кроме того, в некоторых задачах (например, в линейных СОДУ большого размера) вывод точного решения может оказаться сопряжён с трудозатратными вычислениями, и получение приближённого решения с некоторой желаемой точностью будет более оправданным. Поэтому численные методы решения систем дифференциальных уравнений всегда будут необходимым инструментом математического моделирования.

В начале XX века немецкий математик Карл Рунге разработал теорию графических численных методов, показав, как с помощью простейших инженерных инструментов на бумаге производить построения, аналогичные сложным математическим операциям вплоть до приближённого интегрирования скалярных дифференциальных уравнений [1]. Развитие вычислительной техники привело к тому, что теория Рунге оказалась не нужна, однако предложенные им идеи привели к созданию явных одношаговых методов Рунге Кутты, позволявших получать приближённое решение вплоть до четвёртого порядка точности при сравнительно небольших трудозатратах: для точности порядка р требовалось р вычислений правой части СОДУ (р-этапные методы одношаговые методы, Р < 4).

Технологический рывок середины XX века принёс новые возможности ЭВМ и вместе с тем предъявил новые требования к точности численного интегрирования. Безрезультатные годы поиска пятиэтапных методов пятого порядка завершились разработкой теории Джона Бутчера, систематизировавшего процесс нахождения методов Рунге Кутты. Бутчер показал, что для высоких порядков существуют ограничения («барьеры Бутчера»): при р ^ 5 методы р-ого порядка требуют не менее р + 1 вычисление правой части СОДУ [2], а при р ^ не мен ее р + 2 [3; 4]. Эти ограничения а также трудоёмкость получения новых методов высокого порядка стали основными аргументами для поиска альтернатив классическим методам типа Рунге Кутты.

Начиная с 1970-х лет наступает время бурного развития новых методов интегрирования и способов сравнения одних методов с другими. Важными свойствами становятся не только порядок точности, но также возможность автоматического управления длиной шага и устойчивость методов. Эрвин Фель-берг одним из первых предложил идею вложенных методов Рунге Кутты, позволяющих существенно экономить вычисления с помощью автоматического управления длиной шага интегрирования [5]. Благодаря этому методы Рунге Кутты обрели новую популярность, и вскоре на их основе были созданы схемы Дж. Дорманда и П. Принса [6]. Одна из этих схем сейчас является основным интегратором в среде Mat lab [7] (как функция ode45) и в РуШоп-модуле SciPy [8] (как метод RK45). Один из методов Фельберга является основным интегратором в среде Maple [9] (как метод rkf45).

Один из способов построения новых эффективных методов был предложен 14. В. Олемским [10]. Он заключается в преодолении «барьеров Бутчера» с помощью анализа структуры СОДУ, разбиения переменных на группы и модификации классической схемы Рунге Кутты. Оказалось, что для многих задач, особенно в области механики, возможно сократить количество вычислений правой части системы, сохранив порядок точности. На основе структурного подхода было построено несколько вложенных схем, превосходящих по своим характеристикам известные аналоги. Следует отметить, что поиск таких методов связан с решением громоздких систем алгебраических уравнений, поэтому представляет особый интерес модификация теории Бутчера на случай структурных методов.

Также представляет интерес возможность применения идей структурного подхода не только для численного решения СОДУ. К примеру, в последние десятилетия получили распространение методы решения задач с запаздыванием на основе т.н. непрерывных методов Рунге Кутты (в этой области стоит отметить работы М. Ценнаро [11; 12]). Как и в уже описанном случае обыкновенных дифференциальных уравнений, для систем специального вида существует возможность модифицировать классическую схему интегрирования, чтобы сократить объём вычислений.

Целью работы является проведение исследований, направленных на развитие теоретических подходов и эффективных вычислительных методов, служащих базовым инструментом математического и компьютерного моделирования динамических объектов в задачах небесной механики.

Научная новизна: проведён анализ структурных особенностей систем дифференциальных уравнений, описывающих математические модели динамики различных задач небесной механики; построены методы шестого порядка для различного вида структурных особенностей СОДУ, требующие меньших вычислительных затрат по сравнению с известными методами того же класса; теория структурных методов применена к системам уравнений с запаздыванием.

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

Методология и методы исследования. Методологической основой работы служат общая теория явных методов Рунге Кутты, теория помеченных деревьев Дж. Бутчера, теория устойчивости численных методов интегрирования.

Основные положения, выносимые на защиту:

1. Проанализированы структурные особенности класса математических моделей небесной механики и продемонстрирован способ модификации этих моделей, позволяющий строить экономичные численные методы их исследования.

2. Разработан общий алгоритм формирования системы условий порядка для экономичных методов численного интегрирования. Построен метод приведения полной системы условий порядка к системе-следствию с более простой структурой связей. На их основе реализован комплекс программ, позволяющий автоматически строить и упрощать системы условий порядка для произвольных структурных классов.

3. Найдены семейства решений систем условий шестого порядка. На их основе построены вложенные схемы шестого порядка точности, требующие меньших вычислительных затрат по сравнению с уже существующими методами. Это преимущество продемонстрировано с помощью разработанной программы, реализующей предложенную схему вместе с алгоритмом автоматического выбора шага интегрирования.

4. Предложена методика исследования устойчивости структурных методов численного интегрирования. Построены области устойчивости некоторых методов.

5. Построены непрерывные экономичные методы интегрирования для моделей, описываемых системами с запаздыванием, и задач, требующих непрерывного приближения к точному решению.

Достоверность полученных результатов обеспечивается сравнением результатов работы построенных методов с известными существующими методами других авторов а также сравнением вычисленных параметров модели с полученными аналитически.

Апробация работы. Основные результаты работы докладывались на международной научной конференции по механике «VIII Поляховские чтения» (г. Санкт-Петербург, февраль 2018 г.), международной конференции «International Conference on Computational Science and Applications» (г. Триест, Италия, июль 2017 г.) и на XLVI, XLVII и XLVIII международных научных конференциях «Процессы управления и устойчивость» (г. Санкт-Петербург, апрель 2015, 2016 и 2017 гг.)

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

Публикации. Основные результаты но теме диссертации изложены в 8 печатных изданиях, 2 из которых изданы в журналах, рекомендованных ВАК, 5 в тезисах докладов.

Объем и структура работы. Диссертация состоит из введения, пяти глав, заключения и трёх приложений. Полный объём диссертации составляет 165 страниц, включая 20 рисунков и 15 таблиц. Список литературы содержит 45 наименований.

Глава 1. Структурные методы интегрирования

1.1 Классы структурных особенностей систем ОДУ

Исторически первым использованием структуры системы обыкновенных дифференциальных уравнений для построения экономичных схем численного интегрирования можно считать методы Нюстрёма. В случае системы уравнений второго порядка у" = /(х,у,у') алгоритм, предложенный Нюстрёмом, позволяет сократить машинную память, используемую при вычислениях. Однако настоящее преимущество возникает в том случае, когда правая часть не зависит от первой производной у'. Для систем вида

удаётся сократить также и количество вычислений функции /, которые, как правило, занимают основную часть времени вычислительного процесса. Например, в то время, когда классические методы Рунге Кутты для достижения пятого порядка точности требуют минимум шести вычислений правой части, метод Нюстрёма требует только четырёх, что даёт экономию в полтора раза.

Под термином «система ОДУ со структурными особенностями» обычно полагают то, что рассматриваемая система принадлежит к одному из трёх классов А В или С.

К классу А относятся системы с перекрёстной структурой зависимостей:

у" = / {х,у)

(1.1)

1.1.1 Класс А

(1.2)

ж е [Хо,Хк] с К, уа :[Хо,Хк] , 5 = 1,2,

: [Хо,Хк] х ]Г2 ,

/2 : [Хо,Хк] х № ]Г2.

Можно видеть, что этот класс является обобщением систем второго порядка, не зависящих от первой производной. Каждая система вида у" = /(х,у) с помощью введения обозначений у1 = у7 у2 = у' может быть записана в виде системы первого порядка с перекрёстной системой связей

{

у1 = У2' (1.3)

У2 = ¡2(Х,У1).

Важным качеством при этом является то, что, в отличие от систем второго порядка, две группы переменных могут иметь разную размерность.

1.1.2 Класс В

К классу В относятся системы со следующей структурой зависимостей:

{

У[ = М^Яи ... М-Ш+ъ ... ,Уn), г = 1,2, ...^ У'з = I) (x,Уl,... ,Уз-1), 3 = 1 + 1,...,n,

(1.4)

р = ^ ту, х е [Хо,Хк] СК, уа :[Хо,Хк] Мг-, 5 = 1,...,п,

^=1

-»К

/г : [Хо,Хк] х

£ : [Хо,Хк] х ,

г = у г

У=1

п

Г3 = > Гу, ^=3

г = 1,2,...,/, 3 = I + 1,...,п.

Этот класс является обобщением класса А: в системах такого вида неизвестные двух групп могут зависеть не только от противоположной группы, но и от неизвестных своей группы, имеющих меньший индекс.

1.1.3 Класс С

Следующее обобщение — класс С, к нему относятся системы с т.н. «общей» группой уо:

Уо

< У; —

Уг

у'з

Мх,Уо,

/з {х,Уо,

■ ,Уп),

,Уз-1).

. . , Уп) :

г = 1 , 2 ,..., I , 3 = I + 1 ,..., п ,

(1.5)

(1.6) (1.7)

р = ^гу, х е [Хо,Хк] СМ, ^=0

/с : [Хо,Хк] х Мр ,

ь : [Хо,Хк] х Мр-'г ,

у8 :[Хо,Хк] , 5 = 0,... ,п,

ъ : [Хо,Хк] х Мр

г —

грЗ -

£

У=1

п

е

^=3

г\

гл

г - 1,2,...,/, 3 = I + 1,...,п.

Последовательная вложенность классов А, В и С позволяет использовать одни и те же названия для групп переменных ук\ здесь и далее мы будем называть их нулевой (1.5), первой (1.6) и второй (1.7). В литературе также часто называют нулевую группу общей, а первую и вторую обозначают какг-ю и 3-ю.

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

В случаях классов В и С любая группа может быть вырожденной: например, система вида

У'о = fо{t,yо,Уl), У1 = Шмо)

принадлежит классу С (группа 2 отсутствует).

(1.8)

3

и

Е

а)

б)

в

а) методы Нюстрёма (Е обозначает единичную матрицу); б) класс А; в) класс В; г) класс С Рисунок 1.1 Структуры допустимых зависимостей в СОДУ

1.2 Метод интегрирования

Метод интегрирования систем класса С. Будем считать, что нам известно точное решение у8(х), $ = 0,... ,п, системы (1.5)—(1.7) в точке х € [Х0,Хк]. Не умаляя общности рассуждений, для простоты вывода будем считать г8 = 1, в = 0,... ,п.

Для численного интегрирования систем (1.5) (1.7) рассматривается явный одношаговый метод [13]. В предположении достаточной гладкости правой части рассматриваемой системы приближение у8 к точному решению у8(х + К), в = 1,... ,п в точке х + К € [Х0,Хк] ищется в виде:

то

Уо(х + h) ж у(х + h) = уо(х) + h^bovfavih), (1.9)

v=1

т\

уг(х + h) ж уг(х + h) = Уг(х) + h^2bivkiv{h), i = l,...,l, (1.10)

V=1

m\

yj (x + h) ж yj (x + h) = yj (x) + h^^b2vkj,v(h), j = l + !,...,n, (1.11)

V=1

причём ks,w = ks,w(h) вычисляются в строгой последовательности

ko,i,..., kn,i, ко,2,..., kn,2, ko;i, ki¡3,...

(1.12)

по формулам ко,У — /о (х + C-Qyh,

V— 1

Уо + h^2a00v\ih^(h),

Ц=1

V—1 V—1

У1 + h a01v^k1^(h),..., уг—1 + h ,

Ц=1 Ц=1

V —1 V —1

У1+1 + h^a-02V^+1^(h), ...,Уп + h^a02*^n^(h)) , (1.13)

in

Ц=1 Ц=1

ki,V — X + Clvh,

У0 +

Ц=1

VV

У1 + a11V^hAh),..., уг—1 + h a,11V^h—1^(h),

Ц=1 Ц=1

V —1 V —1

У1+1 + ...,Уп + h^a,12V^n,^(h)) , г — l,...,l,

Ц=1 Ц=1

(1.14)

kj,V — fj ^x + C2Vh,

V

У0 + h a20V^k0,^ (h), ц=1

VV

(h),... , Уг—1 + h^2 a21V\ikl,^ (h),

ц=1 ц=1

V V

yi+1 + h ^ ^22V^^/+1,^(h), ...,Уп + h ^ a22V^kj—1^(h)j , j — I + 1,... ,n.

Ц=1 Ц=1

(1.15)

Параметрами метода являются коэффициенты aWlW2Vц (матрицы этих коэффициентов будем обозначать как AWlW2), bWlV, cWlV, w1,w2 G {0,1,2}. Говорят, что метод имеет порядок р, если

\у8(х + h) — у„1 « 0(h^+1), s — 0,...,п.

Строгий порядок вычисления значений к связывает количество этапов по каждой группе: ту ^ т0 ^ ту + 1. Основная идея структурного подхода заключается в алгоритмическом использовании структурных особенностей математической модели для сокращения числа обращений к правым частям первой и второй групп уравнений, поэтому мы будем считать далее т = т0 = ту + 1.

Для того, чтобы метод был явным, необходимо, чтобы матрицы Ау0, А207 Ац, А'ц и А22 были нижнетреуголъными, а матрицы А00, А0у, А02 и Ау2 строго нижнетреу зольными. В этом случае к моменту вычисления каждого из выражений (1.13) (1.15) аргументы в их правых частях будут уже известны. В табл. 1 приведён компактный общий вид коэффициентов явного метода при т = т0 = ту + 1 (для удобства индексы отделены запятыми). Дополнительное требование ацц = аюи = 0 обусловлено часто используемым предположением

сШ1у = аШ1 'Ш1,'Ш2 € {0,1,2} и тем, что аут = 0.

ц

Замечание 1. В случае, когда система состоит только из нулевой группы, метод вырождается в классический метод Рунге Кутты с коэффициентами А = А00, Ь = Ь0, с = с0 и количеством этапов в = т0:

Замечание 2. В силу того, что системы класса В отличаются от систем класса С лишь отсутствием нулевой группы, рассмотренный метод можно при-

0

же касается и систем класса А достаточно вдобавок исключить все параметры, использующие зависимости внутри первой группы и внутри второй. Табличные записи коэффициентов для обеих схем будут содержать коэффициенты ЬУУ7 Ь2у, С-1У, с2у, а\2уЦ) а2\уЦ) а для схемы класса В дополиительно а22уц- Требо-

вания к треуголыюсти матриц для явных методов останутся теми же.

в

у(х + К) ~ у(х + К) = у(х) + К^^ Ъуку(Ь),

(1.16)

у=1

Таблица 1 — Коэффициенты явного метода класса £ при т = т0 = т\ + 1

С0,у «0,0,у,ц. «0,1,у,ц. «0,2,у,ц.

0 0 0 ... 0 0 0 ... 0 0 0 ... 0 &0,1

С0,2 «0,0,2,1 0 ... 0 «0,1,2,1 0 ... 0 «0,2,2,1 0 ... 0 &0,2

С0,т-1 «0,0,т-1,1 «0,0,т—1,2 ... 0 «0,1,т—1,1 «0,1,т—1,2 . . . 0 «0,2,т—1,1 «0,2,т—1,2 . . . 0 &0,т-1

С0,т «0,0,т,1 «0,0,т,2 . . . «0,0,т,,т—1 «0,1,т,1 «0,1,т,2 . . . «0,1,т,,т—1 «0,2,т,1 «0,2,т,2 . . . «0,2,т,т-1

С1,у

0 0 0 ... 0 0 0 ... 0 0 0 ... 0 ^1,1

<Л,2 «1,0,2,1 «1,0,2,2 . . . 0 «1,1,2,1 «1,1,2,2 . . . 0 «1,2,2,1 0 ... 0 &1,2

С1,т-2 «1,0,т-2,1 «1,0,т—2,2 . . . 0 «1,1,т—2,1 «1,1,т—2,2 . . . 0 «1,2,т—2,1 «1,2,т—2,2 . . . 0 &1,т-2

С1,т-1 «1,0,т-1,1 «1,0,т—1,2 ... «1,0,т- 1,т— 1 «1,1,т—1,1 «1,1,т—1,2 . . . «1,1,т— 1,т— 1 «1,2,т—1,1 «1,2,т—1,2 . . . 0 &1,т-1

С2,у «2,0,у,ц.

С2,1 «2,0,1,1 0 ... 0 «2,1,1,1 0 ... 0 «2,2,1,1 0 ... 0 &2,1

С2,2 «2,0,2,1 «2,0,2,2 . . . 0 «2,1,2,1 «2,1,2,2 . . . 0 «2,2,2,1 «2,2,2,2 . . . 0 &2,2

С2,т-2 «2,0,т-2,1 «2,0,т—2,2 . . . 0 «2,1,т—2,1 «2,1,т—2,2 . . . 0 «2,2,т—2,1 «2,2,т—2,2 . . . 0 &2,т-2

С2,т-1 «2,0,т-1,1 «2,0,т—1,2 ... «2,0,т-1,т-1 «2,1,т—1,1 «2,1,т—1,2 . . . «2,1,т-1,т-1 «2,2,т—1,1 «2,2,т—1,2 . . . «2,2,т-1,т-1 &2,т-1

Глава 2. Структурные особенности математических моделей

небесной механики

Поведение многих математических моделей (например, движение материальной точки в гравитационном или электростатическом поле) описывается дифференциальными уравнениями второго порядка вида (1.1). Как было показано в разд. 1.1.1, каждая такая система может быть приведена к перекрёстному виду (1.2) (класс А) с помощью введения новых переменных. Далее речь пойдёт о моделях, описываемых системами более сложного вида.

Основное ограничение, накладываемое на СОДУ классов А и В — отсутствие зависимостей производных от переменных с тем же индексом В уравнениях движения многих моделей небесной механики производные второго порядка зависят от производных первого порядка. Однако, как правило, эти связи обусловлены силами, появляющимися из-за выбора иеинерциальной системы отсчёта, и компоненты вектора ускорений не зависят от соответствующих компонент вектора скоростей. К примеру, сила Кориолиса^ = так = -2т[шху] направлена перпендикулярно вектору скорости V, поэтому в системе координат, вращающейся вместе с радиус-вектором тела, компоненты вектора ак не зависят от соответствующих компонент вектора V. Это даёт возможность переопределить переменные так, чтобы структура зависимостей приобрела нужный вид. Следующие ниже примеры иллюстрируют этот приём.

2.1 Линейные модели движения ИСЗ

Модель Клохесси—Уилтшира [14] (С\¥). Это классическая модель движения объекта в окрестности искусственного спутника Земли (ИСЗ). Спутник А совершает движение по круговой орбите £ вокруг земли Е. Выбирается правая система прямоугольных координат так, чтобы ось Оу была сонаправле-на мгновенной скорости, а ось Ох была направлена вдоль радиусаЕА (рис. 2.1).

Рисунок 2.1 Движение в окрестности ИСЗ

Тогда линеаризованная система уравнений движения спутника В с радиус-вектором р = {х,у,г} будет иметь вид:

о

х = 2шу + 3ш2х, у = — 2ш±,

г = — Ш

(2.1)

Перестановка переменных сводит систему (2.1) к системе класса А перекрёстной). Матрицы зависимостей до и после перестановки приведены на рис. 2.2. Первая группа: щ = х,и2 = у,Щ = вторая группа: VI = х,ю2 = у,у3 = ¿. Производные переменных первой группы зависят только от переменных второй группы, и наоборот.

Модель Чаунера—Хемпеля [15] (ТН). Модификация модели С\¥, в которой основная орбита рассматривается как эллиптическая с эксцентриситетом е. Истинная аномалия 6 служит независимой переменной:

3х

х =

1 + е сое 6 у = — 2±,

+ 2у,

(2.2)

г = —

Модель Швайгарта—Седвика [16] (ББ). Модификация модели С\¥, в которой учитывается влияние сжатия Земли на движение спутников.

х

= 2шсу + (5с2 — 2)ш2х,

< у = —2ысх, г = — д2г + 21д + ф).

(2.3)

X у X X у X

X

У

X

X

У

X

>

X у X X у X

X

У

X

X

У

X

Рисунок 2.2 Перестановка переменных для моделей С\¥, ТН и ББ

Системы ТН и ББ имеют ту же структуру зависимостей, что и модель С\¥, поэтому аналогичной перестановкой переменных сводятся к системам классаА.

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

2.2 Модели ограниченной задачи трёх тел

Ограниченная задача трёх тел часто рассматривается в аналогичной системе координат, ось Ох которой проходит через два основных тела М\ и М2, ось Ох направлена вдоль вектора п угловой скорости их обращения друг относительно друга, а ось Оу выбирается так, чтобы образовывать правую тройку. Пусть тела имеют массы ту т2 и т. Введём обозначения для координат тел: М\(х\, 0, 0) М2(х2, 0, 0) и М (х,у,х) и для отрез ков: М\М = г\ =

\/(х — х\)2 + у2 + М2М = г2 = ^(х — х2)2 + у2 + х2. В этом случае уравнения движения будут выглядеть [17] следующим образом:

( п . 2 ди

х = 2пу + п х + -7—,

ох

2 ди , ч

\у = —2пх + п2у + -—-, (2.4)

оу

.. _ ди

^ дх'

где и = /т(у1 + ^) ^силовая функция, зависящая только от координат.

Перестановка переменных сводит эту систему к системе класса В. Матрицы зависимостей до и после перестановки приведены на рис. 2.4. Первая группа: щ = х,и2 = = у, вторая группа: уг = у,У2 = х,юз = ¿.

Модель движения вокруг точки либрации Ь1 [18]. Рассматривается движение космического аппарата (КА) в окрестности точки либрации Ьг системы Солнце Земля. В качестве системы координат выбирается геоцентрическая вращающаяся система, аналогичная изображённой на рис. 2.1, в которой осъОх направлена на Солнце, масштаб расстояния выбран так, чтобы точка Ьг имела координаты (1, 0, 0), а масштаб времени^так, чтоб период обращения Земли вокруг Солнца составлял 2п. В этом случае линейное приближение уравнений движения К А описывается системой:

Хг = Ж2 + У1, Х2 = -Хг + У2, Хз = Уз,

3хг _ , ,

\Уг = —л^з + 2хг + У2, (2.5)

II ^ II

3^2

У2 = - ^2 - Уг,

ЦхЦ3

3хз

Уз = -^"ич - хз, ЦхЦ3

где (хг,х2,хз) — вектор координат КА, (уг, у2, уз) — вектор импульсов. Перестановка переменных сводит её к системе класса В. Матрицы зависимостей до и после перестановки приведены на рис. 2.5. Первая группа: иг = хг,и2 = хз,из = у2) вторая группа: уг = х2,у2 = уг,уз = уз.

Модель Аренсторфа [19]. Рассматривается плоское движение космического аппарата с координатами (хг,х2) в гравитационном иоле, создаваемом Землей (0, 0) и Луной (1, 0).

где

, 0 . ,хг + ц хг - ц

хг = хг + 2^2 - ц —=;--ц—-—

иг ^2

0. ,Х2 Х2 Х2 = ^2 - 2хг - ц — - ц—,

иг ^2

Вг = ((жг + ц)2 + ^2)з/2 , ^2 = ((жг - Ц)2 + х^ ц = 0,012277471, ц' = 1 - ц.

(2.6)

-1 -0.5 0 0.5 1

Рисунок 2.3 Орбита Аренсторфа

Перестановка переменных сводит эту систему к системе класса В. Матрицы зависимостей до и после перестановки приведены на рис. 2.6. Первая группа: щ = Х\,и2 = ¿2, вторая группа: = х2,у2 = У\- При должных начальных условиях КА движется по периодической орбите (см. рис. 2.3), которая получила имя Ричарда Аренсторфа и стала основной для программы «Аполлон».

2.3 Задача об электродинамическом тросе на орбите Земли

Тело, движущееся в магнитном поле, испытывает действие силы Лоренца

р = +[у х В]),

X у X X у X

X

У

X

X

У

х

>

X X у У X X

X

X

У

У

X

X

Рисунок 2.4 Перестановка переменных для системы (2.4)

Х1 Х2 Хз У1 У2 Уз

Х\

Х'1

Хз

У1

У2

Уз

>

Х1 Хз У2 У1 Уз

Х\

Хз

У2

Х'2

У\

Уз

Рисунок 2.5 — Перестановка переменных для модели Ь\

Х\ Х2 Х\ х2

Х\

Х'2

XI

Х2

>

Х\ Х2 Х2 Х\

Х\ Х2

Х2 Х\

Рисунок 2.6 Перестановка переменных для модели Аренсторфа

е ф ре Рф

е

ф

Ре

Рф

>

е рф ф Ре

е рф

ф Ре

Рисунок 2.7 Перестановка переменных для моделей электродиманического

троса

и, как и в случае с силой Кориолиса, магнитная составляющая перпендикулярна вектору скорости. Это объясняет, почему системы уравнений, описывающие модели, учитывающие силу Лоренца или силу Ампера, удаётся свести к классу В. К примеру, задача моделирования системы из двух ИСЗ, соединённых электродинамическим тросом в геомагнитном поле, описывается системой [20]:

е =

Ре

СОЙ2ф

- 1,

Ф = рф.

3

ре = —- cos2 ф sin 2е + те(е,ф), 2

(2.7)

Рф

2

рее

3

tg ф — - cos2 е sin 2ф + тФ(е,ф).

cos2 ф 2

Перестановка переменных сводит её к системе класса B. Матрицы зависимостей до и после перестановки приведены на рис. 2.7. Первая группа: щ = е,и2 = рф, вторая группа: V\ = ф,^2 = ре.

В работе [21] подобная модель описана системой

Ре + 1

е=

СОв2ф

ф

Ре

= рф-,

3

—- cos2 ф sin2е + де(е,ф), 2

(2.8)

Рф =

2

ре

3

2 tg ф — - cos е sin 2ф + ^ф(е,ф).

cos2 ф 2

к которой можно применить ту же перестановку переменных.

2.4 Другие модели, имеющие структурные особенности 2.4.1 Сведение уравнений высших порядков к классу В

Однородное дифференциальное уравнение

у(п) + ап-1(х)у(п-1) + ... + а1(х)уГ + ао(х)у = 0 (2.9)

с помощью замены у = и ехр ^о МОжет быть избавлено от производной

(п — 1)-го порядка. Если в получившемся уравнении

и(п) + Ьп—2(х)и(п—2) + ... + Ь1(х)и' + Ъо(х)и = 0 (2.10)

производные чётных и нечётных порядков отнести к двум группам, то соответствующая система первого порядка будет принадлежать к классу В. В самом деле, если обозначить ук = и(к\ то при к < п ук зависит от (т.е. от неизвестной из другой группы), а уп будет зависеть от производных другой группы и производных своей группы, имеющих меньший порядок.

Результат замены представлен на рис. 2.8. В зависимости от чётности п параметр уп—1 попадёт в первую или вторую группу. Соответствующая строчка матрицы зависимостей будет заполнена целиком, кроме диагонального элемента, поскольку уп—1 не зависит от уп—1 (в выр. 2.10 не зависит от и(п—1^). Это значит, что новая система принадлежит классу В.

Vо

щ

■Ьп—2 Уп-1

Рисунок 2.8 Замена переменных для однородного уравнения высокого

порядка

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

У{) У1 У2 ... Уп-1

У0 У2 ... У1 Уз ...

Уо

щ

Уз

>

силы, зависящие от ускорений, редки. Примером может служить сила Абрага-ма—Лоренца, действующая на точечный заряд со стороны поля, вызванного его собственным неравномерным движением:

р 2 (¡2 (!а

" 3с3

где а — вектор ускорения. При учёте этой силы закон движения становится системой уравнений третьего порядка [22].

Модель Курамото—Сивашинского [23]. Г. И. Сивашинский разработал модель распространения пламени, которая в далынейшем была сведена И. Курамото к уравнению, описывающему фазовую турбулентность в реакционно диффузионных системах:

щ + У2и + У2и + 1 |Уи|2 = 0. (2.11)

2

С помощью замен и(х^) = —(ж) и у = Ь(х) уравнение (2.11) сводится

к обыкновенному дифференциальному уравнению:

2

у'" + у' = с2 — \, (2.12)

где

ний (2.12) посвящено несколько работ. Поскольку уравнение имеет вид (2.10) (вторая производная отсутствует), замена и\ = у,и2 = у'',У\ = у',и2 = у'" приводит его к системе класса В.

Трёхслойная балка. Рассматривается трёхслойная балка составлена из трёх различных материалов. Было показано [24], что отклонение^ подчиняется уравнению третьего порядка

—^ — к2 + а = 0, 2.13

ах3 ах

где к2 и а ^физические параметры системы, зависящие от упругости слоёв.

Как и в предыдущем случае, отсутствие второй производной позволяет свести (2.13) к системе класса В.

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

с1 Ж п(1> X , с1 X , О- л \

тИ* + ? Ж + кл +кГ = °, (2'14)

где т — масса клапана, / — сила тренпя, к — коэффициент упругости пружины клапана, I — момент инерции турбины. [25]

Это однородное уравнение третьего порядка, включающее производную второго порядка. Следовательно, оно может быть сведено к виду (2.10) и далее

В

2.5 Выделение структурных особенностей

Список литературы диссертационного исследования кандидат наук Коврижных Николай Александрович, 2018 год

Список литературы

1. Runge, С. Graphical Methods: A Course of Lectures Delivered in Columbia University, New York, October, 1909, to January, 1910 / C. Runge. Columbia University Press, 1912. - (... Graphical Methods).

2. Butcher, J. C. On Runge Kutta Processes of High Order / J. C. Butcher. -1964. - May.

3. Butcher, J. C. On the Attainable Order of Runge Kutta Methods / J. C. Butcher // Math. Comput. - 1965. Vol. 19. - P. 408 417.

4. Butcher, J. C. The non-existence of ten Stage eight Order Explicit Runge Kutta Methods / J. C. Butcher // BIT. - 1985. - Vol. 25. -P. 521 - 540.

5. Fehlberg, E. Classical fifth-, sixth-, seventh-, and eighth-order Runge Kutta formulas with stepsize control / E. Fehlberg // NASA Technical Report 287. -1968.

6. Dormand, J. R. A Family of Embedded Runge Kutta Formulae / J. R. Dor-mand, P. J. Prince // Journal of Computational and Applied Mathematics. -1980. - Mar. - Vol. 6. - P. 19 - 26.

7. Choose an ODE Solver [Электронный ресурс]. URL: https : / / www . mathworks. com/help/matlab/math/choose- an-ode- solver. html (дата обр. 11.08.2018).

8. scipy.integrate.solve_ivp [Электронный ресурс]. URL: https://docs.scipy. org/doc/sinpy/reference/generated/scipy.integrate.solve_ivp.html (дата обр. 13.08.2018).

9. dsolve/numeric/rkf45 [Электронный ресурс]. URL: https : / / www . maplesoft. com / support / help / maple / view. aspx ? path — dsolve % 2Frkf45 (дата обр. 06.09.2018).

10. Олемскощ И. В. Численный метод интегрирования систем обыкновенных дифференциальных уравнений / И. В. Олемской // Математические методы анализа управляемых процессов. 1986. С. 157 160.

11. 0wren, В. Derivation of efficient continuous explicit Runge Kutta methods /

B. О wren, M. Zennaro // SIAM J. Sci. Stat. Comput. - 1992. - Vol. 13, no. 6. P. 1488 - 1501.

12. В ell en, A. Numerical methods for delay differential equations / A. Bellen, M. Zennaro. 1st ed. - Oxford : Oxford Science Publications, Clarendon Press, 2003.

13. Олемской, И. В. Структурный подход в задаче конструирования явных одношаговых методов / 14. В. Олемской // Журнал вычислительной математики и математической физики. 2003. Т. 43, вып. 7. С. 961 974.

14. Clohessy, W. Н. Terminal Guidance System for Satellite Rendezvous / W. H. Clohessy, R. S. Wiltshire // J. Astronaut. Sci. - 1960. - Vol. 27, no. 9. P. 653 - 678.

15. Tschauner, J. Rendezvous Zu Einem In Elliptischer Bahn Umlaufenden Ziel / J. Tschauner, P. Hempel // Astronáutica Acta. 1965. T. 11, № 2.

C. 104 109.

16. Schweighart, S. A. High-Fidelity Linearized J2 Model for Satellite Formation Flight / S. A. Schweighart, R. J. Sedwick // Journal of Guidance, Control, and Dynamics. 2002. - Vol. 25, no. 6. - P. 1073 - 1080.

17. Лукьянов, Л. Лекции но небесной механике / Л. Лукьянов, Г. Ширмин. Алматы: Эверо, 2009. 227 с.

18. Шиманчук, Д. В. Построение траектории возвращения в окрестность кол-линеарной точки либрации системы Солнце Земля / Д. В. Шиманчук, А. С. Шмыров // Вести. С.-Петерб. ун-та, Сер. 10. 2013. Вып. 2.

С. 76 85.

19. Хайрер, Э. Решение обыкновенных дифференциальных уравнений. Нежёсткие задачи / Пер. с англ. 14. А. Кульчицкой, С. С. Филиппова под ред. С. С. Филиппова. / Э. Хайрер, С. Нёрсетт, Г. Ваннер. М.: Мир, 1990. 512 с.

20. La/rsen, М. В. Control by damping injection of electrodynamic tether system in an inclined orbit / M. B. Larsen, M. Blanke // 2009 American Control Conference. 2009. - P. 4824 - 4829.

21. Lars en, M. B. Passivity-Based Control of Rigid Electrodynamic Tether / M. B. Larsen, M. Blanke // Journal of Guidance, Control, and Dynamics. -2011. - Vol. 34, no. 1. - P. 118-127.

22. Ландау, Л. Теория ноля / Л. Ландау, Е. Лифшиц. Наука; Глав. ред. физико-математической лит-ры, 1967. (Теоретическая физика : в десяти томах / Л.Д. Ландау и Е.М. Лифшиц).

23. С. Troy, W. The existence of steady solutions of the Kuramoto-Sivashinsky equation / W. C.Troy //. Vol. 82. - 1989. P. 269 - 313.

24. Krajcinovic, D. Sandwich Beam Analysis / D. Krajcinovic //. Vol. 39. 1972. - P. 773 - 778.

25. Gregus, M. Third Order Linear Differential Equations / M. Gregus. D. Reidel Publishing Company, 1987. (Mathematics and its Applications).

26. Олемской, И. В. Модификация алгоритма выделения структурных особенностей / И. В. Олемской // Вести. С.-Петерб. ун-та, Сер. 10. 2006. Вып. 2. С. 46 54.

27. Еремин, А. С. Разработка явного одношагового вложенного метода для систем структурно разделенных обыкновенных дифференциальных уравнений : дис. ... кандидата физико-математических наук: 01.01.07 / А. С. Еремин. С.-Петерб. гос. ун-т, Санкт-Петербург, 2009. 91 с.

28. Коврижмых, Н. А. Об оптимальном выборе параметра регуляризации в алгоритме Левенберга Марквардта / Н. А. Коврижных // Процессы управления и устойчивость. 2015. Т. 2. С. 426 431.

29. Коврижных, Н. А. Вложенный шестиэтапный метод шестого порядка точности интегрирования систем структурно разделённых дифференциальных уравнений / Н. А. Коврижных // Процессы управления и устойчивость. 2016. Т. 3. С. 183 187.

30. Олемской, И. В. Семейство шестиэтапных методов шестого порядка / И. В. Олемской, Н. А. Коврижных // Вести. С.-Петерб. ун-та. 2018. Т. 14, вып. 8. С. 215 229.

31. Kovrizhnykh, N. A. On a Two Families of Efficient Fifth Order Schemes for Solving ODE Systems / N. A. Kovrizhnykh, A. S. Eremin // AIP Conference Proceedings. - 2018. - Vol. 1959, no. 1. - P. 030014.

32. Олемской, И. В. Методы интегрирования систем структурно разделенных дифференциальных уравнений / И. В. Олемской. СПб.: Изд-во С.-Петерб. ун-та, 2009. 180 с.

33. Сравнительное исследование преимуществ структурных методов численного решения обыкновенных дифференциальных уравнений / В. П. Бубнов [и др.] // Труды СПИИРАН. 2017. Т. 53, вып. 4. С. 51 72.

34. Tsitouras, С. Cheap Error Estimation for Runge Kutta methods / C. Tsi-touras, S. N. Papakostas // Siam Journal on Scientific Computing. - 1999. -Vol. 20, issue 6. - P. 2067- 2088.

35. Verner, J. Strategies for Deriving New Explicit Runge Kutta Pairs / J. Verner // Annals of Numerical Mathematics. - 1994. - Vol. 1.

P. 225 - 244.

36. El-Mikkawy, M. E. A. A General Four-Parameter Non-FSAL Embedded Runge Kutta Algorithm of Orders 6 and 4 in Seven Stages / M. E. A. El-Mikkawy, M. M. M. Eisa // Applied Mathematics and Computation. 2003. - Vol. 143, no. 2. P. 259 - 267.

37. Фалейчик, Б. В. Одношаговые методы численного решения задачи / Б. В. Фалейчик. Минск: БГУ, 2010. 42 с.

38. Т. Е. Hull [et al.] // SIAM J. Numer. Anal. - 1972. - Vol. 9(4). P. 603 - 637.

39. Хайрер, Э. Решение обыкновенных дифференциальных уравнений. Том 2: Жёсткие и дифференциально-алгебраические задачи / Пер. с англ. Е. Л. Старостина, И. А. Кульчицкой, А. В. Тыглияна и С. С. Филиппова под ред. С. С. Филиппова. / Э. Хайрер, Г. Ваннер. М.: Мир, 1999. 685 с.

40. Целищев, С. О. Об устойчивости структурного метода решения систем ОДУ / С. О. Целищев, А. С. Еремин // Процессы управления и устойчивость: Труды 42-й международной научной конференции аспирантов и студентов / под ред. А. С. Ерёмина, Н. В. Смирнова. 2011. С. 207 212.

41. Olemskoy, I. V. Embedded Methods of Order Six for Special Systems of Ordinary Differential Equations / I. V. Olemskoy, A. S. Eremin, N. A. Kovrizh-nykh // Appl. Math. Sci. - 2017. - Vol. 11(1). P. 31 38.

42. Calvo, M. A New Embedded Pair of Runge Kutta Formulas of Orders 5 and 6 / M. Calvo, J. I. Montijano, L. Randez // Computers Math. Applic. -1990. - Vol. 20, no. 1. - P. 15 - 24.

43. Винничек, H. H. Исследование устойчивости структурных методов интегрирования систем обыкновенных дифференциальных уравнений / Н. Н. Винничек, Н. А. Коврижных // Процессы управления и устойчивость. 2017. Т. 4. С. 149 153.

44. Олемскощ И. В. Вложенные методы пятого порядка / 14. В. Олемской // Вести. С.-Петерб. ун-та, Сер. 10. 2004. Вып. 2. С. 82 93.

45. Eremin, A. S. Continuous Extensions for Structural Runge Kutta Methods / A. S. Eremin, N. A. Kovrizhnykh // Computational Science and Its Applications ICCSA 2017. - Cham : Springer International Publishing, 2017.

P. 363 - 378. - (Lecture Notes in Computer Science ; 10405).

Список рисунков

1.1 Структуры допустимых зависимостей в СОДУ..........................11

2.1 Движение в окрестности ИСЗ ............................................16

2.2 Перестановка переменных для моделей С\¥, ТН и ББ..................17

2.3 Орбита Аренсторфа........................................................19

2.4 Перестановка переменных для системы (2.4)............................20

2.5 Перестановка переменных для модели Ь\................................20

2.6 Перестановка переменных для модели Аренсторфа....................20

2.7 Перестановка переменных для моделей электродиманического троса 21

2.8 Замена переменных для однородного уравнения высокого порядка . 22

3.1 Использование упрощающих предположений............................29

4.1 Зависимость глобальной погрешности от количества вычислений

для модели 1................................................................52

4.2 Зависимость глобальной погрешности от количества вычислений

для модели 2................................................................53

4.3 Зависимость глобальной погрешности от количества вычислений

для модели 3................................................................54

1 2 3

5.1 Области устойчивости методов порядка 5................................60

5.2 Области устойчивости методов порядка 6................................61

5.3 Области устойчивости (а^ ,2 = п)......................................62

5.4 Области устойчивости (^п ^ а^ ,2 ^ п)............................62

Список таблиц

1 Коэффициенты явного метода класса С при т = т0 = т\ + 1 . . . . 14

2 Размеры системы условий порядка........................................25

3 Примеры построения условий порядка ..................................26

4 Коэффициенты метода ЯКВ6(4){7Р}....................................43

5 Коэффициенты метода 11КС6(4){8,7,7}Р................................48

1

7 Точность при фиксированном количестве шагов для модели 2 ... . 55

8 Точность при фиксированном количестве шагов для модели 3 . . . . 55

9 Коэффициенты схемы С11КОХ4с13и(2){4Р,3} ..........................67

10 Коэффициенты схемы СЯК0X4(3) {5Р} ................................68

11 Коэффициенты С11КОТЗс12и{2}..........................................68

12 Коэффициенты С11КОСЗс12и{3,2}........................................69

13 Коэффициенты схемы С11КОСЗ(2){4Р,ЗР}..............................69

14 Коэффициенты С11КОС4с13и{4,3}........................................70

15 Сравнение непрерывных методов..........................................70

Приложение А Текст программы построения системы условий порядка

Листинг А.1 Класс помеченных деревьев # -*- coding: utf-8 -*-

from copy import deepcopy greeks = ['\\nu '\\mu '\\xi

10

15

20

25

Wpsi

Wphi ']

class TreeParameters(object) max.order = 6

dependencies = {

[0 , 1 , 2] , [0 , 1 , 2] , [0, 1, 2]>

: {0: 1 , 1 : 1 , 2 : 1>,

: {0: 0 , 1 : 0 , 2 : 1>,

: {0: 0 , 1 : 0 , 2 : 0»

strictly_triangular 0 1 2

char_a = {}

for i in range(0, len(dependencies)): char_a[i] = {}

for j in range(0, len(dependencies [i])) : char_a[i][dependencies[i][j]] = \

' A{i}{j} '.format(i=i, j=dependencies[i][j])

char_b = [ 'B{i} format(i = i) for i in range(0, len(

dependencies))] char_c = [ 'C{i} format(i = i) for i in range(0, len( dependencies))]

t

up_root = {

0: {0: 1, 1: 2, 2: 2}, 1 : {0 : 1 , 1 : 2 , 2 : 1},

40

45

50

55

60

65

70

2 {0 1 , 1 2 , 2 2»

root2 = {

0 {0 3 , 1 3 , 2 3},

1 {0 3 , 1 3 , 2 1>,

2 {0 3 , 1 3 , 2 2»

leaf = {

0 {0 2 , 1 1 , 2 1>,

1 {0 2 , 1 1 , 2 1>,

2 {0 2 , 1 1 , 2 1»

def tex_term(t, h, i, j=Mone): cl = str (i) al = greeks[h-1]

s = '\sum\limits_{ ' + al + '} '

if t [0] == 'b ' :

return s + 'b_{' + cl + al + '}'

if t [0] == ' a ' : c2 = str (j) aO = greeks[h-2]

return s + 'a_{' + cl + c2 + aO + al + '}

if t [0] == ' c ' :

return 'c_{' + cl + al + '}'

class Tree(object) :

def __init__(self , source, parent = Mone) : self.color = source [0] self . children = [] self.parent = parent

for child in source [1:]:

self.children.append(Tree(child))

def __str__(self) :

if len(self.children) > 0:

s = ')( ' .join(str(ch) for ch in self.children)

85

90

95

100

105

return str(self . color) + '(' + s + ')' else :

return str ( self . color)

def clone(self ) :

tree = deepcopy(self ) tree.index () return tree

def get_matlab_equation(self ) : if self.parent is None:

if len(self.children) == 0:

return ' sum ({s}) '.format(s = TreeParameters.char_b [self.color])

else :

return TreeParameters.char_b[self.color] + '*((' + ').*('.j oin (

child.get_matlab_equation() for child in self.children ) + ')) '

else :

if len(self.children) == 0:

return TreeParameters.char_c[self.parent.color] +

else :

return TreeParameters.char_a[self.parent.color][ self.color] + '*((' + ').*('.join( ch.get_matlab_equation() for ch in self, children ) + ')) '

def get_tex_equation(self):

ch_a = tuple(filter(lambda x: x.get_order() > 1, self, children))

ch_c = tuple(filter(lambda x: x.get_order() == 1, self.

children)) h = self.get_height()

if self.parent is None:

s = '$' + tex_term( 'b', h, self.color) if len(ch_c) :

s += ' ' + tex_term ( ' c ' , h, self.color) + ' ~ ' + str(len(ch_c) )

120

125

130

135

140

if len(ch_a) == 1:

s += ch_a [0] .get_tex_equation() if len(ch_a) > 1:

s += '(' + ')\\cdot('.join(ch.get_tex_equation() for ch in ch_a) + ') '

return s + '= \\frac{l}{ ' + str(self.get_weight ()) + '}$ '

else :

s = tex_term( ' a', h, self.parent.color, self.color) if len(ch_c) :

s += ' ' + tex_term ( ' c ' , h, self.color) + ' ~ ' + str(len(ch_c) ) if len(ch_a) == 1:

s += ch_a[0] .get_tex_equation () if len(ch_a) > 1:

s += '(' + ').('.join(ch.get_tex_equation() for ch in ch_a) + ') ' return s

def get_leaf.count(self) : n — 0

for ch in self.children:

if ch.get_order() == 1: n += 1

return n

def get_order(self) :

return 1 + sum(child.get_order() for child in self, children)

def get_weight(self) :

w = self.get_order() for child in self.children : w *= child.get_weight() return w

def index(self, parent=Mone): self.parent = parent for child in self.children : child.index(self)

def get_height(self):

155

160

165

170

175

180

185

h = 1

p = self.parent while p is not None: h += 1

p = p•parent return h

def get_depth(self) : d = 1

if len(self.children) > 0:

child_depth = max(list(map(lambda x: x.get_depth(),

self.children))) if d < child_depth + 1: d = child_depth + 1

return d

def sort(self ) :

self.index(self.parent) for ch in self.children: ch.sort ()

self.children.sort(key=lambda x: x.get_matlab_equation() )

def grow(self) : forest = []

for color in TreeParameters.dependencies[self.color]: tree = self . clone ()

tree.children.append(Tree([color], self))

tree . sort ()

forest.append(tree)

for i, child in enumerate(self.children): if type(child) == Tree:

child_forest = child.grow() else :

child_forest = Tree([child]).grow()

for child_tree in child_forest: tree = self.clone()

tree.children[i] = child_tree.clone()

tree.sort ()

forest.append (tree)

200

205

210

map(Tree.index , forest)

return forest

def find_leaf (self , a, c) :

if len(self.children) == 0: return -1

if self.parent is not None:

if self.parent.color == a and self.color == c:

if self.get_height() > 1 and self.get_leaf.count () == len(self.children) <= TreeParameters. up_leaf [a] [c] : if not TreeParameters.strictly_triangular[ self.parent.color] [self.color] : return 1 if self.parent.parent is None: return 1

res = []

for child in self.children :

resl = child.find_leaf(a, c) res.append(resl)

res = max(res) return res

Листинг A.2 Скрипт для получения системы-следствия # -*- coding: utf-8 -*-

from tree import TreeParameters, Tree

def main () :

output = ( # 31ex ', 'matlab ' ,

)

10

forestO = {1: [Tree ( [0]) , Tree ( [1]) , Tree([2])]>

for order in range (2, TreeParameters.max.order + 1) forest = []

25

30

35

40

45

50

for tree in forestO[order-1]: forest += tree.growO

forestl = []

strings = [str(t) for t in forest] for i in range(len(forest)): is_chosen = 1

for j in range(i + 1, len (forest)) : if strings [i] == strings [j]: is_chosen = 0 break if is_chosen == 1:

forestl.append(forest[i])

forestO [order] = forestl

for order in ranged, Tr eeParamet ers . max.order + 1):

eq = [t.get_matlab_equation() for t in forestO [order]] forestl = []

for i in range(len(forestO [order])) : b = 1

for j in range(i + 1, len(forestO [order])) : if eq [i] == eq [ j] : b = 0 break if b == 1:

forestl.append(forestO [order] [i])

forestO [order] = forestl

full_system_size = 0 result_system_size = 0

for order in ranged, Tr eeParamet ers . max.order + 1): for tree in forestO [order] : tree.preamble = ' ' tree.up = ' '

# root presumptions

for a in TreeParameters.up_root[tree.color].keys() if TreeParameters.up_root[tree.color] [a] > 0: ch_a0 = 0

65

70

75

80

85

90

95

ch_al = 0

for child in tree.children:

if child.get_order() > 1: if child.color == a:

ch_a0 += 1 else :

ch_al += 1 if ch_al == 0 and \

ch_a0 == 1 and \ tree.get_leaf.count() < TreeParameters.up.root [ tree.color] [a] : tree.preamble = '

tree.up += '(y.n. B{c}{a}format(entree .color , a=a) +\

+ str(len(tree.children) -1) + ') '

# leaf presumptions

for a in TreeParameters.up.leaf.keys():

for c in TreeParameters.up.leaf.keys(): if TreeParameters.up.leaf[a][c] > 0: res = tree.find.leaf(a , c) if 0 <= res :

tree.preamble = '

tree.up += '(y.n. A{c}{a}format (c = tree . color , a=a) + '.1)'

for tree in forestO [order] : full.system.size += 1

if tree . preamble != ' °/0

result.system.size += 1 if 'tex ' in output :

print(tree.get.tex.equation () , '$(#', str ( result.system.size), ')$\\\\') if 'matlab' in output:

print ('F (end + 1) = {eq} - l/{w>; {num} '

.format(eq=tree.get_matlab_equation(), w=tree.get.weight(), num=full_system_size))

else

if 'tex ' in output :

105

pass

if 'matlab' in output:

# print ('% F(end+1) = {eq} - l/{w}; % {num}

{up} J

# .format (eq-tree.get_matIab_equation

O,

# w-tree . get_weight () ,

# num-full_system_size,

# up-tree.up)) pass

print(result_system_size, '/', full_system_size)

if __name_. main ()

.mam.

Приложение Б Текст программы ode46b

Листинг Б.1 Функция численного интегрирования систем класса B function varargout = ode46b (ode , tspan , yO , options , varargin)

solver.name = 'ode45';

X Stats nsteps = 0; nfailed = 0; nfevals = 0;

10

15

20

25

% Output

FcnHandlesUsed = isa(ode,'function_handle ' ) ; output_sol = (FcnHandlesUsed && (nargout==1));

sol = [] ; f 3d = [] ; if output_sol

sol.solver = solver.name;

sol.extdata.odefun = ode;

sol.extdata.options = options;

sol.extdata.varargin = varargin;

X Handle solver arguments for 1=1:4

[neq, tspan , ntspan, next, tO, tfinal, tdir , yO , fO , odeArgs , odeFcn , ...

options, threshold, rtol , normcontrol , normy , hmax , htry

, htspan , dataType] = ... odearguments(FcnHandlesUsed, solver.name, ode, tspan, yO , options , varargin);

nfevals = nfevals + 1;

% Handle the output if nargout > 0

outputFcn = odeget(options,'OutputFcn',[],'fast') ;

else

40

45

50

55

60

65

70

outputFcn = odeget(options , 'OutputFcnSodeplot ,'fast ');

outputArgs = {}; if isempty(outputFcn)

haveOutputFcn = false;

else

haveOutputFcn = true;

outputs = odeget(optionsl ,'OutputSel1:neq , 'fast ') ; if isa(outputFcn,'function_handle ')

X With MATLAB 6 syntax pass additional input arguments

to outputFcn. outputArgs = varargin;

end

refine = max (1 , odeget ( opt ions Ref ine 4 f ast ')) ; % Handle the event function

[haveEventFcn , eventFcn , eventArgs , valt , teout , yeout , ieout] = ... odeevents(FcnHandlesUsed,odeFcn,tO,yO,opt ions,varargin) ;

t = tO; y = yO;

% Allocate memory if we're generating output. nout = 0;

tout = [] ; yout = [] ; if nargout > 0

chunk = min(max (100 , 50*refine) , refine + floor((2~11)/neq)) ; tout = zeros(1,chunk,dataType) ; yout = zeros(neq,chunk,dataType); f3d = zeros(neq,7,chunk,dataType) ;

nout = 1;

tout(nout) = t ;

yout (: ,nout) = y ;

X Initialize method parameters. pow = 1/5;

[Bll B12 B21 B22 A1 A2 El E2] = coeff ; f(neq,7) = zeros(1,1,dataType);

hmin = 16* eps(t) ;

X Compute an initial step size h using y'(t). absh = min(hmax, htspan);

rh = norm(f0 ./ max(abs(y) ,threshold), inf) / (0.8 * rtol~pow); if absh * rh > 1 absh = 1 / rh;

end

85 absh = max (absh, hmin); f (: ,1) = f0;

structure = feval(odeFcn, 's'); 90 nyl = structure('1');

ny2 = structure('2')+nyl;

% Initialize the output function. if haveOutputFcn

f eval(outputFen ,[t tfinal],y( outputs) , 'init' ,outputArgs{:}) ;

end

% THE MAIN LOOP

100 done = false ; while "done

X By default, hmin is a small number such that t+hmin is

only slightly % different than t. It might be 0 if t is 0. 105 hmin = 16*eps(t);

absh = min(hmax, max(hmin, absh)); % couldn't limit absh

until new hmin h = tdir * absh;

% Stretch the step if within 10% of tfinal-t. if l.l*absh >= absCtfinal - t) h = tf inal - t ; absh = abs(h); done = true ;

end

% LOOP FOR ADVANCING ONE STEP.

125

130

135

140

145

150

nofailed = true ; while true

hAl = h * Al; hA2 = h * A2 ; hBll = h * Bll; hB12 = h * B12; hB21 = h * B21 ; hB22 = h * B22 ; for i-2 : 6

for j-1:nyl yt = [

y(

Z no failed attempts

1:ny1)+f( 1:nyl , :)*hBll(: ,i-1) y(nyl+1:ny2)+f(nyl+1:ny2,:)*hB12(:,i-1) ] ;

f(j,i) = feval (odeFcn , t + hAl(i-l), yt , j, odeArgs{ :}) ;

end

for j=nyl+l:ny2 yt = [

y( 1:ny1)+f( 1:nyl ,:)*hB21(: , i-1) y(nyl + 1:ny2)+f(nyl + 1:ny2 , :)*hB22 (: ,i-1) ] ;

f(j,i) = feval (odeFcn , t + hA2(i-l), yt , j, odeArgs{ :}) ;

end

end

tnew = t + hAl(6); if done