Разработка экономичных схем интегрирования структурно-разделенных систем обыкновенных дифференциальных уравнений

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

Разработка экономичных схем интегрирования структурно-разделенных систем обыкновенных дифференциальных уравнений

Работа посвящена конструированию структурных методов типа Рунге-Кутты, способных преодолевать барьер Бутчера для методов шестого порядка (классические методы должны иметь минимум 7 этапов). В явном виде представлено семипараметрическое семейство шестиэтапных методов шестого порядка. При конкретных значениях параметров построена вложенная расчётная схема, позволяющая оценивать локальную погрешность и автоматически регулировать шаг интегрирования. С помощью этой схемы функция MATLAB ode45 (пятый порядок точности) модифицирована для достижения шестого порядка при том же количестве этапов. Произведено сравнение на моделировании задач небесной механики и показаны преимущества построенной схемы.

Организация и управление

Дипломы

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 5be56f567966e104e4a123f7

UUID: 7d0f1cb0-c640-0136-3cd0-525400005860

Язык: Русский

Опубликовано: больше 5 лет назад

Просмотры: 8

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

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 684,3 КБ

Поделиться работой

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ Коврижных Николай Александрович Разработка экономичных схем интегрирования структурно разделённых систем обыкновенных дифференциальных уравнений Образовательная программа «Математическая кибернетика» Направление подготовки 02.06.01 — «Компьютерные и информационные науки» Выпускная квалификационная работа Научный руководитель: д. ф.-м. н., проф. Олемской Игорь Владимирович Рецензент: д. т. н., проф. Голоскоков Дмитрий Петрович Санкт-Петербург — 2018

2 Оглавление Стр. Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Обзор литературы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Глава 1. Структурные методы интегрирования . 1.1 Классы структурных особенностей систем ОДУ 1.1.1 Класс A . . . . . . . . . . . . . . . . . . . 1.1.2 Класс B . . . . . . . . . . . . . . . . . . . 1.1.3 Класс C . . . . . . . . . . . . . . . . . . . 1.1.4 Обзор моделей . . . . . . . . . . . . . . . 1.2 Метод интегрирования . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 8 8 9 10 11 12 . . . . . . . . . . . . . . . 16 16 18 . . . . . . . . . . 18 20 . . . . . . . . 22 22 23 24 Выводы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Заключение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 Список литературы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 Список рисунков . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 Глава 2.1 2.2 2.3 2. Схемы шестого порядка класса B . . . . . . . . . Условия порядка . . . . . . . . . . . . . . . . . . . . . . . Семейство схем класса B . . . . . . . . . . . . . . . . . . Интегрирование с переменным шагом. Оценка локальной погрешности . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Схема 𝑅𝐾𝐵6(4){7𝐹 } . . . . . . . . . . . . . . . . . . . . Глава 3.1 3.2 3.3 3. Численное исследование моделей Модель 1: точка либрации 𝐿1 . . . . . . . Модель 2: орбита Аренсторфа . . . . . . Результаты моделирования . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3 Стр. Список таблиц . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 Приложение А. Семейство методов шестого порядка класса B . 35 Приложение Б. Текст программы ode46b . . . . . . . . . . . . . . . 40

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

5 Начиная с 1970-х лет наступает время бурного развития новых методов интегрирования и способов сравнения одних методов с другими. Важными свойствами становятся не только порядок точности, но также возможность автоматического управления длиной шага и устойчивость методов. Эрвин Фель берг одним из первых предложил идею вложенных методов Рунге—Кутты, позволяющих существенно экономить вычисления с помощью автоматического управления длиной шага интегрирования [5]. Благодаря этому методы Рунге— Кутты обрели новую популярность, и вскоре на их основе были созданы схе мы Дж. Дорманда и П. Принса [6], одна из которых в качестве функции 𝑜𝑑𝑒45 сейчас является основным интегратором среды MATLAB.

6 Постановка задачи Барьеры Бутчера являются одним из основных недостатков методов Рун ге—Кутты. Способ их преодоления был предложен И. В. Олемским [7]. Он заключается в анализе структуры СОДУ, разбиении переменных на группы и модификации классической схемы Рунге—Кутты. Для многих задач, особенно в области механики, возможно сократить количество вычислений правой части системы, сохранив порядок точности. На основе структурного подхода было построено несколько вложенных схем, превосходящих по своим характеристи кам известные аналоги. К примеру, метод пятого порядка точности требует в полтора раза меньше обращений к правой части (4 вместо 6). Целью работы является построение экономичного метода интегрирова ния, преодолевающего барьер Бутчера для 𝑝 = 6 (требующего менее семи вычислений правой части СОДУ). Кроме того, метод должен сопровож даться алгоритмом оценки локальной погрешности на шаге. Подтверждение соответствия заявленному порядку производится путём сравнения с уже суще ствующими методами того же порядка на известных моделях.

7 Обзор литературы Несмотря на долгую историю, методы Рунге—Кутты по-прежнему при влекают много внимания. Популярна идея разбиения переменных на несколько групп с последующим применением к каждой группе отдельной вычислитель ной схемы. К примеру, ранние работы были связаны с разбиением жёстких задач на «жёсткую» и «нежёсткую» компоненты и построению явно-неявных методов [8]. Существуют работы, посвящённые разбиению СОДУ на «быстрые» и «медленные» компоненты для лучшего контроля динамических процессов в больших системах [9], а также получению схем с определёнными численными свойствами [10; 11]. Схожие идеи встречаются при решении систем уравнений в частных производных [12]. В работе [13] предложена идея непрерывных расширений структурных методов, дающих приближение к решению на всём интервале инте грирования, и продемонстрировано их применение для интегрирования систем с запаздыванием.

8 Глава 1. Структурные методы интегрирования 1.1 Классы структурных особенностей систем ОДУ Исторически первым использованием структуры системы обыкновенных дифференциальных уравнений для построения экономичных схем численного интегрирования можно считать методы Нюстрёма. В случае системы уравнений второго порядка 𝑦 ′′ = 𝑓 (𝑥,𝑦,𝑦 ′ ) алгоритм, предложенный Нюстрёмом, позволяет сократить машинную память, используемую при вычислениях. Однако насто ящее преимущество возникает в том случае, когда правая часть не зависит от первой производной 𝑦 ′ . Для систем вида 𝑦 ′′ = 𝑓 (𝑥,𝑦) (1.1) удаётся сократить также и количество вычислений функции 𝑓 , которые, как правило, занимают основную часть времени вычислительного процесса. Напри мер, в то время, когда классические методы Рунге—Кутты для достижения пятого порядка точности требуют минимум шести вычислений правой части, метод Нюстрёма требует только четырёх, что даёт экономию в полтора раза. Под термином «система ОДУ со структурными особенностями» обычно полагают то, что рассматриваемая система принадлежит к одному из трёх классов A, B или C. 1.1.1 Класс A К классу A относятся системы с перекрёстной структурой зависимостей: {︃ 𝑦1′ = 𝑓1 (𝑥,𝑦2 ), (1.2) 𝑦2′ = 𝑓2 (𝑥,𝑦1 ),

9 𝑥 ∈ [𝑋0 ,𝑋𝑘 ] ⊂ R, 𝑦𝑠 : [𝑋0 ,𝑋𝑘 ] −→ R𝑟𝑠 , 𝑠 = 1,2, 𝑓1 : [𝑋0 ,𝑋𝑘 ] × R𝑟2 −→ R𝑟1 , 𝑓2 : [𝑋0 ,𝑋𝑘 ] × R𝑟1 −→ R𝑟2 . Можно видеть, что этот класс является обобщением систем второго по рядка, не зависящих от первой производной. Каждая система вида 𝑦 ′′ = 𝑓 (𝑥,𝑦) с помощью введения обозначений 𝑦1 = 𝑦, 𝑦2 = 𝑦 ′ может быть записана в виде системы первого порядка с перекрёстной системой связей {︃ 𝑦1′ = 𝑦2 , (1.3) 𝑦2′ = 𝑓2 (𝑥,𝑦1 ). Важным качеством при этом является то, что, в отличие от систем второго порядка, две группы переменных могут иметь разную размерность. 1.1.2 Класс B К классу B относятся системы со следующей структурой зависимостей: {︃ 𝑦𝑖′ = 𝑓𝑖 (𝑥,𝑦1 , . . . ,𝑦𝑖−1 ,𝑦𝑙+1 , . . . ,𝑦𝑛 ), 𝑖 = 1,2, . . . ,𝑙, (1.4) 𝑦𝑗′ = 𝑓𝑗 (𝑥,𝑦1 , . . . ,𝑦𝑗−1 ), 𝑗 = 𝑙 + 1, . . . ,𝑛, ρ= 𝑛 ∑︁ 𝑟ν , 𝑥 ∈ [𝑋0 ,𝑋𝑘 ] ⊂R, 𝑦𝑠 :[𝑋0 ,𝑋𝑘 ] −→ R𝑟𝑠 , 𝑠 = 1, . . . ,𝑛, ν=1 𝑓𝑖 : [𝑋0 ,𝑋𝑘 ] × R ρ−^ 𝑟𝑖 𝑗 𝑟𝑖 −→R , 𝑓𝑗 : [𝑋0 ,𝑋𝑘 ] × Rρ−¯𝑟 −→R𝑟𝑗 , 𝑖 𝑟^ = 𝑟¯𝑗 = 𝑙 ∑︁ ν=𝑖 𝑛 ∑︁ 𝑟ν , 𝑖 = 1,2, . . . ,𝑙, 𝑟ν , 𝑗 = 𝑙 + 1, . . . ,𝑛. ν=𝑗 Этот класс является обобщением класса A: в системах такого вида неиз вестные двух групп могут зависеть не только от противоположной группы, но и от неизвестных своей группы, имеющих меньший индекс.

10 1.1.3 Класс C Следующее обобщение — класс C, к нему относятся системы с т.н. «об щей» группой 𝑦0 : ⎧ ⎪ 𝑦 ′ = 𝑓0 (𝑥,𝑦0 , . . . ,𝑦𝑛 ), ⎪ ⎨ 0 𝑦𝑖′ = 𝑓𝑖 (𝑥,𝑦0 , . . . ,𝑦𝑖−1 ,𝑦𝑙+1 , . . . ,𝑦𝑛 ), ⎪ ⎪ ⎩𝑦 ′ = 𝑓 (𝑥,𝑦 , . . . ,𝑦 ), 𝑗 0 𝑗−1 𝑗 ρ= 𝑛 ∑︁ 𝑟ν , (1.5) 𝑖 = 1,2, . . . ,𝑙, (1.6) 𝑗 = 𝑙 + 1, . . . ,𝑛, (1.7) 𝑦𝑠 :[𝑋0 ,𝑋𝑘 ] −→ R𝑟𝑠 , 𝑥 ∈ [𝑋0 ,𝑋𝑘 ] ⊂R, 𝑠 = 0, . . . ,𝑛, ν=0 𝑓0 : [𝑋0 ,𝑋𝑘 ] × Rρ −→R𝑟0 , 𝑓𝑖 : [𝑋0 ,𝑋𝑘 ] × R ρ−^ 𝑟𝑖 𝑗 𝑟𝑖 −→R , 𝑓𝑗 : [𝑋0 ,𝑋𝑘 ] × Rρ−¯𝑟 −→R𝑟𝑗 , 𝑖 𝑟^ = 𝑟¯𝑗 = 𝑙 ∑︁ ν=𝑖 𝑛 ∑︁ 𝑟ν , 𝑖 = 1,2, . . . ,𝑙, 𝑟ν , 𝑗 = 𝑙 + 1, . . . ,𝑛. ν=𝑗 Последовательная вложенность классов A, B и C позволяет использовать одни и те же названия для групп переменных 𝑦𝑘 : здесь и далее мы будем на зывать их нулевой (1.5), первой (1.6) и второй (1.7). В литературе также часто называют нулевую группу общей, а первую и вторую обозначают как 𝑖-ю и 𝑗-ю. Нулевая группа — единственная, в которой производная какого-либо па раметра системы может зависеть от него самого. В случаях классов B и C любая группа может быть вырожденной: на пример, система вида {︃ ′ 𝑦0 = 𝑓0 (𝑡,𝑦0 , 𝑦1 ), 𝑦1′ = 𝑓0 (𝑡,𝑦0 ) принадлежит классу C (группа 2 отсутствует).

11 𝐸 а) б) в) г) а) методы Нюстрёма (𝐸 обозначает единичную матрицу); б) класс A; в) класс B; г) класс C Рисунок 1.1 — Структуры допустимых зависимостей в СОДУ 1.1.4 Обзор моделей Поведение многих математических моделей описывается дифференциаль ными уравнениями второго порядка вида вида 1.1. Как уже было показано, каждая такая система может быть приведена к перекрёстному виду 1.2 (класс A) с помощью введения новых переменных. Существуют также и модели, опи сываемые системами класса A, но не являющиеся системами второго порядка. Примером такой задачи может являться движение неуправляемого космическо го аппарата (КА) вокрут точек либрации 𝐿1 и 𝐿2 [14]. Более сложные задачи относительного движения нескольких тел также часто могут быть приведены к виду 1.4 (класс B). Например: 1. Модель Хилла—Клохесси—Уилтшира: движение объектов в окрестно сти искусственного спутника Земли [15]. 2. Модель Чаунера—Хемпеля: движение двух спутников по эллиптиче ской орбите [16]. 3. Модель Швайгарта—Седвика: движение КА в геопотенциальном по ле [17]. 4. Орбита Аренсторфа: периодическое движение КА в гравитационном поле Земли и Луны [18]. Все эти модели описываются системами второго порядка, однако при соотвест вующем переопределении переменных (менее очевидном, чем описанное выше для класса A) могут быть отнесены к классу B. Пример такого переопределе ния для орбиты Аренсторфа будет рассмотрен ниже при проведении численных экспериментов.

12 Очевидно также, что любая система ОДУ, описывающая какую-либо модель, может быть представлена как система класса C, содержащая лишь ну левую группу. Однако представление системы как имеющей только нулевую группу не принесёт результата по сравнению с классическими методами Рун ге—Кутты, поскольку преимущество структурных методов заключается именно в экономии на вычислениях правых частей первой и второй группы. Многие встречающиеся в реальной жизни системы имеют структуру, не позволяющую сразу выделить первую или вторую группу уравнений, однако с помощью простого переопределения переменных удаётся привести их к необхо димому виду. К примеру, если зависимости производных представляют собой двудольный граф, она является системой класса A. Существует также алго ритм, позволяющий переопределять параметры системы таким способом, чтобы первая и вторая группа имели максимальный возможный размер [19]. 1.2 Метод интегрирования Метод интегрирования систем класса C. Будем считать, что нам из вестно точное решение 𝑦𝑠 (𝑥), 𝑠 = 0, . . . ,𝑛, системы 1.5—1.7 в точке 𝑥 ∈ [𝑋0 ,𝑋𝑘 ]. Не умаляя общности рассуждений, для простоты вывода будем считать 𝑟𝑠 = 1, 𝑠 = 0, . . . ,𝑛. Для численного интегрирования систем 1.5—1.7 рассматривается явный одношаговый метод [20]. В предположении достаточной гладкости правой части рассматриваемой системы приближение 𝑦˜𝑠 к точному решению 𝑦𝑠 (𝑥 + ℎ), 𝑠 = 1, . . . ,𝑛 в точке 𝑥 + ℎ ∈ [𝑋0 ,𝑋𝑘 ] ищется в виде: 𝑦0 (𝑥 + ℎ) ≈ 𝑦˜(𝑥 + ℎ) = 𝑦0 (𝑥) + ℎ 𝑦𝑖 (𝑥 + ℎ) ≈ 𝑦˜𝑖 (𝑥 + ℎ) = 𝑦𝑖 (𝑥) + ℎ 𝑚0 ∑︁ 𝑏0ν 𝑘0,ν (ℎ), ν=1 𝑚1 ∑︁ 𝑦𝑗 (𝑥 + ℎ) ≈ 𝑦˜𝑗 (𝑥 + ℎ) = 𝑦𝑗 (𝑥) + ℎ 𝑏1ν 𝑘𝑖,ν (ℎ), ν=1 𝑚1 ∑︁ ν=1 𝑏2ν 𝑘𝑗,ν (ℎ), (1.8) 𝑖 = 1, . . . ,𝑙, 𝑗 = 𝑙 + 1, . . . ,𝑛, (1.9) (1.10)

13 причём 𝑘𝑠,𝑤 ≡ 𝑘𝑠,𝑤 (ℎ) вычисляются в строгой последовательности 𝑘0,1 , . . . , 𝑘𝑛,1 , 𝑘0,2 , . . . , 𝑘𝑛,2 , 𝑘0,3 , 𝑘1,3 , . . . (1.11) по формулам (︁ 𝑘0,ν = 𝑓0 𝑥 + 𝑐0ν ℎ, 𝑦0 + ℎ ν−1 ∑︁ 𝑎00νµ 𝑘0,µ (ℎ), µ=1 𝑦1 + ℎ 𝑦𝑙+1 + ℎ ν−1 ∑︁ 𝑎01νµ 𝑘1,µ (ℎ), . . . , 𝑦𝑖−1 + ℎ ν−1 ∑︁ µ=1 µ=1 ν−1 ∑︁ ν−1 ∑︁ 𝑎02νµ 𝑘𝑙+1,µ (ℎ), . . . , 𝑦𝑛 + ℎ µ=1 𝑎01νµ 𝑘𝑙,µ (ℎ), )︁ 𝑎02νµ 𝑘𝑛,µ (ℎ) , (1.12) µ=1 (︁ 𝑘𝑖,ν = 𝑓𝑖 𝑥 + 𝑐1ν ℎ, 𝑦0 + ℎ 𝑦1 + ℎ 𝑦𝑙+1 + ℎ ν ∑︁ µ=1 ν ∑︁ 𝑎10νµ 𝑘0,µ (ℎ), 𝑎11νµ 𝑘1,µ (ℎ), . . . , 𝑦𝑖−1 + ℎ ν ∑︁ µ=1 µ=1 ν−1 ∑︁ ν−1 ∑︁ 𝑎12νµ 𝑘𝑙+1,µ (ℎ), . . . , 𝑦𝑛 + ℎ µ=1 𝑎11νµ 𝑘𝑖−1,µ (ℎ), )︁ 𝑎12νµ 𝑘𝑛,µ (ℎ) , 𝑖 = 1, . . . ,𝑙, µ=1 (1.13) (︁ 𝑘𝑗,ν = 𝑓𝑗 𝑥 + 𝑐2ν ℎ, 𝑦0 + ℎ 𝑦1 + ℎ 𝑦𝑙+1 + ℎ ν ∑︁ µ=1 ν ∑︁ µ=1 ν ∑︁ µ=1 𝑎20νµ 𝑘0,µ (ℎ), 𝑎21νµ 𝑘1,µ (ℎ), . . . , 𝑦𝑖−1 + ℎ 𝑎22νµ 𝑘𝑙+1,µ (ℎ), . . . , 𝑦𝑛 + ℎ ν ∑︁ µ=1 ν ∑︁ 𝑎21νµ 𝑘𝑙,µ (ℎ), )︁ 𝑎22νµ 𝑘𝑗−1,µ (ℎ) , 𝑗 = 𝑙 + 1, . . . ,𝑛. µ=1 (1.14) Параметрами метода являются коэффициенты 𝑎𝑤1 𝑤2 νµ (матрицы этих ко эффициентов будем обозначать как 𝐴𝑤1 𝑤2 ), 𝑏𝑤1 ν , 𝑐𝑤1 ν , 𝑤1 , 𝑤2 ∈ {0,1,2}. Говорят,

14 что метод имеет порядок 𝑝, если |𝑦𝑠 (𝑥 + ℎ) − 𝑦˜𝑠 | ≈ 𝑂(ℎ𝑝+1 ), 𝑠 = 0, . . . ,𝑛. Строгий порядок вычисления значений 𝑘 связывает количество этапов по каждой группе: 𝑚1 6 𝑚0 6 𝑚1 + 1. Основная идея структурного подхо да заключается в алгоритмическом использовании структурных особенностей математической модели для сокращения числа обращений к правым частям первой и второй групп уравнений, поэтому мы будем считать далее 𝑚 = 𝑚0 = 𝑚1 + 1. Для того, чтобы метод был явным, необходимо, чтобы матрицы 𝐴10 , 𝐴20 , 𝐴11 , 𝐴21 и 𝐴22 были нижнетреугольными, а матрицы 𝐴00 , 𝐴01 , 𝐴02 и 𝐴12 — строго нижнетреугольными. В этом случае к моменту вычисления каждого из выражений 1.12—1.14 аргументы в их правых частях будут уже известны. В табл. 1 приведён компактный общий вид коэффициентов явного метода при 𝑚 = 𝑚0 = 𝑚1 + 1 (для удобства индексы отделены запятыми). Дополнительное требование 𝑎1111 = 𝑎1011 = 0 обусловлено часто используемым предположением ∑︀ 𝑐𝑤1 ν = 𝑎𝑤1 𝑤2 νµ , 𝑤1 ,𝑤2 ∈ {0,1,2} и тем, что 𝑎1211 = 0. µ Замечание 1. В случае, когда система состоит только из нулевой груп пы, метод вырождается в классический метод Рунге—Кутты с коэффициентами 𝐴 = 𝐴00 , 𝑏 = 𝑏0 , 𝑐 = 𝑐0 и количеством этапов 𝑠 = 𝑚0 : 𝑦(𝑥 + ℎ) ≈ 𝑦˜(𝑥 + ℎ) = 𝑦(𝑥) + ℎ (︃ 𝑘ν = 𝑓 𝑥 + 𝑐ν ℎ, 𝑦 + ℎ 𝑠 ∑︁ ν=1 𝑠 ∑︁ 𝑏ν 𝑘ν (ℎ), (1.15) )︃ 𝑎νµ 𝑘µ (ℎ) , µ=1 Замечание 2. В силу того, что системы класса B отличаются от систем класса C лишь отсутствием нулевой группы, рассмотренный метод можно при менить и к ним, просто исключив из всех формул параметры с индексом 0. То же касается и систем класса A: достаточно вдобавок исключить все параметры, использующие зависимости внутри первой группы и внутри второй. Табличные записи коэффициентов для обеих схем будут содержать коэффициенты 𝑏1ν , 𝑏2ν , 𝑐1ν , 𝑐2ν , 𝑎12νµ , 𝑎21νµ , а для схемы класса B дополнительно 𝑎11νµ и 𝑎22νµ . Требо вания к треугольности матриц для явных методов останутся теми же.

Таблица 1 — Коэффициенты явного метода класса C при 𝑚 = 𝑚0 = 𝑚1 + 1 𝑐0,ν 𝑎0,0,ν,µ 0 0 𝑐0,2 .. . 𝑐0,𝑚−1 𝑐0,𝑚 𝑎0,1,ν,µ 𝑎0,0,2,1 .. . 0 0 .. . ... ... .. . 0 0 .. . 𝑎0,0,𝑚−1,1 𝑎0,0,𝑚,1 𝑎0,0,𝑚−1,2 𝑎0,0,𝑚,2 ... ... 0 𝑐1,ν 𝑎0,0,𝑚,𝑚−1 0 𝑎0,1,2,1 .. . 0 0 .. . ... ... .. . 0 0 .. . 𝑎0,1,𝑚−1,1 𝑎0,1,𝑚,1 𝑎0,1,𝑚−1,2 𝑎0,1,𝑚,2 ... ... 0 𝑎1,0,ν,µ 0 0 0 𝑐1,2 .. . 𝑎1,0,2,1 .. . 𝑐1,𝑚−2 𝑐1,𝑚−1 𝑎1,0,𝑚−2,1 𝑎1,0,𝑚−1,1 𝑎1,0,2,2 .. . 0 0 .. . 𝑎1,0,𝑚−2,2 𝑎1,0,𝑚−1,2 ... ... 0 𝑎1,0,𝑚−1,𝑚−1 0 0 𝑎1,1,2,1 .. . 𝑎1,1,𝑚−2,1 𝑎1,1,𝑚−1,1 𝑐2,𝑚−2 𝑐2,𝑚−1 𝑎2,0,𝑚−2,1 𝑎2,0,𝑚−1,1 0 𝑎0,2,2,1 .. . 0 0 .. . ... ... .. . 0 0 .. . 𝑏0,1 𝑏0,2 .. . 𝑎0,2,𝑚−1,1 𝑎0,2,𝑚,1 𝑎0,2,𝑚−1,2 𝑎0,2,𝑚,2 ... ... 0 𝑏0,𝑚−1 𝑏0,𝑚 𝑎1,1,2,2 .. . 0 0 .. . 𝑎1,1,𝑚−2,2 𝑎1,1,𝑚−1,2 ... ... 0 𝑎1,1,𝑚−1,𝑚−1 0 𝑎2,0,2,2 .. . 0 0 .. . 𝑎2,1,1,1 𝑎2,1,2,1 .. . 𝑎2,0,𝑚−2,2 𝑎2,0,𝑚−1,2 ... ... 0 𝑎2,1,𝑚−2,1 𝑎2,1,𝑚−1,1 𝑎2,0,𝑚−1,𝑚−1 0 𝑏1,ν 𝑎1,2,2,1 .. . 0 0 .. . ... ... .. . 0 0 .. . 𝑏1,1 𝑏1,2 .. . 𝑎1,2,𝑚−2,1 𝑎1,2,𝑚−1,1 𝑎1,2,𝑚−2,2 𝑎1,2,𝑚−1,2 ... ... 0 0 𝑏1,𝑚−2 𝑏1,𝑚−1 𝑎2,1,ν,µ ... ... .. . 𝑎0,2,𝑚,𝑚−1 𝑎2,2,ν,µ 𝑎2,1,2,2 .. . ... ... .. . 0 0 .. . 𝑎2,2,1,1 𝑎2,2,2,1 .. . 𝑎2,1,𝑚−2,2 𝑎2,1,𝑚−1,2 ... ... 0 𝑎2,2,𝑚−2,1 𝑎2,2,𝑚−1,1 𝑎2,1,𝑚−1,𝑚−1 0 𝑏2,ν 𝑎2,2,2,2 .. . ... ... .. . 0 0 .. . 𝑏2,1 𝑏2,2 .. . 𝑎2,2,𝑚−2,2 𝑎2,2,𝑚−1,2 ... ... 0 𝑏2,𝑚−2 𝑏2,𝑚−1 𝑎2,2,𝑚−1,𝑚−1 15 𝑎2,0,1,1 𝑎2,0,2,1 .. . 𝑏0,ν 𝑎1,2,ν,µ ... ... .. . 𝑎2,0,ν,µ 𝑐2,1 𝑐2,2 .. . 𝑎0,1,𝑚,𝑚−1 0 𝑎1,1,ν,µ ... ... .. . 𝑐2,ν 𝑎0,2,ν,µ

16 Глава 2. Схемы шестого порядка класса B 2.1 Условия порядка Коэффициенты схемы интегрирования типа Рунге—Кутты должны удо влетворять системе условий порядка. Эти условия представляют собой ал гебраические уравнения, получающиеся при приравнивании соответствующих членов ряда Тейлора точного 𝑦 и приближённого 𝑦˜ решений вплоть до 𝑝-го, где 𝑝 — порядок метода. Размеры системы условий порядка для классических и структурных методов типа Рунге—Кутты приведены в таблице 2. Таблица 2 — Размеры системы условий порядка количество условий порядка класс метода 𝑅𝐾 1 2 4 8 17 37 A 2 4 8 16 34 74 B 2 4 10 28 88 292 C 3 6 18 66 276 1224 порядок метода 1 2 3 4 5 6 Быстрый рост количества уравнений с ростом порядка вынуждает авто матизировать процесс их получения. Джон Бутчер развил теорию помеченных деревьев, с помощью которой можно однозначно сопоставить каждый элемен тарный дифференциал Φ некоему древовидному графу τ(Φ). В классической теории для получения полного набора условий требуется выписать все поме ченные деревья порядка не выше 𝑝 (порядком дерева называется количество его вершин) и соответствующие им уравнения. Примеры соответствия деревьев и уравнений приведены в таблице 3. В работе [21] предложена модификация этой теории, предназначенная для методов класса C и представлены условия порядка до пятого включительно. В случае структурных методов класса C помимо получения полного набора де ревьев необходимо сопоставить индекс 0, 1 или 2 каждой вершине, имеющей потомков. Большим количеством возможных расстановок индексов и объясня ется возрастание размера системы при переходе к структурным методам.

17 Таблица 3 — Примеры построения условий порядка τ γ(τ) условие порядка ∑︀ 𝑏0ν = 1 0 1 ν 1 0 1 1 2 1 0 1 𝑏0ν 𝑐20ν = 1 3 ∑︀ 6 ∑︀ ν µ 5 ∑︀ 𝑏2ν 𝑐42ν = 36 ∑︀ 24 ∑︀ 48 ∑︀ 120 ∑︀ 2 0 1 2 3 2 2 𝑏1ν 𝑐1ν = 2 ∑︀ ν ν ν ν ν 𝑏1ν ∑︀ 𝑎11νµ 𝑐1µ = 1 6 1 5 ∑︀ ∑︀ 𝑏2ν ( 𝑎20νµ 𝑐20µ ) · ( 𝑎21νµ 𝑐1µ ) = µ µ 𝑏2ν 𝑐12ν ∑︀ 𝑏0ν 𝑐0ν ∑︀ ν µ 𝑎22νµ 𝑐32µ = 𝑎01νµ 𝑐1µ µ ∑︀ 1 24 𝑎10µξ 𝑐0ξ = ξ 0 1 1 2 ν 𝑏2ν ∑︀ µ 𝑎21νµ ∑︀ ξ 1 36 𝑎11µξ 𝑐31ξ = 1 120 1 48

18 Отметим, что систему условий для методов класса B того же порядка можно получить из неё, отбросив все уравнения, соответствующие деревьям, содержащим вершины с индексом 0 (поскольку используются только индек сы 1 и 2). Аналогично, систему условий порядка для методов класса A можно по лучить из последней отбрасыванием деревьев, имеющих рёбра вида 1 − 1 и 2 − 2. Таким образом останутся только те деревья, рёбра которых соединяют вершины с разными индексами. 2.2 Семейство схем класса B Система условий шестого порядка для методов класса B состоит из 292 алгебраических уравнений. Было получено семипараметрическое семейство ре шений этой системы со свободными параметрами 𝑐12 = α1 𝑐13 = α2 𝑐14 = α3 𝑐23 = α4 𝑐24 = α5 𝑎1132 = α6 и 𝑎2243 = α7 . В прил. А представлены значения ненулевых параметров метода. 2.3 Интегрирование с переменным шагом. Оценка локальной погрешности Регулирование длины шага интегрирования является основным инстру ментом при численном решении задачи Коши. Малая длина шага, как правило, даёт большую точность, однако приводит и к большим вычислительным затра там. Неоправданно высокие вычислительные затраты в свою очередь могут привести к тому, что основной составляющей глобальной погрешности станет накопившаяся ошибка округления. Кроме того, во многих моделях интервал интегрирования можно раз делить на «быстрые» и «медленные» участки: в первом случае переменные меняются «быстро», и для удержания локальной погрешности в заданных рамках приходится делать малые шаги, во втором — «медленно», и можно уве личить длину шага, чтобы ускорить вычислительный процесс.

19 Вложенные методы. Обозначим результат применения метода Рунге— Кутты 1.15 на одном шаге длины ℎ как 𝑦˜ = 𝑅𝐾(𝑥0 , 𝑦0 , 𝑓, 𝐴, 𝑏, ℎ), (2.1) где 𝑦0 и 𝑓 — параметры задачи Коши, 𝐴 и 𝑏 — коэффициенты конкретной схемы Рунге—Кутты, ℎ — длина шага итегрирования. Главный параметр, в соответствии с которым изменяют длину шага — это локальная погрешность, которая на каждом шаге не должна превышать некоей заданной величины 𝑡𝑜𝑙. Э. Фельберг предложил оценивать локальную погреш ность, вычисляя одновременно два приближения: 𝑦˜1 = 𝑅𝐾(𝑥0 , 𝑦0 , 𝑓, 𝐴, 𝑏, ℎ) порядка 𝑝 и 𝑦˜2 = 𝑅𝐾(𝑥0 , 𝑦0 , 𝑓, 𝐴, 𝑑, ℎ) порядка 𝑞 < 𝑝, то есть с помощью двух схем Рунге—Кутты, отличающихся лишь весовыми коэффициентами 𝑏 и 𝑑. Главные члены локальных погрешностей для этих приближений будут соответ ственны равны Δ1 = 𝐶1 ℎ𝑝+1 и Δ2 = 𝐶2 ℎ𝑞+1 , где 𝐶1 и 𝐶2 — некоторые скаляры. Матрицы коэффициентов 𝐴 для обеих схем одинаковы, поэтому оба при ближения можно получить, не затрачивая лишнее время на вычисления правых частей СОДУ. Подобные пары получили название вложенных методов. По скольку 𝑦˜1 имеет больший порядок, главный член их разности Δ = 𝑦˜1 − 𝑦˜2 будет равен 𝐶2 ℎ𝑞+1 и может быть использован в качестве оценки локальной погреш ности приближённого решения 𝑦˜2 . В дальнейших работах (Дорманд, Принс) в качестве приближения стали использовать 𝑦˜1 , поскольку его порядок выше, а величина Δ приобрела значение уже не оценки локальной погрешности, а некоего контрольного члена, применяемого для регулирования длины шага. Другим отличительным свойством вложенных методов типа Дорманда— Принса является их свойство FSAL (First Same As Last). 𝑠-этапную схему инте ∑︀ грирования дополняют ещё одним этапом 𝑘𝑠+1 = 𝑓 (𝑥 + ℎ, 𝑦0 + ℎ 𝑠µ=1 𝑏µ 𝑘µ (ℎ)). Фактически 𝑘𝑠+1 совпадает с результатом первого вычисления правой части, которое будет сделано на следующем шаге. А значит, на всех шагах начиная со второго можно по-прежнему совершать 𝑠 вычислений правой части, и вычис лительные затраты почти не растут.

20 2.4 Схема 𝑅𝐾𝐵6(4){7𝐹 } На основе решения системы условий порядка, полученного при значени ях свободных параметров α1 = 29 , α2 = 61 , α3 = 12 , α4 = 61 , α5 = 21 , α6 = 0, α7 = 23 построена семиэтапная расчетная схема метода 𝑅𝐾𝐵6(4){7𝐹 }. Это название означает: – 𝑅𝐾 — метод типа Рунге—Кутты; – 𝐵 — метод класса B; – 6 — порядок точности приближения на шаге; – 4 — порядок точности вложенного метода—оценщика; – 7 — количество используемых значений правой части СОДУ на шаге; – 𝐹 — метод обладает свойством FSAL, фактически на каждом шаге на чиная со второго требуется только шесть обращений к правой части СОДУ. Значения параметров были выбраны из двух соображений: минимизация невяз ки по условиям седьмого порядка (это позволяет уменьшить норму главного члена локальной погрешности [22]) и краткость записи. Коэффициенты метода представлены в таблице 4. Замечание. Благодаря тому, что все коэффициенты 𝑎𝑤1 𝑤2 1µ равны нулю, построенный метод избегает существенного недостатка, которым могут обла дать другие вложенные структурные методы со свойством FSAL. В случае, когда среди коэффициентов 𝑎𝑤1 𝑤2 1µ есть ненулевые, нельзя утверждать, что в выражениях 1.12—1.14 обращения к правым частям СОДУ на смежных шагах происходят с одинаковыми значениями аргументов. Это приводит к необходимо сти модифицировать алгоритм адаптивного выбора шага интегрирования [23].

21 𝑐1ν 0 2 9 1 6 1 2 5 6 1 1 Таблица 4 — Коэффициенты метода 𝑅𝐾𝐵6(4){7𝐹 } 𝑎11νµ 𝑎12νµ 7 150 1 9 1 12 −1 44 7 36 −3 7 7 150 1 9 0 0 0 0 0 𝑐2ν 0 2 9 1 6 1 2 5 6 1 1 𝑏1ν 𝑑1ν 1 12 9 22 0 9 8 27 100 5 44 5 9 −5 28 11 30 1 12 27 56 27 100 7 150 2 9 5 48 37 176 −635 432 29 4 7 150 1 16 243 176 −167 16 1377 28 0 𝑎21νµ 1 9 7 48 −31 176 73 144 −39 28 7 150 1 9 3 16 −81 176 15 16 −81 28 0 −1 6 45 44 −5 4 279 56 27 100 5 44 5 9 −5 28 11 30 −12 11 100 9 −1425 28 27 100 44 27 −11 2 11 30 27 28 27 100 7 150 𝑎22νµ 1 12 27 56 27 100 7 150 1 9 7 48 −185 1584 1031 3888 −29 63 7 150 1 9 3 16 −123 880 −53 144 15 7 0 −1 6 2 3 65 324 −103 168 27 100 13 200 0 0 27 100 11 30 27 100 7 150 183 800 33 80 183 800 7 300 1 24 0 𝑏2ν 𝑑2ν 89 990 317 486 −139 252 11 30 1 12 27 56 27 100 7 150 7 150 13 200 0 0 27 100 11 30 27 100 7 150 183 800 33 80 183 800 7 300 1 24 0

22 Глава 3. Численное исследование моделей Для проверки качества построенного метода было произведено сравнение с явными одношаговыми методами шестого порядка Цитураса [24], Вернера [25] и Эль-Миккави [26] а также методом Дорманда—Принса пятого порядка [27]. Помимо набора коэффициентов схема интегрирования должна также об ладать алгоритмом выбора длины шага. Для объективной оценки качества работы конкретных методов за основу был взят алгоритм из реализации метода Дорманда—Принса в среде MATLAB — функции 𝑜𝑑𝑒45, основного интегратора в этой среде [27]. Таким образом, различия между методами заключались лишь в затратах на вычисления правых частей СОДУ и достигаемой точности на ша ге. На основе метода 𝑅𝐾𝐵6(4){7𝐹 } реализована модификация 𝑜𝑑𝑒46𝑏 той же функции, позволяющая обращаться к группам системы независимо. Аналогич ная модификация функций MATLAB и их внедрение для методов класса A были представлены в работе [28]. Текст программы 𝑜𝑑𝑒46𝑏 приведён в прил. Б. Интерфейс функции 𝑜𝑑𝑒46𝑏 аналогичен интерфейсу 𝑜𝑑𝑒45. Дополнитель но от функции вычисления значений правой части СОДУ (аргумент под именем 𝑜𝑑𝑒) требуется возможность вернуть информацию о структуре СОДУ (количе ство уравнений первой и второй групп) и выдавать значение правой части для конкретных уравнений СОДУ. Для исследования были выбраны модели, описывающие движение косми ческого аппарата в окрестности точки либрации 𝐿1 системы Солнце—Земля и движение КА по орбите Аренсторфа. 3.1 Модель 1: точка либрации 𝐿1 Плоское движение неуправляемого КА во вращающейся системе коор динат в окрестности точки либрации 𝐿1 = (1,1) системы Солнце–Земля

23 описывается [14] системой обыкновенных дифференциальных уравнений: ⎧ ⎪ 𝑥˙ 1 = 𝑥2 + 𝑦1 , ⎪ ⎪ ⎪ ⎪ ⎨𝑥˙ = −𝑥 + 𝑦 , 2 1 2 (3.1) ⎪ 𝑦 ˙ = 8(𝑥 − 1) + (𝑦 − 1), ⎪ 1 1 2 ⎪ ⎪ ⎪ ⎩ 𝑦˙ = −4𝑥 − 𝑦 . 2 2 1 При начальных условиях √ 7−3 ε,0), 2 𝑦(0) = (0,1 + ε) 𝑥(0) = (1 + КА движется по периодической орбите вокруг 𝐿1 . При переопределении пара метров: 𝑥˜1 = 𝑥1 , 𝑥˜2 = 𝑦2 , 𝑥˜3 = 𝑥2 , 𝑥˜4 = 𝑦1 СОДУ приобретает перекрестную 1 . структуру. Для тестирования было выбрано значение ε = 100 Поскольку полученная система является линейной, существует возмож ность на каждом шаге сравнивать численное решение с точным. 3.2 Модель 2: орбита Аренсторфа Плоское движение космического аппарата с координатами (𝑥1 ,𝑥2 ) в гра витационном поле, создаваемом Землей (0,0) и Луной (1,0) описывается [18] системой обыкновенных дифференциальных уравнений: ⎧ 𝑥 +µ 𝑥1 − µ′ ⎪ ⎨𝑥¨1 = 𝑥1 + 2𝑥˙ 2 − µ′ 1 −µ , 𝐷1 𝐷2 (3.2) 𝑥 𝑥 2 2 ⎪ ′ ⎩𝑥¨2 = 𝑥2 − 2𝑥˙ 1 − µ −µ , 𝐷1 𝐷2 где (︁ )︁3/2 (︁ )︁3/2 2 2 ′ 2 2 𝐷1 = (𝑥1 + µ) + 𝑥2 , 𝐷2 = (𝑥1 − µ ) + 𝑥2 , µ = 0.012277471, µ′ = 1 − µ. При начальных условиях 𝑥(0) = (0.994, 0), 𝑥(0) ˙ = (0, 2.00158510637908252240537862224)

24 КА движется по орбите с периодом 𝑇𝑝𝑒𝑟 = 17.0652165601579625588917206249. При переопределении параметров: 𝑥˜1 = 𝑥1 , 𝑥˜2 = 𝑥˙ 2 , 𝑥˜3 = 𝑥2 , 𝑥˜4 = 𝑥˙ 1 система приобретает структуру, описываемую классом B. Поскольку орбита периодична, глобальной погрешностью численного при ближения считается отклонение от исходного положения. 3.3 Результаты моделирования На рис. 3.1 и 3.2 приведён график зависимости глобальной погрешности 𝐸𝑟𝑟 от трудоёмкости 𝑁𝑓 𝑒𝑣𝑎𝑙𝑠 (количества обращений к правой части СОДУ). 13 12 11 -log10 Err 10 9 8 ode45 ode56tsit1999 ode56vern1994 ode56elm2003 ode46b 7 6 5 1.6 1.8 2 2.2 2.4 2.6 2.8 3 log10 Nfevals Рисунок 3.1 — Зависимость глобальной погрешности от количества вычислений для модели 1

25 10 9 8 -log10 Err 7 6 5 ode45 ode56tsit1999 ode56vern1994 ode56elm2003 ode46b 4 3 2 3 3.2 3.4 3.6 3.8 4 4.2 log10 Nfevals Рисунок 3.2 — Зависимость глобальной погрешности от количества вычислений для модели 2 На рис. 3.3 и 3.4 приведён график изменения длины шага 𝑇𝑠𝑡𝑒𝑝 во время интегрирования для конкретных значений глобальной погрешности. Поскольку её нельзя выбрать заранее, для каждого метода были подобраны такие значения 𝑡𝑜𝑙, при которых глобальная погрешность примерно равна 10−9 для модели 1 и 10−5 для модели 2. В табл. 5 и 6 приведены значения глобальных погрешностей 𝐸𝑟𝑟, достига емых при конкретных количествах шагов интегрирования. Помимо некоторого выигрыша в точности метод 𝑜𝑑𝑒46𝑏 требует и меньших вычислительных затрат на каждом шаге по сравнению с методами шестого порядка. В сравнении с 𝑜𝑑𝑒45 метод требует таких же вычислительных затрат, однако имеет больший порядок точности и, как следствие, существенно меньшую глобальную погрешность.

26 0.35 ode45 ode56tsit1999 ode56vern1994 ode56elm2003 ode46b 0.3 Tstep 0.25 0.2 0.15 0.1 0.05 0 0.5 1 1.5 2 2.5 3 T Рисунок 3.3 — Изменение длины шага для модели 1 Таблица 5 — Точность при фиксированном количестве шагов модели 1 Кол-во −𝑙𝑜𝑔10 (𝐸𝑟𝑟) шагов 𝑜𝑑𝑒45 𝑜𝑑𝑒56𝑡𝑠𝑖𝑡1999 𝑜𝑑𝑒56𝑣𝑒𝑟𝑛1994 𝑜𝑑𝑒56𝑒𝑙𝑚2003 20 7,2431 9,3053 9,4704 9,1986 30 8,1387 10,4688 10,5491 10,2473 40 8,7382 11,3921 11,3264 11,0563 для Таблица 6 — Точность при фиксированном количестве шагов модели 2 Кол-во −𝑙𝑜𝑔10 (𝐸𝑟𝑟) шагов 𝑜𝑑𝑒45 𝑜𝑑𝑒56𝑡𝑠𝑖𝑡1999 𝑜𝑑𝑒56𝑣𝑒𝑟𝑛1994 𝑜𝑑𝑒56𝑒𝑙𝑚2003 400 4,0095 5,2410 5,8257 5,1803 500 4,4443 5,8332 6,5114 5,7427 600 4,8206 6,2936 6,6131 6,2454 для 𝑜𝑑𝑒46𝑏 9,7187 10,7620 11,6226 𝑜𝑑𝑒46𝑏 6,4222 6,9794 7,4493

27 ode45 ode56tsit1999 ode56vern1994 ode56elm2003 ode46b 0.3 0.25 Tstep 0.2 0.15 0.1 0.05 0 0 2 4 6 8 10 12 14 T Рисунок 3.4 — Изменение длины шага для модели 2 16

28 Выводы Для того, чтобы метод имел 𝑝-й порядок точности, график зависимости логарифма погрешности от логарифма длины шага должен вести себя примерно так же, как прямая линия с углом наклона arctg 𝑝. Видно, что на рис. 3.1 и 3.2 наклоны графиков для схем шестого порядка примерно равны. Это значит, что построенная схема действительно имеет шестой порядок. Выигрыш, достигаемый функцией 𝑜𝑑𝑒46𝑏 (табл. 5, табл. 6) при равном количестве шагов, объясняется тем, что свободные параметры α1 , . . . , α7 были выбраны из соображений минимизации членов седьмого порядка в разложении методической погрешности на шаге. Это приводит к тому, что на большинстве шагов главный член локальной погрешности по модулю меньше, чем у мето дов-конкурентов. Алгоритмическое же преимущество функции 𝑜𝑑𝑒46𝑏 достигается тем, что на каждом шаге совершается лишь шесть вычислений правой части СОДУ, тогда как у методов-конкурентов того же порядка — семь или восемь. На рис. 3.4 отчётливо видны «быстрые» и «медленные» области (более «медленные» участки соответствуют более «прямым» участкам траектории кос мического аппарата). Также видно, что метод 𝑜𝑑𝑒46𝑏 позволяет совершать более длинные шаги, ускоряя процесс интегрирования.

29 Заключение 1. В работе в явном виде представлено семипараметрическое семейство шестиэтапных методов шестого порядка класса B. 2. При фиксированных значениях параметров представлена расчётная схема численного интегрированияа 𝑅𝐾𝐵6(4){7𝐹 }. Найденная схемы экономичны в плане меньшего количества вычислений правой части СОДУ по сравнению с уже существующими методами. 3. На базе расчётной схемы создана функция, аналогичная встроенному интегратору MATLAB 𝑜𝑑𝑒45 в части автоматического выбора шага. При том же количестве этапов представленная функция 𝑜𝑑𝑒46𝑏 имеет на один порядок точности больше. 4. Проведено численное исследование различных моделей механики с помощью полученных методов численного интегрирования. Продемон стрирована эффективность в сравнении с известными классическими методами, в том числе с 𝑜𝑑𝑒45.

30 Список литературы 1. Runge, C. 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. Олемской, И. В. Численный метод интегрирования систем обыкновенных дифференциальных уравнений / И. В. Олемской // Математические ме тоды анализа управляемых процессов. — 1986. — С. 157—160. 8. Hofer, E. A partially implicit method for large stiff systems of ODEs with only few equations introducing small time-constants / E. Hofer // SIAM J. Numer. Anal. –– 1976. –– Vol. 13, no. 5. –– P. 645––663. 9. Sandu, A. Multirate generalized additive Runge-Kutta methods / A. Sandu, M. Günther // Numer. Math. –– 2016. –– Vol. 133, no. 3. –– P. 497––524. 10. McLachlan, R. High order multisymplectic Runge-Kutta methods / R. McLach lan, B. Ryland, Y. Sun // SIAM J. Sci. Comput. –– 2014. –– Vol. 36, no. 5. –– A2199––A2226.

31 11. Wang, D. Parametric symplectic partitioned Runge-Kutta methods with ener gy-preserving properties for Hamiltonian systems / D. Wang, A. Xiao, X. Li // Comput. Phys. Comm. –– 2013. –– Vol. 184, no. 2. –– P. 303––310. 12. Ketcheson, D. Spatially partitioned embedded Runge-Kutta methods / D. Ketcheson, C. MacDonald, S. Ruuth // SIAM J. Numer. Anal. –– 2013. –– Vol. 51, no. 5. –– P. 2887––2910. 13. Eremin, A. S. Continuous Extensions for Structural Runge—Kutta Methods / A. S. Eremin, N. A. Kovrizhnykh // Computational Science and Its Applica tions – ICCSA 2017. –– Cham : Springer International Publishing, 2017. –– P. 363––378. –– (Lecture Notes in Computer Science ; 10405). 14. Шиманчук, Д. В. Построение траектории возвращения в окрестность кол линеарной точки либрации системы Солнце–Земля / Д. В. Шиманчук, А. С. Шмыров // Вестн. С.-Петерб. ун-та, Сер. 10. — 2013. — Т. 2. — С. 76—85. 15. Clohessy, W. H. Terminal Guidance System for Satellite Rendezvous / W. H. Clohessy, R. S. Wiltshire // J. Astronaut. Sci. –– 1960. –– Vol. 27, no. 9. –– P. 653––678. 16. Tschauner, J. Rendezvous Zu Einem In Elliptischer Bahn Umlaufenden Ziel / J. Tschauner, P. Hempel // Astronautica Acta. — 1965. — Т. 11, № 2. — С. 104—109. 17. 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. 18. Хайрер, Э. Решение обыкновенных дифференциальных уравнений. Нежёсткие задачи / Пер. с англ. И. А. Кульчицкой, С. С. Филиппова под ред. С. С. Филиппова. / Э. Хайрер, С. Нёрсетт, Г. Ваннер. — М.: Мир, 1990. — 512 с. 19. Олемской, И. В. Модификация алгоритма выделения структурных особен ностей / И. В. Олемской // Вестн. С.-Петерб. ун-та, Сер. 10. — 2006. — Т. 2. — С. 46—54. 20. Олемской, И. В. Структурный подход в задаче конструирования явных одношаговых методов / И. В. Олемской // Журнал вычислительной ма тематики и математической физики. — 2003. — Т. 43, № 7. — С. 961—974.

32 21. Еремин, А. С. Разработка явного одношагового вложенного метода для систем структурно разделенных обыкновенных дифференциальных уравнений : дис. ... кандидата физико-математических наук: 01.01.07 / А. С. Еремин. — С.-Петерб. гос. ун-т, Санкт-Петербург, 2009. — 91 с. 22. 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. 23. Олемской, И. В. Методы интегрирования систем структурно разделен ных дифференциальных уравнений / И. В. Олемской. — СПб.: Изд-во С.-Петерб. ун-та, 2009. — 180 с. 24. Tsitouras, C. 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. 25. Verner, J. Strategies for Deriving New Explicit Runge—Kutta Pairs / J. Verner // Annals of Numerical Mathematics. –– 1994. –– Vol. 1. –– P. 225––244. 26. El-Mikkawy, M. E. A. A General Four-Parameter Non-FSAL Embed ded Runge–Kutta Algorithm of Orders 6 and 4 in Seven Stages / M. E. A. El-Mikkawy, M. M. M. Eisa // Applied Mathematics and Com putation. –– 2003. –– Vol. 143, no. 2. –– P. 259––267. 27. Choose an ODE Solver [Электронный ресурс]. — URL: https : / / www . mathworks.com/help/matlab/math/choose- an- ode- solver.html (дата обр. 19.01.2018). 28. Сравнительное исследование преимуществ структурных методов числен ного решения обыкновенных дифференциальных уравнений / В. Бубнов [и др.] // Труды СПИИРАН. — 2017. — Т. 53, № 4. — С. 51—72.

33 Список рисунков 1.1 Структуры допустимых зависимостей в СОДУ . . . . . . . . . . . . . 3.1 Зависимость глобальной погрешности для модели 1 . . . . . . . . . . . . . . Зависимость глобальной погрешности для модели 2 . . . . . . . . . . . . . . Изменение длины шага для модели 1 Изменение длины шага для модели 2 3.2 3.3 3.4 от количества вычислений . . . . . . . . . . . . . . . . от количества вычислений . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 . . 24 . . . . . . 25 26 27

34 Список таблиц 1 Коэффициенты явного метода класса C при 𝑚 = 𝑚0 = 𝑚1 + 1 . . . . 15 2 3 4 Размеры системы условий порядка . . . . . . . . . . . . . . . . . . . . Примеры построения условий порядка . . . . . . . . . . . . . . . . . Коэффициенты метода 𝑅𝐾𝐵6(4){7𝐹 } . . . . . . . . . . . . . . . . . 16 17 21 5 6 Точность при фиксированном количестве шагов для модели 1 . . . . Точность при фиксированном количестве шагов для модели 2 . . . . 26 26

35 Приложение А Семейство методов шестого порядка класса B Параметры метода 𝑐12 = α1 , 𝑐13 = α2 , 𝑐14 = α3 , 𝑐23 = α4 , 𝑐24 = α5 , 𝑎1132 = α6 и 𝑎2243 = α7 могут быть выбраны в качестве свободных параметров, подчиненных лишь некоторым очевидным ограничениям. В силу громоздкого вида общего решения ограничимся последовательным выражением параметров метода друг через друга. 𝑐16 = 𝑐26 = 1 𝑎1166 = 𝑎2166 = 𝑎2266 = 0 3𝑐13 − 2𝑐23 𝑐22 = 15𝑐13 𝑐23 − 10𝑐23 + 2 Введём вспомогательные переменные ξµ ≡ ξ2 = + 6 ∑︁ 6 ∑︁ (µ = 2,3,4). + 𝑐24 𝑐225 − 𝑐24 𝑐325 − 𝑐23 𝑐24 𝑐25 + 𝑐23 𝑐24 𝑐225 − 𝑐23 𝑐325 − 𝑐325 ) 𝑏25 𝑐22 (𝑐24 − 𝑐22 )(𝑐23 − 𝑐22 ) 5 + 15𝑐24 𝑐22 − 8𝑐24 − 8𝑐22 − 120(𝑐23 − 𝑐22 )𝑐23 (𝑐24 − 𝑐23 ) + 𝑐24 𝑐225 − 𝑐24 𝑐325 + 𝑐425 − 𝑐325 − 𝑐24 𝑐25 𝑐22 + 𝑐225 𝑐22 − 𝑐325 𝑐22 𝑏25 − (𝑐23 − 𝑐22 )𝑐23 (𝑐24 − 𝑐23 ) 𝑏1ν 𝑐1ν 𝑎12ν4 = ν=5 2 𝑐25 𝑐23 𝑐22 ν=µ+1 𝑏1ν 𝑐1ν 𝑎12νµ , 5 + 15𝑐24 𝑐23 − 8𝑐23 − 8𝑐24 + 120𝑐22 (𝑐24 − 𝑐22 )(𝑐23 − 𝑐22 ) 𝑏1ν 𝑐1ν 𝑎12ν3 = − ν=4 𝑐24 𝑐225 𝑐22 ξ4 = + 𝑏1ν 𝑐1ν 𝑎12ν2 = ν=3 4 (𝑐25 + 𝑐23 𝑐225 ξ3 = − 6 ∑︁ ∑︀6 5 + 15𝑐23 𝑐22 − 8𝑐23 − 8𝑐22 + 120(𝑐24 − 𝑐23 )(𝑐24 − 𝑐22 )𝑐24 + 𝑐23 𝑐225 − 𝑐23 𝑐325 + 𝑐425 − 𝑐325 − 𝑐25 𝑐23 𝑐22 + 𝑐225 𝑐22 − 𝑐325 𝑐22 𝑏25 (𝑐24 − 𝑐23 )(𝑐24 − 𝑐22 )𝑐24 Для 𝑞 = 1,2: 5𝑐𝑞4 𝑐𝑞3 − 3𝑐𝑞4 − 3𝑐𝑞3 + 2 10𝑐𝑞4 𝑐𝑞3 − 5𝑐𝑞4 − 5𝑐𝑞3 + 3 50𝑐2𝑞4 𝑐2𝑞3 − 40𝑐2𝑞4 𝑐𝑞3 + 5𝑐2𝑞4 − 40𝑐𝑞4 𝑐2𝑞3 + 35𝑐𝑞4 𝑐𝑞3 − 5𝑐𝑞4 + 5𝑐2𝑞3 − 5𝑐𝑞3 + 1 𝑏𝑞1 = 60𝑐𝑞4 𝑐𝑞3 (5𝑐𝑞4 𝑐𝑞3 − 3𝑐𝑞4 − 3𝑐𝑞3 + 2) 𝑐𝑞5 =

36 5𝑐2𝑞4 − 5𝑐𝑞4 + 1 𝑏𝑞3 = 60𝑐𝑞3 (1 − 𝑐𝑞3 )(𝑐𝑞4 − 𝑐𝑞3 )(3𝑐𝑞4 + 6𝑐𝑞3 − 2 − 10𝑐𝑞4 𝑐𝑞3 + 10𝑐𝑞4 𝑐2𝑞3 − 5𝑐2𝑞3 ) 5𝑐2𝑞3 − 5𝑐𝑞3 + 1 𝑏𝑞4 = 60𝑐𝑞4 (𝑐𝑞4 − 1)(𝑐𝑞4 − 𝑐𝑞3 )(6𝑐𝑞4 + 3𝑐𝑞3 − 2 − 10𝑐𝑞4 𝑐𝑞3 + 10𝑐2𝑞4 𝑐𝑞3 − 5𝑐2𝑞4 ) 1 × 𝑏𝑞5 = 60(5𝑐𝑞4 𝑐𝑞3 − 3𝑐𝑞4 − 3𝑐𝑞3 + 2)(3𝑐𝑞4 + 6𝑐𝑞3 − 2 − 10𝑐𝑞4 𝑐𝑞3 + 10𝑐𝑞4 𝑐2𝑞3 − 5𝑐2𝑞3 ) (10𝑐𝑞4 𝑐𝑞3 − 5𝑐𝑞4 − 5𝑐𝑞3 + 3)5 × (5𝑐𝑞4 𝑐𝑞3 − 2𝑐𝑞4 − 2𝑐𝑞3 + 1)(6𝑐𝑞4 + 3𝑐𝑞3 − 2 − 10𝑐𝑞4 𝑐𝑞3 + 10𝑐2𝑞4 𝑐𝑞3 − 5𝑐2𝑞4 ) 50𝑐2𝑞4 𝑐2𝑞3 − 60𝑐2𝑞4 𝑐𝑞3 + 15𝑐2𝑞4 + 75𝑐𝑞4 𝑐𝑞3 − 20𝑐𝑞4 − 60𝑐𝑞4 𝑐2𝑞3 + 15𝑐2𝑞3 − 20𝑐𝑞3 + 6 𝑏𝑞6 = 60(5𝑐𝑞4 𝑐𝑞3 − 2𝑐𝑞4 − 2𝑐𝑞3 + 1)(1 − 𝑐𝑞4 )(1 − 𝑐𝑞3 ) {︁(︀ 𝑐14 {︁(︀ }︁ 𝑎1243 = 20𝑐13 − )︀ 2 2 6𝑐23 (5𝑐13 − 5𝑐13 + 1) 15𝑐13 − 10 𝑐23 + 4𝑐23 − 3𝑐13 (𝑐13 − 𝑐14 ) }︁ )︀ 2 2 2 − 30𝑐13 + 150𝑐14 𝑐13 − 100𝑐14 𝑐13 − 6 + 20𝑐14 𝑐23 − 45𝑐14 𝑐13 + 5𝑐13 + 20𝑐14 𝑐13 − 4𝑐14 𝑐23 (3𝑐13 − 2𝑐23 )(15𝑐13 𝑐23 − 10𝑐23 + 2) 2𝑐12 (45𝑐213 − 30𝑐13 + 4) 𝑐24 (15𝑐13 𝑐24 𝑐23 + 2𝑐23 − 10𝑐24 𝑐23 − 3𝑐13 + 2𝑐24 )(𝑐23 − 𝑐24 ) = 6𝑐14 (𝑐13 − 𝑐14 )(5𝑐223 − 5𝑐23 + 1) 15𝑐213 𝑐23 − 10𝑐13 𝑐23 + 2𝑐23 − 𝑐13 = 45𝑐213 − 30𝑐13 + 4 {︁(︀ = 600𝑐323 − 270𝑐213 − 2025𝑐23 𝑐313 − 6750𝑐323 𝑐313 − 3900𝑐323 𝑐13 𝑎2132 = 𝑎2144 𝑎2233 𝑎2244 − 520𝑐223 − 9000𝑐213 𝑐223 − 1200𝑐13 𝑐23 + 9000𝑐323 𝑐213 + 6750𝑐223 𝑐313 + 120𝑐13 )︀ + 160𝑐23 + 3600𝑐223 𝑐13 + 3150𝑐213 𝑐23 − 16 𝑐324 + (︀ + 630𝑐213 𝑐23 − 180𝑐213 − 120𝑐13 𝑐23 − 900𝑐213 𝑐223 + 2025𝑐223 𝑐313 − 120𝑐223 𝑐13 + + 405𝑐313 − 1350𝑐23 𝑐313 − 320𝑐323 − 1350𝑐323 𝑐213 + 36𝑐13 − )︀ − 24𝑐23 + 160𝑐223 + 1200𝑐323 𝑐13 𝑐224 + (︀ + 90𝑐213 𝑐23 − 135𝑐23 𝑐313 − 60𝑐323 𝑐13 + 24𝑐13 𝑐23 + 180𝑐213 𝑐223 − 120𝑐223 𝑐13 + 40𝑐323 )︀ (︀ − 18𝑐213 − 8𝑐223 𝑐24 + 5400𝑐323 𝑐13 − 720𝑐323 + 96𝑐223 + 4050𝑐223 𝑐313 − 1620𝑐213 𝑐223 + 1680𝑐423 − 27000𝑐523 𝑐213 − 1200𝑐523 − 360𝑐223 𝑐13 − 72𝑐13 𝑐23 − 810𝑐23 𝑐313 + 540𝑐213 𝑐23 − 20250𝑐423 𝑐313 + 20250𝑐523 𝑐313 − 8100𝑐323 𝑐213 − 14400𝑐423 𝑐13 }︁ )︀ 4 2 5 + 32400𝑐23 𝑐13 + 10800𝑐23 𝑐13 𝑎2243 × −1 6𝑐24 ((15𝑐13 𝑐23 + 2 − 10𝑐23 )𝑐24 + 2𝑐23 − 3𝑐13 )(5𝑐223 + 5𝑐23 + 1)(45𝑐213 − 30𝑐13 + 4) 1 𝑎1121 = 𝑎1122 = 𝑐12 2 ×

37 1 𝑎1155 = (1 − 𝑐15 ) 2 1 𝑏15 𝑎1165 = (1 − 𝑐15 ) 2 𝑏16 𝑐13 − 𝑐12 1 𝑎1131 = 𝑐13 − 𝑎1132 2 𝑐13 1 𝑐12 𝑎1133 = 𝑐13 − 𝑎1132 2 𝑐13 𝑏13 (𝑐15 − 𝑐13 )(1 − 𝑐13 ) · 𝑎1132 𝑎1142 = 𝑏14 (𝑐15 − 𝑐14 )(𝑐14 − 1) 𝑏13 (𝑐14 − 𝑐13 )(1 − 𝑐13 ) · 𝑎1132 𝑎1152 = − 𝑏15 (𝑐15 − 𝑐14 )(𝑐15 − 1) 𝑏13 (𝑐14 − 𝑐13 )(𝑐13 − 𝑐15 ) 𝑎1162 = 𝑎1132 · 𝑏16 (1 − 𝑐14 )(1 − 𝑐15 ) 30𝑐13 𝑐214 + 30𝑐314 𝑐12 − 30𝑐314 𝑐13 − 30𝑐214 𝑐12 − 6𝑐14 𝑐13 + 6𝑐12 𝑐14 𝑎1132 + 𝑎1141 = 6(5𝑐213 − 5𝑐13 + 1)𝑐213 𝑐313 − 𝑐13 𝑐214 − 20𝑐313 𝑐14 + 5𝑐214 𝑐213 + 15𝑐413 𝑐14 + 3𝑐14 𝑐213 + 6(5𝑐213 − 5𝑐13 + 1)𝑐213 2𝑐14 − 10𝑐14 𝑐13 − 𝑐13 + 15𝑐213 𝑐14 𝑎1144 = 6(5𝑐213 − 5𝑐13 + 1) 𝑏14 (𝑐14 − 1)(2𝑎1144 + 𝑐14 − 1) 𝑎1154 = 2𝑏15 (1 − 𝑐15 ) (︀ )︀ 𝑏14 2𝑎1144 (𝑐14 − 𝑐15 ) + (2𝑐15 − 𝑐14 − 1)(1 − 𝑐14 ) 𝑎1164 = 2𝑏16 (𝑐15 − 1) 𝑐215 (𝑎1152 𝑐12 + 𝑎1154 𝑐14 + 𝑎1155 𝑐15 ) 𝑎1153 = − 2𝑐13 𝑐13 1 𝑎1162 𝑐12 + 𝑎1164 𝑐14 + 𝑎1165 𝑐15 𝑎1163 = − 2𝑐13 𝑐13 𝑤 ∑︁ 𝑎11𝑤1 = 𝑐1𝑤 − 𝑎11𝑤ν , 𝑤 = 5,6 ν=2 𝑏25 (1 − 𝑐25 ) 𝑏16 𝑐13 (2𝑐22 − 𝑐13 ) = 2𝑐22 2 𝑐 = 13 2𝑐22 ξ4 − 𝑏24 (1 − 𝑐24 ) = 𝑏15 (𝑐15 − 1) 𝑎1265 = 𝑎1231 𝑎1232 𝑎1254

38 𝑏24 𝑐15 (1 − 𝑐24 ) − ξ4 𝑏15 (𝑐15 − 1) (1 − 𝑐14 )𝑎1243 𝑏14 + ξ3 + (𝑐23 − 1)𝑏23 𝑎1253 = 𝑏15 (𝑐15 − 1) (𝑐15 − 𝑐14 )𝑎1243 𝑏14 + 𝑐15 (𝑐23 − 1)𝑏23 + ξ3 𝑎1263 = 𝑏16 (1 − 𝑐15 ) (𝑐14 − 1)𝑏14 𝑐23 𝑎1243 ξ2 (1 − 𝑐13 )𝑐213 𝑏13 + (1 − 𝑐14 )𝑐214 𝑏14 𝑎1252 = + + 𝑏15 𝑐22 (𝑐15 − 1) 𝑏15 (𝑐15 − 1) 2𝑏15 𝑐22 (𝑐15 − 1) ξ2 (𝑐15 − 𝑐13 )𝑐213 𝑏13 + (𝑐15 − 𝑐14 )𝑐214 𝑏14 (𝑐15 − 𝑐14 )𝑏14 𝑐23 𝑎1243 − − 𝑎1262 = 𝑏16 𝑐22 (𝑐15 − 1) 𝑏16 (𝑐15 − 1) 2𝑏16 𝑐22 (𝑐15 − 1) 2 (𝑐23 − 𝑐22 )𝑎1243 𝑐 𝑎1241 = − 14 + 𝑐14 𝑐22 2𝑐22 2 𝑐 𝑐23 𝑎1243 + 14 𝑎1242 = − 𝑐22 2𝑐22 𝑤−1 ∑︁ 𝑎12𝑤1 = 𝑐1𝑤 − 𝑎12𝑤ν , 𝑤 = 5,6 𝑎1264 = ν=2 𝑐22 (2𝑐12 − 𝑐22 ) 𝑐12 2 𝑐 = 22 𝑐12 (𝑐12 − 𝑐13 ) 𝑐23 (2𝑐13 − 𝑐23 ) = 𝑎2132 + 𝑐13 2𝑐13 2 𝑐12 𝑐 𝑎2132 = 23 − 2𝑐13 𝑐13 𝑏15 (1 − 𝑐15 )2 = 2𝑏25 (1 − 𝑐25 ) 𝑏15 (1 − 𝑐15 )(1 + 𝑐15 − 2𝑐25 ) = 2𝑏26 (1 − 𝑐25 ) 𝑏23 𝑎2132 (𝑐23 − 1)(𝑐25 − 𝑐23 ) =− 𝑏24 (𝑐24 − 1)(𝑐25 − 𝑐24 ) 𝑏23 𝑎2132 (𝑐23 − 1)(𝑐24 − 𝑐23 ) = 𝑏25 (𝑐25 − 1)(𝑐25 − 𝑐24 ) 𝑏23 𝑎2132 (𝑐25 − 𝑐23 )(𝑐24 − 𝑐23 ) =− 𝑏26 (1 − 𝑐25 )(1 − 𝑐24 ) {︁ 1 (1 − 𝑐14 )2 }︁ = (1 − 𝑐24 )𝑎2144 𝑏24 − 𝑏14 𝑏25 (𝑐25 − 1) 2 {︁ 1 (1 − 𝑐14 )(1 + 𝑐14 − 2𝑐25 ) }︁ = (𝑐25 − 𝑐24 )𝑎2144 𝑏24 + 𝑏14 𝑏26 (1 − 𝑐25 ) 2 𝑎2121 = 𝑎2122 𝑎2131 𝑎2133 𝑎2155 𝑎2165 𝑎2142 𝑎2152 𝑎2162 𝑎2154 𝑎2164

39 𝑐12 𝑐14 𝑐224 − 𝑎2142 − 𝑎2144 𝑎2143 = 2𝑐13 𝑐13 𝑐13 5 ∑︁ }︀ 1 {︀ 𝑐22𝑤 − 𝑎21𝑤ν 𝑐1ν , 𝑎21𝑤3 = 𝑐13 2 𝑤 = 5,6 ν=2ν̸=3 𝑎21𝑤1 = 𝑐2𝑤 − 𝑤 ∑︁ 𝑎21𝑤ν , 𝑤 = 4,5,6 ν=2 𝑐22 𝑎2221 = 𝑎2222 = 2 𝑏23 (𝑐23 − 1)(2𝑎2233 + 𝑐23 − 1) 𝑏24 (𝑐24 − 1)𝑎2243 + 𝑎2253 = 2𝑏25 (1 − 𝑐25 ) 𝑏25 (1 − 𝑐25 ) 𝑏23 (𝑐25 − 𝑐23 ) 𝑏23 (1 − 𝑐23 )(1 + 𝑐23 − 2𝑐25 ) 𝑏24 (𝑐25 − 𝑐24 ) 𝑎2243 + 𝑎2233 + 𝑎2263 = 𝑏26 (1 − 𝑐25 ) 𝑏26 (1 − 𝑐25 ) 2𝑏26 (1 − 𝑐25 ) {︁ }︁ 𝑏24 (𝑐24 − 1) 𝑎2254 = 2𝑎2244 + 𝑐24 − 1 2𝑏25 (1 − 𝑐25 ) {︁ }︁ 𝑏24 𝑎2264 = 2(𝑐25 − 𝑐24 )𝑎2244 + (1 − 𝑐24 )(1 + 𝑐24 − 2𝑐25 ) 2𝑏25 (1 − 𝑐25 ) 1 𝑎2255 = (1 − 𝑐25 ) 2 𝑏25 (1 − 𝑐25 ) 𝑎2265 = 2𝑏26 µ 1 𝑐22µ ∑︁ 𝑎22µ1 = 𝑐2µ − { + (𝑐22 − 𝑐2ν )𝑎22µν }, µ = 3,4,5,6 𝑐22 2 ν=3 µ 𝑎22µ2 1 𝑐22µ ∑︁ { − 𝑐2ν 𝑎22µν }, = 𝑐22 2 ν=3 µ = 3,4,5,6

40 Приложение Б Текст программы ode46b Листинг Б.1 Функция численного интегрирования систем класса B function varargout = ode46b ( ode , tspan , y0 , options , varargin ) solver_name = ’ ode45 ’; 5 % Stats nsteps = 0; nfailed = 0; nfevals = 0; 10 % Output FcnHandlesUsed = isa ( ode , ’ function_handle ’) ; output_sol = ( FcnHandlesUsed && ( nargout ==1) ) ; sol = []; f3d = []; 15 if output_sol sol . solver = solver_name ; sol . extdata . odefun = ode ; sol . extdata . options = options ; sol . extdata . varargin = varargin ; 20 end % Handle solver arguments for i =1:4 [ neq , tspan , ntspan , next , t0 , tfinal , tdir , y0 , f0 , odeArgs , odeFcn , ... 25 options , threshold , rtol , normcontrol , normy , hmax , htry , htspan , dataType ] = ... odearguments ( FcnHandlesUsed , solver_name , ode , tspan , y0 , options , varargin ) ; end nfevals = nfevals + 1; 30 % Handle the output if nargout > 0 outputFcn = odeget ( options , ’ OutputFcn ’ ,[] , ’ fast ’) ; else

41 outputFcn = odeget ( options , ’ OutputFcn ’ , @odeplot , ’ fast ’) ; 35 end outputArgs = {}; if isempty ( outputFcn ) haveOutputFcn = false ; else 40 haveOutputFcn = true ; outputs = odeget ( options1 , ’ OutputSel ’ ,1: neq , ’ fast ’) ; if isa ( outputFcn , ’ function_handle ’) % With MATLAB 6 syntax pass additional input arguments to outputFcn . outputArgs = varargin ; 45 end end refine = max (1 , odeget ( options , ’ Refine ’ ,4 , ’ fast ’) ) ; % Handle the event function 50 [ haveEventFcn , eventFcn , eventArgs , valt , teout , yeout , ieout ] = ... odeevents ( FcnHandlesUsed , odeFcn , t0 , y0 , options , varargin ) ; t = t0 ; y = y0 ; 55 % Allocate memory if we ’ re generating output . nout = 0; tout = []; yout = []; if nargout > 0 60 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 ) ; 65 nout = 1; tout ( nout ) = t ; yout (: , nout ) = y ; end 70 % Initialize method parameters . pow = 1/5; [ B11 B12 B21 B22 A1 A2 E1 E2 ] = coeff ; f ( neq ,7) = zeros (1 ,1 , dataType ) ; 75

42 hmin = 16* eps ( t ) ; % Compute an initial step size h using y ’( t ) . absh = min ( hmax , htspan ) ; 80 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 ny1 = structure ( ’1 ’) ; ny2 = structure ( ’2 ’) + ny1 ; % Initialize the output function . if haveOutputFcn 95 feval ( outputFcn ,[ t tfinal ] , y ( outputs ) , ’ init ’ , outputArgs {:}) ; end % THE MAIN LOOP 100 done = false ; while ~ done 105 110 % 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. 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 1.1* absh >= abs ( tfinal - t ) h = tfinal - t ; absh = abs ( h ) ; done = true ; end 115 % LOOP FOR ADVANCING ONE STEP .

43 120 125 130 nofailed = true ; % no failed attempts while true hA1 = h * A1 ; hA2 = h * A2 ; hB11 = h * B11 ; hB12 = h * B12 ; hB21 = h * B21 ; hB22 = h * B22 ; for i =2:6 for j =1: ny1 yt = [ y( 1: ny1 ) + f ( 1: ny1 ,:) * hB11 (: ,i -1) y ( ny1 +1: ny2 ) + f ( ny1 +1: ny2 ,:) * hB12 (: ,i -1) ]; f (j , i ) = feval ( odeFcn , t + hA1 (i -1) , yt , j , odeArgs {:}) ; end for j = ny1 +1: ny2 yt = [ y( 1: ny1 ) + f ( 1: ny1 ,:) * hB21 (: ,i -1) y ( ny1 +1: ny2 ) + f ( ny1 +1: ny2 ,:) * hB22 (: ,i -1) ]; f (j , i ) = feval ( odeFcn , t + hA2 (i -1) , yt , j , odeArgs {:}) ; end 135 140 end 145 150 tnew = t + hA1 (6) ; if done tnew = tfinal ; end % Hit end point exactly . ynew = y + [ f( 1: ny1 ,:) * hB11 (: ,6) f ( ny1 +1: ny2 ,:) * hB22 (: ,6) ]; f (: ,7) = feval ( odeFcn , tnew , ynew , odeArgs {:}) ; nfevals = nfevals + 6; 155 % Estimate the error . NNrejectStep = false ;

44 160 165 170 175 180 185 190 err = absh * norm (([ f (1: ny1 ,:) * E1 ; f ( ny1 +1: ny2 ,:) * E2 ]) ... ./ max ( max ( abs ( y ) , abs ( ynew ) ) , threshold ) , inf ) ; % Accept the solution only if the weighted error is no more than the % tolerance rtol . Estimate an h that will yield an error of rtol on % the next step or the next try at taking this step , as the case may be , % and use 0.8 of this value to avoid failures . if err > rtol % Failed step nfailed = nfailed + 1; if absh <= hmin warning ( message ( ’ MATLAB : ode45 : IntegrationTolNotMet ’ , sprintf ( ’% e ’ , t ) , sprintf ( ’% e ’ , hmin ) ) ) ; solver_output = odefinalize ( solver_name , sol ,... outputFcn , outputArgs ,... 0 , [ nsteps , nfailed , nfevals ] ,... nout , tout , yout ,... haveEventFcn , teout , yeout , ieout ,... { f3 , idxNonNegative }) ; if nargout > 0 varargout = solver_output ; end return ; end if nofailed nofailed = false ; if NNrejectStep absh = max ( hmin , 0.5* absh ) ; else absh = max ( hmin , absh * max (0.1 , 0.8*( rtol / err ) ^ pow ) ) ; end else absh = max ( hmin , 0.5 * absh ) ; end h = tdir * absh ;

45 done = false ; 195 else % Successful step NNreset_f7 = false ; 200 break ; end end nsteps = nsteps + 1; 205 210 215 220 225 230 if output_sol nout = nout + 1; if nout > length ( tout ) tout = [ tout , zeros (1 , chunk , dataType ) ]; % requires chunk >= refine yout = [ yout , zeros ( neq , chunk , dataType ) ]; f3d = cat (3 , f3d , zeros ( neq ,7 , chunk , dataType ) ) ; end tout ( nout ) = tnew ; yout (: , nout ) = ynew ; f3d (: ,: , nout ) = f ; end if done break end % If there were no failures compute a new h . if nofailed % Note that absh may shrink by 0.8 , and that err may be 0. temp = 1.25*( err / rtol ) ^ pow ; if temp > 0.2 absh = absh / temp ; else absh = 5.0* absh ; end end % Advance the integration one step . t = tnew ;

46 235 y = ynew ; 240 if NNreset_f7 % Used f7 for unperturbed solution to interpolate . % Now reset f7 to move along constraint . f (: ,7) = feval ( odeFcn , tnew , ynew , odeArgs {:}) ; nfevals = nfevals + 1; end f (: ,1) = f (: ,7) ; % Already have f ( tnew , ynew ) 245 end solver_output = odefinalize ( solver_name , sol ,... outputFcn , outputArgs ,... 0 , [ nsteps , nfailed , nfevals ] ,... 250 nout , tout , yout ,... haveEventFcn , teout , yeout , ieout ,... { f3d , false }) ; 255 if nargout > 0 varargout = solver_output ; end end 260 function [ B11 , B12 , B21 , B22 , A1 , A2 , E1 , E2 ] = coeff B11 = [ [ 1/9 , 1/9 , 0, 0, 0, 0 , 0] [ 1/12 , 0, 1/12 , 0, 0, 0 , 0] 265 [ -1/44 , 0, 9/22 , 5/44 , 0, 0 , 0] [ 7/36 , 0, 0, 5/9 , 1/12 , 0 , 0] [ -3/7 , 0, 9/8 , -5/28 , 27/56 , 0 , 0] [ 7/150 , 0 , 27/100 , 11/30 , 27/100 , 7/150 , 0] ] ’; 270 B12 = [ [ 2/9 , 0, 0, 0, 0, 0, [ 5/48 , 1/16 , 0, 0, 0, 0, [ 37/176 , 243/176 , -12/11 , 0, 0, 0, [ -635/432 , -167/16 , 100/9 , 44/27 , 0, 0, 275 [ 29/4 , 1377/28 , -1425/28 , -11/2 , 27/28 , 0, [ 7/150 , 0, 27/100 , 11/30 , 27/100 , 7/150 , ] ’; 0] 0] 0] 0] 0] 0]

47 A1 = sum ( B11 ) ; 280 B21 = [ [ 1/9 , 1/9 , 0, 0, 0, 0 , 0] [ 7/48 , 3/16 , -1/6 , 0, 0, 0 , 0] [ -31/176 , -81/176 , 45/44 , 5/44 , 0, 0 , 0] [ 73/144 , 15/16 , -5/4 , 5/9 , 1/12 , 0 , 0] 285 [ -39/28 , -81/28 , 279/56 , -5/28 , 27/56 , 0 , 0] [ 7/150 , 0 , 27/100 , 11/30 , 27/100 , 7/150 , 0] ] ’; B22 = [ [ 1/9 , 1/9 , 0, 0, 0, 0, 290 [ 7/48 , 3/16 , -1/6 , 0, 0, 0, [ -185/1584 , -123/880 , 2/3 , 89/990 , 0, 0, [ 1031/3888 , -53/144 , 65/324 , 317/486 , 1/12 , 0, [ -29/63 , 15/7 , -103/168 , -139/252 , 27/56 , 0, [ 7/150 , 0, 27/100 , 11/30 , 27/100 , 7/150 , 295 ] ’; A2 = sum ( B22 ) ; E1 = [11/25 , 0 , -99/100 , 11/10 , -99/100 , -14/25 , 1] ’/24; E2 = E1 ; 300 end 0] 0] 0] 0] 0] 0]

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв