Содержание к диссертации
Введение
1. Обзор литературы 7
1.1. Одношаговыс методы 8
1.2. Жесткие уравнение и понятие устойчивости 13
1.3. Итерационные процессы 16
1.4. Автоматический контроль точности вычислений 21
2. Е-методы и их устойчивость 32
2.1. Пример построения Е-метода 32
2.2. Е-методы со старшими производными 34
2.3. А-устойчивость Е-методов со старшими производными . 40
2.4. Формулировка комбинированных Е-методов 43
2.5. Исследование сходимости комбинированных Е-методов . 50
3. Эффективная реализация Е-методов и контроль точности 61
3.1. Применение предикторов 61
3.2. Итерации Ньютона общего вида 71
3.3. Локально-глобальный контроль точности 83
А. Библиотека программ 105
АД. Функции и структуры для численного решения обыкновенных дифференциальных уравнений 105
А.1.1. Заголовочный файл numerics 105
А.1.2. Структура odestruct 106
А.1.3. Структура solveropt 107
А.1.4. Структура statistics 109
А.1.5. Заголовочный файл emeth.h 110
А.1.6. Пример использования функций библиотеки 111
Литература
- Жесткие уравнение и понятие устойчивости
- Автоматический контроль точности вычислений
- Формулировка комбинированных Е-методов
- Функции и структуры для численного решения обыкновенных дифференциальных уравнений
Введение к работе
Известно, что многие математические модели динамического типа в области биологии, медицины, механики, электротехники, химической кинетики и т. д. представляют собой системы обыкновенных дифференциальных уравнений (ОДУ). Очень часто аналитическое решение указанных систем найти крайне сложно или вообще невозможно. Поэтому для практического исследования таких моделей используют разного рода численные методы.
Естественно, что наибольший интерес представляют высокоточные методы, имеющие высокий порядок сходимости. В виду того, что построение данных методов представляет собой довольно трудоемкий процесс, актуальным является разработка методик, позволяющих облегчить процесс формулировки методов и определения порядка их сходимости. В диссертации указанную проблему предлагается решить с помощью явного вывода класса методов, использующих старшие производные.
Отметим также, что практическое использование численных методов невозможно без внесения определенных ошибок в получаемые результаты. Поэтому крайне актуальной задачей является исследование методик оценки локальной и глобальной ошибок для полученного класса одно-шаговых методов и построение на их основе алгоритмов, позволяющих достичь любую разумную наперед заданную точность вычислений в автоматическом режиме.
Таким образом, целью диссертации является разработка, обоснование и тестирование эффективных численных методов, использующих старшие производные, для решения обыкновенных дифференциальных уравнений. При этом важнейшей задачей, имеющей большое прикладное значение, является обеспечение заданной точности вычислений в автоматическом режиме.
Диссертация имеет как теоретическую, так и практическую ценность. Полученные в ней теоретические результаты представляют интерес для численного решения нелинейных дифференциальных уравнений с контролем глобальной ошибки, что в прикладном аспекте создает плодотворную базу для разработки автоматизированных систем математического моделирования в различных областях науки. При этом способность обеспечивать заданную пользователем точность вычислений в автоматическом режиме существенно повышает достоверность математического моделирования на практике. ^ Приведем краткое содержание диссертации, состоящей из трех глав и приложения.
Первая глава, носит вводный, вспомогательный характер и обозначает основные достижения в области исследования одношаговых методов, включая алгоритмы вычисления и управления их локальной и глобальной ф ошибками в процессе счета.
В параграфе 1.1 изложены базовые результаты, полученные в области одногошаговых методов для численного решения ОДУ. Приведена последовательность развития теории указанных методов от явных к методам со старшими производными. Параграф 1.2 посвящен рассмотрению жестких задач и связанных с ними аспектов теории устойчивости одношаговых методов. Необходимость применения итерационных алгоритмов для реализации неявных методов стало причиной детального исследования теории итерационных процессов в параграфе 1.3. Наконец, параграф 1.4 посвящен такому важному аспекту как построение алгоритмов для автоматического контроля точности на практике.
Вторая глава, состоящая из шести параграфов, посвящена решению задачи построения исследуемого класса методов и рассмотрению основ-
Л ных аспектов их практической реализации. Именно глава 2 содержит основные теоретические результаты данной диссертации.
В параграфах 2.1 и 2.2 приведена методика построения методов со старшими производными на основе использования интерполяционных полиномов Эрмита, выведены явные формулы для вычисления коэффициентов методов, показана сходимость всего класса методов и выведены оценки их погрешности. Параграф 2.3 посвящен исследованию устойчивости нового класса численных методов.
В параграфах 2.4 и 2.5 рассмотрены различные аспекты практической реализации указанных методов в силу их неявного характера. Во-первых, на основе метода простых итераций и итераций Ньютона были построены три класса комбинированных алгоритмов. Во-вторых, исследован вопрос о сходимости полученных методов и выведены оценки их погрешности. В-третьих, приведены результаты численных экспериментов, подтверждающих теоретические выкладки.
Третья глава, состоящая из трех параграфов, посвящена решению задачи эффективной реализации предложенных в диссертации методов, а
4 также обеспечению заданной точности вычислений на основе автоматического управления размером шага интегрирования.
Во-первых, в параграфе 3.1 рассмотрен вопрос о применении предикторов для вычисления стартовых значений итерационных процессов, что в практическом плане приводит к сокращению числа итераций, требуемых для достижения максимального порядка сходимости. Затем, в параграфе 3.2 выведена асимптотическая оценка погрешности для методов ньютоновского типа, что позволило разработать четвертый тип комбинированных алгоритмов, построенный с использованием упрощенных итераций Ньютона и направленный на сокращение вычислительных затрат для каждой отдельной итерации. В завершении, в параграфе 3.3 ранее изученные способы оценки и автоматического контроля локальной и глобальной ошибок адаптированы для Е-методов со старшими производными. Все теоретические выводы подтверждены результатами численных экспериментов.
В приложении, состоящем из одного параграфа, дано описание библиотеки программ, разработанной и использованной для проведения всех практических расчетов в рамках работы над диссертацией. Библиотека является набором процедур, написанных на языке C++ и ориентированных на решение систем ОДУ с заданной точностью.
Итак, основные результаты диссертации следующие:
Построено семейство Л-устойчивых одношаговых коллокационных методов, использующих старшие производные, для нахождения численного решения систем ОДУ, в том числе и жестких задач.
Исследована сходимость и точность указанных численных методов.
Разработана методика корректной реализации Е-методов со старшими производными на практике.
Исследован вопрос применения предикторов и выведен рекомендуемый тип предикторов с соответствующей коррекцией алгоритмов.
Выведена асимптотическая оценка погрешности для методов ньютоновского типа и предложена на ее основе схема упрощенных итераций Ньютона, оптимизирующих реализацию разработанных в диссертации алгоритмов.
Предложенные ранее способы оценки и автоматического контроля локальной и глобальной ошибок адаптированы для Е-методов со старшими производными.
Кроме того, в прикладном аспекте важной является разработка библиотеки программ, позволяющей решать системы обыкновенных дифференциальных уравнений с любой разумной наперед заданной точностью.
Таким образом, в диссертации решен ряд трудных и актуальных задач теории численных методов для математических моделей динамического типа. ^ Основные результаты диссертации доложены: на XXV-XXVI конфе- ренции молодых ученых механико-математического факультета МГУ им. М.В. Ломоносова (Москва, 2003-2004), на пятой международной научно-технической конференции "Математическое моделирование физических, экономических, технических, социальных систем и процес-
4 сов" (Ульяновск, 2003), на 2nd International Conference on Computer Science and its Applications (Assisi, Italy, 2004), на 4th International Conference on Computer Science (Krakow, Poland, 2004), на научно-исследовательских семинарах Вычислительного центра РАН, Института вычислительной математики РАН, механико-математического факультета МАИ, механико-математического факультета МГУ им. М.В. Ломоносова и механико-математического факультета УлГУ.
По теме диссертации опубликовано 10 работ, в том числе [22]-[23], [32]—[35], [116, 117]. При этом результаты из работ, подготовленных в соавторстве, либо использованы частично, в соответствии с личным вкладом каждого автора ([22, 24, 25, 23]), либо приведены в переработанном виде ([116, 117]).
Автор считает приятной обязанностью выразить искреннюю благо- дарность своему научному руководителю д.ф.-м.н., доценту Куликову Г.Ю. за постоянное внимание и помощь в работе, а также академику РАН Бахвалову Н.С., профессорам Семушину И.В. и Горбунову В.К. за плодотворное обсуждение результатов, полученных на разных этапах подготовки диссертации.
Кроме того, необходимо отметить, что развитие данной тематики в России было бы невозможно без финансовой поддержки со стороны: Российского фонда фундаментальных исследований (проект № 98-01-00006, 1998-2000; проект № 01-01-00066, 2001-н. вр.).
Жесткие уравнение и понятие устойчивости
Раздел 1.2 отражает еще одно преимущество применения НРК-методов по сравнению с ЯРК-методами. Дело в том, что существует ряд задач, для которых явные методы работают только в случае очень малых шагов интегрирования или не работают вообще. Такие задачи принято называть жесткими. В качестве примера жесткой задачи можно привести следующее одномерное уравнение (см. [69]): x (t) = -5Q(x{t) - cos t), х{0) = 2500/2501
При применении явного метода Эйлера и увеличении длины шага больше 2/50 численное решение уходит слишком далеко от точного решения, причем колебания возникают все более сильные.
В качестве других примеров можно привести пример Робертсона (см. [133]) системы, описывающей химические реакции; уравнение Ван-дер-Поля; уравнение диффузии; ряд примеров из механики о движении упругого стержня (см. [99]).
Впервые анализ явлений неустойчивости был выполнен в работе Куранта, Фридрихса и Леви [26] для гиперболических уравнений. Позднее анализ устойчивости проводился многими авторами, ранняя работа по данному вопросу — Гийу и Лаго [92]. Однако стандартным подходом является исследование свойств устойчивости на основе применения заданного численного метода к тестовому уравнению Далквиста (см., например, [71], [72]) x {t) = \x(t), t 0, х(0) = х, (1.2.1) где константа А Є С — множеству комплексных чисел. Применяя численный метод для решения (1.2.1), приходим к выражению xk+l = R{z)xh (1.2.2) где z = тХ.
Нетрудно видеть, что для ЯРК-методов R{z), которая называется функцией устойчивости, представляет собой многочлен с вещественными коэффициентами над полем комплексных чисел. Например, все ЯРК-методы с порядком р равным числу стадий s обладают функцией устойчивости вида Множество 5 = {геС;ДЙ 1} (1.2.3) называется областью абсолютной устойчивости. Если область абсолютной устойчивости содержит левую полуплоскость Re (z) 0, то численный метод называют А-устойчивым (или абсолютно устойчивым), а функцию R(z) — А-приемлемой (см., например, [8] и [44]). В силу (1.2.2) и (1.2.3), из Л-устойчивости следует невозрастание численного решения при любом размере шага интегрирования г, что также характерно для точного решения задачи (1.2.1) при Re (Л) 0.
В силу специфики функции R(z), область устойчивости ЯРК-методов не распространяется на всю левую полуплоскость, поэтому для задач, матрица Якоби которых имеет большие собственные значения А, размер шага интегрирования должен выбираться очень малым, чтобы обеспечить попадание =г А в область устойчивости. Это объясняет существование задач, для которых ЯРК-методы работают плохо.
Еще жестче проблема устойчивости становиться при использовании переменного шага интегрирования (например, при реализации методик автоматического управления длиной шага интегрирования).
Применив методику контроля, записанную в самом общем виде, к уравнению (1.2.1) приводим к следующему вычислительному процессу: хк+1 = R(rk\)xk, еггк — Е(тк\)хк, тк+{ = 7 ( ), (1.2.4) СТ"Гк где еггк — оценка погрешности, B(z) = R{z) — R{z), a = 1/(р-\- 1), е — желаемая точность метода, задаваемая пользователем, и р — порядок R. В [100, 101, 102] было показано, что приведенный выше процесс имеет неподвижную точку, которая лежит на границе области устойчивости. Матрица Якоби отображения (1.2.4) относительно хк и тк в этой неподвижной точке будет / і Re(4z) \ И- -«-( )) (1-2-5) Определение 1.7. Если спектральный радиус радиус матрицы С удовлетворяет условию Р{С) 1, то метод называют SC-устойчивым (Step-Control stable) в точке z = TkX. Таким образом, число SC-устойчивых ЯРК - методов определяется их областью устойчивости и является еще более ограниченным.
Большая работа по поиску ЯРК-методов с наиболее широкими областями SC-устойчивости приведана в работе Хайма и Холла [94]. Результаты показали, что чем больше желаемая область SC-устойчивости, тем больше становится константа погрешности метода. Таким образом, возможен только компромиссный вариант, ограничивающий область задач, для которых применение ЯРК-методов эффективно.
Тем не менее существует целый ряд практических задач для которых применение ЯРК-методов является приоритетным. В основном это задачи большой размерности, обычно не очень жесткие, и с собственными числами, про которые известно, что они лежат в определенной области. В частности, много работ посвящено построению и изучению методов с областями устойчивости, расширенными вдоль отрицательной вещественной оси, которые поэтому подходят для интегрирования по времени систем параболических дифференциальных уравнений с частными производными. Наиболее существенный вклад в это направление внесли работы Лебедева с соавт. [27, 29, 30, 120]. Тем не менее, применимость ЯРК-методов все равно ограничена с точки зрения устойчивости, что и послужило еще одной причиной для развития многошаговых и неявных одношаговых методов.
Многошаговые методы явились первыми численными алгоритмами, предложенными для решения жестких дифференциальных уравнений (см. [69]). С появлением книги Гира [90] основанные на этих методах программы приобрели наибольшую популярность и стали наиболее широко использоваться для всех жестких расчетов. Однако исследование устойчивости для многошаговых методов привело к знаменитой теореме Далквиста [72] о том, что А-устойчивые многошаговые методы не могут иметь высокий порядок сходимости. В настоящий момент работы по преодолению данного барьера ведуться либо в направлении изучения методов с несколько ослабленными требованиями устойчивости, либо D построении и изучении новых классов одношаговых методов. В любом случае первые методы гарантированно действуют только для ограниченного круга задач, а применение последних приводит к значительному увеличению вычислительных затрат. Таким образом, наиболее приемлемой альтернативой поиска высокоточных, А-устойчивых и оптимальных по вычислительным затратам методов являются исследования в области неявных одношаговых методов.
Автоматический контроль точности вычислений
Другая особенность современных методов численного интегрирования обыкновенных дифференциальных уравнений заключается в наличии алгоритмов автоматического изменения размера шага в зависимости от поведения решения задачи (1.1.1). Данная степень свободы позволяет в некотором смысле оптимизировать вычислительной процесс как по точности, так и по времени счета. Именно этот вопрос представляет собой суть настоящего раздела.
По определению, локальной ошибкой метода в точке & называется разность между точным значением решения задачи (1.1.1) x(tk) и его приближенным значением полученным с помощью того или иного численного метода Хк- Если предположить, что указанный метод имеет порядок s, то справедливо следующее разложение локальной ошибки ф) - xk = 4+i( i)rs+1 + 0(TS+\ (1.4.1) где ф8+х(ііі-і) называют главным членом локальной ошибки.
Практически все известные на сегодняшний момент численные методы с автоматическим выбором шага интегрирования основаны на вычислении и контроле именно ф3+1 (см., например, [1, 5, 8, 9, 14, 31, 40, 42, 43, 44, 105]). И несмотря на то, что такой способ построения "оптимального" по затратам машинного времени разбиения отрезка интегрирования для заданной точности решения поставленной задачи имеет целый ряд недостатков, на которых мы остановимся позднее, из-за простоты практической реализации он остается весьма популярным.
В целом существует два принципиально разных способа оценки и контроля локальной ошибки: правило Рунге и применение численных методов разных порядков (см., например, [5, б, 31, 40, 42, 43, 46, 97, 98]). Оба способа опираются на следующее предположение об асимптотическом представлении локальной ошибки метода порядка s (см., например, [1, 5, 6, 43]): Дя +і = Е Й( И+1 + 0{т+2). (1.4.2) i=s
Правило Рунге основывается на разложении локальной погрешности (1.4.2) при S = s и использует два значения: Xk+i - численное решение задачи (1-1.1) с начальным условием x(tk) = #ь полученное за один шаг размера 7 и x(tf. + 1)- численное решение задачи (1.1.1) с начальным условием x(tk) = Xk, полученное за два шага размера т /2 на отрезке
Тогда для оценки главного члена локальной ошибки при достаточно малом 7 справедливо (см., например, [1, 5, 43]): +1 = 2 ї к+1 +о+м+2)- с1-4-3) Впервые оценка (1.4.3) была выведена Ричардсоном в [131] и [132].
В итоге, для вычисления главного члена локальной ошибки правило Рунге требует трех шагов численного метода, что может служить серьезным препятствием для его использования в том случае, когда накладные расходы на один шаг численного интегрирования велики.
Другой подход к оценке локальной ошибки основан на использовании численных методов разных порядков. Здесь для получения оценки проводят два интегрирования из точки Ц с шагом 7 методами, имеющими разные порядки. Пусть мы для определенности применяем методы порядков s и р, причем s р. Предположим, что приближенное решение Xk+i вычислено методом порядка s, a x +i — методом порядка р. Допустим также, что правая часть задачи (1.1.1) / + 1 непрерывно дифференцируема в некоторой окрестности точного решения этой задачи.
Тогда, используя асимптотическое разложение локальной ошибки (1.4.2) для каждого из двух численных методов, легко видеть, что локальная ошибка метода порядка s удовлетворяет равенству Ахш = xk+i - Xjfc+i + 0(т-+1). (1.4.4)
На практике чаще всего ограничиваются оценкой только главного члена локальной ошибки, хотя изложенный выше способ позволяет делать и более точные оценки. Заметим, что способ, опирающийся на два метода разных порядков, требует только двух интегрирований исходной задачи (1.1.1) на каждом шаге, что в некоторых случаях дает заметный выигрыш во времени по сравнению с правилом Рунге.
Примерно последние 40 лет предпринималось немало усилий с целью повышения эффективности именно второго способа оценки локальной ошибки. Для явных одношаговых методов это привело к созданию вложенных формул Рунге-Кутты (см., например, [43, 66, 75, 82, 86, 87, 97, 126, 136, 148]), где для интегрирования методами порядков s и s + 1 применяют одни и те же стадийные величины, что значительно снижает затраты машинного времени.
Покажем, как использовать оценки (1.4.3) и (1.4.4) для управления размером шага на практике. Обычно задают некоторый допуск ел и требуют, чтобы на каждом шаге интегрирования выполнялось неравенство ЛхА+іЦ єл. (1.4.5)
Из сказанного выше вытекает, что осуществление условия (1.4.5) всегда можно обеспечить надлежащим выбором размера шага интегрирования г . Типичную процедуру управления локальной ошибкой при численном интегрировании задачи (1.1.1) можно записать в виде следующего алгоритма. Пусть в к-оїі точке сетки нам известны величины Xk — x(tk) и 7. Тогда а) вычисляем величины Xk+\ и Axf.+1 i б) проверяем справедливость условия (1.4.5) для Дя +і; в) если оценка локальной ошибки Ах +і не удовлетворяет (1.4.5), то отбрасываем величину #&+і, уменьшаем шаг интегрирования 7 и переходим к пункту а); г) если условие (1.4.5) выполнено, то принимаем величину Xk+i в каче стве приближенного решения задачи (1.1.1) в точке tk+i, вычисляем размер нового шага интегрирования Тк+\ и переходим к пункту а) в следующей точке сетки гит.
В качестве размера шага интегрирования 7 выбирают максимальную величину, при которой выполняется условие (1.4.5). Так как на практике обычно ограничиваются оценкой только главного члена локальной погрешности, условие (1.4.5) в этом случае принимает вид Шік)И+і ел. (1.4.6) Пусть на к + 1-ом шаге соотношение (1.4.6) не выполнено. Тогда для вычисления оптимальной длины шага приходим к следующему простому выражению (см., например,[65]): Чет) (1А7)
Если же на к + 1-ом шаге неравенство (1.4.6) имеет место, то для вычисления нового шага интегрирования также используют формулу (1.4.7), но с заменой т на т +і Как мы видим, формула (1-4.7) для вычисления оптимального размера шага интегрирования не учитывает влияние старших членов в разложении локальной ошибки (1.1.7). Поэтому этот факт, а также ошибки округления могут приводить к значительному количеству отброшенных шагов и, следовательно, дополнительным затратам машинного времени. Для преодоления этого недостатка обычно вводят гарантийный множитель 7 Є (0.5,0.9) и формула (1-4.7) для практического вычисления размера шага интегрирования в этом случае принимает вид
Формулировка комбинированных Е-методов
Итак, в каждой точке сетки (т.е. при фиксированном k = 0,1,... , Л" — 1) рассмотрим следующую задачу: Хк+і = ОгкХш (2.4.1) относительно вектора неизвестных Xk+i = ((sifc+1/2) Лхк+і) ) Є R m, где отображение Gk : R2m —» R2m задано формулой І , г { [г) (г) , (г) (г) Л . (0) \ к + г J2 т \а\ 9к + Ч Якії) + та2 9к+і/2і лгу, _ г=0 . + г Е т (Ь Л,« + іВДЛ + г6 V1/2 , \ г=0 ч J!
В общем случае задача (2.4.1) не разрешима аналитически и, таким образом, требует привлечения того или иного итерационного процесса.
Обозначим через xk+i = xk+i(N) значение приближенного решения задачи (2.4.1) в точке tk+u полученное после N итераций некоторого итерационного процесса. Тогда, применяя для решения (2.4.1) простые или ньютоновские итерации (с тривиальным или нетривиальным предиктором), как в [21] или [15, 16], приходим к трем комбинированным алгоритмам. Е-метод с простыми итерациями (ЕПИ-метод): ХІ+1 = GlXl-\, і = 1,2,..., N, (2.4.2a) ХШ = Х?, к = 0,1,...,К -1, (2.4.26) где отображение G% означает тоже самое, что и в формуле (2.4.1), но с заменой хк на хк = хк. Последнее равенство подразумевает, что в качестве приближенного решения задачи (2.4.1) на к-ом шаге мы берем последние т координат вектора Xj?, полученного после N итераций процесса (2.4.2а). Е-метод с итерациями Ньютона (ЕН-метод): Xl+1 = Xt+\ - dFaxt )-1 F[Xt+\, і = 1,2,..., N, (2.4.3а) Хк+1 = Xj t к = 0,1,..., А - 1, (2.4.36) -о -к+І где Fk Е т — GTk (Еът — единичный оператор в R2m) и дРк(Хк+\) есть матрица Якоби отображения AJ, вычисленная в точке Хк _\. Е-метод с модифицированными итерациями Ньютона (ЕМН-метод): Х к+1 = ХД - dmX»k+1)-lFXRu і = 1,2,..., ЛГ, (2.4.4а) Хш = X", А = 0,1,..., А - 1, (2.4.46)
Итак, мы построили три комбинированных Е-метода с различными итерационными процессами. Понаблюдаем теперь на практике какое влияние оказывает количество итераций на скорость сходимости базового Е-метода. В качестве теста выберем сначала задачу из [43, с. 146]: x[(t) = 2tx2(t)l/5x4{t), x 2(t) = 10t ехр{5(я3(0 - l)}z4(t), (2.4.5a) x t) = 2ts4(t), x 4{t) = -2( ln{xi(t)}, (2.4.56) с начальными условиями ж,-(0) = 1, г — 1,2,3,4. Эта задача имеет точное решение xi(t) = exp{sin(f2)}, x2{t) = cxp{5sin(f2)}, я3(і) = sin(2) + 1, xA{t) = cos{t2). [ } Кроме того, рассмотрим также реальную задачу из небесной механики, т.е. ограниченную задачу трех тел: х[ — хз, х 2 = Х4 (2.4.7а) х4 = х2- 2х3 - (/2 - / 2т; .2 , 2.з/2- (2"4-7в) («і + / 2)2 + жі) ((жі - т)2 + ) Если в качестве начальных условий задачи (2.4.7) взять a;i(0) = 0.994, х2(0) = 0, а;3(0) = 0, х4(0) = -2.00158510637908252240, где щ = 1 - д2 и //2 = 0.012277471, то решение будет периодическим с периодом Т = 17.065216560157962558891.
Проведем интегрирование задач (2.4.5) и (2.4.7) всеми тремя комбинированными методами (2.4.2)-(2.4.4) на отрезке [0,10] для первой задачи (или [0,Х] для второй) при различных соотношениях числа итераций N и размера шага т (для задачи (2.4.7) данные приведены в зависимости от числа точек сетки К). Вычислим глобальные ошибки каждого из методов с помощью формул для точного решения (2.4.6), или используя периодичность решения задачи (2.4.7), и определим тем самым реальный порядок Е-методов. В качестве основы для комбинированных алгоритмов возьмем Е-методы шестого и восьмого порядков, т.е. (2.1.6) и (2.2.9) соответственно.
Результаты эксперимента соберем в таблицах 2.1-2.12. В табл. 2.1-2.3 приведены данные применения алгоритмов (2.4.2)-(2.4.4) для первой задачи и метода шестого порядка, затем - для первой задачи и метода восьмого порядка. Аналогично, в табл. 2.7-2.12 представлены результаты экспериментов для второй задачи. Здесь и далее прочерки означают расходимость метода при соответствующих значениях параметров N и т. Кроме того, надо иметь в виду, что ошибки округления оказывают существенное вляние на точность вычислений при высокой скорости сходимости и малом шаге.
Из представленных данных можно сделать следующие выводы. Во-первых, порядок сходимости Е-методов действительно зависит от количества итераций в итерационном процессе, который используется для разрешения соответствующей нелинейной системы. Во-вторых, теоретический (в данном случае шестой или восьмой) порядок сходимости достигается только при достаточно большом числе итераций Б точке сетки, причем для каждого из рассмотренных выше итерационных процессов это число свое.
Таким образом, крайне важно теоретически определить порядок методов (2.4.2)-(2.4.4) и найти тем самым оптимальное число итераций для каждого комбинированного алгоритма, что является предметом исследования в следующем параграфе. Здесь же только отметим, что доказанная в параграфе 2.3 Л-устойчивость Е-методов (2.2.2) естественным образом переносится на ЕЇІИ-, ЕН- и ЕМН-метод при достаточно большом числе итераций, а значит, построенные алгоритмы (2.4.2)-(2.4.4) приспособлены для решения жестких задач.
Функции и структуры для численного решения обыкновенных дифференциальных уравнений
В заголовочном файле numerics.h определены вспомогательные функции и базовые типы, которые применяются во всех функциях-решателях библиотеки. Подробное описание указанного файла, используемого в версии пакета программ ИНТЕГРАТОР 1.0.
В силу специфики построения методов со старшими производными в заголовочном файле numerics.h были модифицирован тип, используемый для описания функции, стоящей в правой части уравнения (1.1.1), и ее производных g (x,t) typedef dblmatrix ( rhfunc )(const dblmatrixfe, long double, int); Для задания матрицы Яко б и правой части и ее производных определен следующий тип: typedef dblmatrix ( jacobian)(const dblmatrix& , long double, int);
Таким образом пользователь должен будет обеспечить описание интегрируемой задачи в приведенном выше формате. Функция типа rhf unc должна возвращать матрицу содержащую по столбцам от нулевого до заданного в параметре типа int значение соотвсствующих производных правой части задачи (1.1.1). Функция типа Jacob і an должна возвращать соответственно ряд состыкованных матриц Якоби производных правой части задачи (1.1.1) от нулевого до заданного в параметре типа int. struct odestruct С odestruct(rhfunc _RHF, jacobian _JAC, const dblmatrixfe _X0, long double _T0, long double _TF); rhfunc rhf; jacobian jac; dblmatrix xstart; long double tstart; long double tfinal; ; Описание:
Структура odestruct содержит пять полей данных для представления задачи Коши для системы обыкновенных дифференциальных уравнений вида (1.1). rhf данное поле служит для задания указателя на функцию вы числения правой части уравнения (1.1.1) и ее производных (значение по умолчанию 0); jac данное поле служит для задания указателя на функцию вы числения матриц Якоби правой части уравнения (1.1.1) и ее старших производных (значение по умолчанию 0); m xstart данное поле задает вектор столбец начальных значений XQ (значение по умолчанию dblamtrix(0,0) - пустая матрица); tstart данное поле задаст начало отрезка интегрирования t.Q (значение по умолчанию 0); tfinal данное поле задает конец отрезка интегрирования Q + Т (значение по умолчанию 0).
Для создания структуры odestruct можно использовать либо конструктор по умолчанию, либо конструктор odestruct::odestruct(rhfunc _RHF, jacobian _JAC, const dblmatrixfe _X0, long double _T0, long double _TF);
Первый конструктор создает пустую структуру odestruct, а второй создает структуру odestruct и копирует содержимое передаваемых параметров в соответствующие поля. Структура solveropt Определение: struct solveropt int order; mttype method; int iter; long double TOL; long double gammal; long double InitH; long double MaxH; int disp; int freq; int stnum; ostream tOUTstream; ostream xOUTstream; ostream eOUTstream; ; Описание: Структура solveropt служит для представления опций функции-решателя и содержит четыре группы таких опций. Первая группа опций предназначена для выбора порядка метода со старшими производными, метода решения систем нелинейных уравне ний и установки числа итераций для этого метода. order данное поле задает порядок метода со старшими производ ными (значение по умолчанию 4); method данное поле задает метод решения систем нелинейных уравнений; пакет ИНТЕГРАТОР 1.1 поддерживает че тыре метода: solveropt: :mt_NEWT0N — метод Ньютона, solveropt: :mt_mKEWTOK — модифицированный метод Нью тона и solveropt: :mt_sMEWTON — обобщенные итерации Ньютона и solveropt: :mt_sITER — метод простой итера ции (значение по умолчанию solveropt: :mt_mNEWTON); iter данное поле определяет количество итераций для итераци онного метода, заданного полем method (значение по умолчанию 1).
