Козлов Олег Степанович
Скворцов Леонид Маркович
Ходаковский Виктор Владимирович

Решение дифференциальных и дифференциально-алгебраических уравнений в программном комплексе «МВТУ»

Аннотация: Рассматриваются методы численного решения обыкновенных дифференциальных уравнений (ОДУ) и дифференциально-алгебраических уравнений (ДАУ), реализованные в программном комплексе (ПК) «МВТУ». Обсуждаются вопросы выбора метода и задания его параметров при решении задач различного типа. Описаны новые явные и неявные методы, эффективные для жестких задач. Один из неявных методов позволяет также решать системы ДАУ индексов 2 и 3. Приведены результаты решения тестовых задач в сравнении с известными методами

  1. Введение
  2. Общие положения
    1. Постановка задачи
    2. Параметры интегрирования
    3. Управление шагом интегрирования
  3. Явные методы
    1. Классические методы
    2. Модифицированные классические методы
    3. Адаптивные методы
  4. Неявные методы
    1. Метод Гира
    2. Адаптивный метод
    3. Диагонально неявные методы Рунге-Кутты
  5. Примеры решения тестовых задач
    1. Жесткие задачи
    2. Локально-неустойчивые задачи
  6. Решение дифференциально-алгебраических уравнений (ДАУ)
    1. Метод пространства состояний
    2. Метод ε-вложения
    3. Системы ДАУ высших индексов
  7. Приложение
  8. Литература и Интернет

Введение

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

Таким образом, для профессиональной работы с любым моделирующим ПК пользователь должен обладать некоторыми знаниями о реализованных в нем численных методах. Цель настоящей статьи – дать пользователю необходимые сведения о решателях ПК «МВТУ». Статья может быть полезной для пользователей и разработчиков других аналогичных ПК.

1. Общие положения

1.1. Постановка задачи

Моделирование процессов в непрерывных динамических системах сводится к численному решению задачи Коши

(1.1)

x' = f(t, x),   x(t0) = x0,   t0t ≤ (t0 + T),

где x ­– n-мерный вектор переменных состояния, fn-мерная вектор-функция правых частей, t – независимая переменная, T – величина интервала интегрирования. Численное решение задачи (1.1) получаем в виде последовательности векторов xm, m = 1, 2, … аппроксимирующих истинное решение на сетке , где hi – величина i-го шага интегрирования.

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

Практика показала, что реальные процессы, как правило, описываются жесткими системами ОДУ. Достаточно часто в технических приложениях встречаются быстро осциллирующие системы, описывающие высокочастотные колебания, а также локально-неустойчивые системы, в решении которых кратковременные участки с расходящимся процессом чередуются с более продолжительными стабильными участками. Перечисленные типы задач предъявляют совершенно разные требования к методам интегрирования. При интегрировании жестких систем необходимо обеспечить быстрое затухание жестких составляющих, для этого используют неявные A- или L‑устойчивые методы [4, 6]. Такие методы подавляют все составляющие решения, соответствующие большим по модулю собственным значениям (если только шаг не выбран очень малым), поэтому они плохо приспособлены для решения осциллирующих и локально-неустойчивых систем. Для интегрирования осциллирующих систем следует применять специальные методы, обеспечивающие правильный характер огибающей колебательного решения. Специальные методы следует использовать также и для решения жестких локально-неустойчивых систем.

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

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

1.2. Параметры интегрирования

Выбор метода и установка его параметров в ПК «МВТУ» осуществляется в диалоговом окне «Параметры расчета». Это окно имеет 6 закладок, но только две из них (Основные и Дополнительные) имеют отношение к методам интегрирования, а остальные позволяют устанавливать различные режимы моделирования.

При открытии окна «Параметры расчета» активной является закладка Основные, где задаются следующие параметры:

Метод интегрирования
Можно выбрать один из 16 методов (Эйлера, Рунге-Кутты классический, Рунге-Кутты модифицированный, Мерсона классический, Мерсона модифицированный, Адаптивный 1, …, Адаптивный 5, Адаптивный неявный, Диагонально неявный, Гира, Неявный Эйлера, DIRK33, DIRK44). По умолчанию используется Адаптивный 1. Метод Эйлера может иметь только фиксированный шаг, а все остальные методы могут использоваться как с фиксированным, так и с переменным шагом.
Время интегрирования T
Длительность интервала модельного времени, на котором выполняется моделирование (по умолчанию ­10).
Минимальный шаг интегрирования hmin
Этот параметр ограничивает снизу шаг интегрирования и по умолчанию равен 10-6. Если значение hmin оказывается недостаточно малым для расчета с заданной точностью, то выдается сообщение «Заданная точность не обеспечивается». В таких случаях следует уменьшить минимальный шаг либо снизить требования к точности. Уменьшение значения hmin обычно мало сказывается на времени счета, однако нежелательно задавать этот параметр меньше 10-15 T.
Максимальный шаг интегрирования hmax
Ограничение сверху на шаг интегрирования, а для метода Эйлера – величина шага (по умолчанию 0.1). Если задать hmax = hmin, то интегрирование будет выполняться с фиксированным шагом. Слишком малое значение hmax приводит к неоправданному увеличению времени счета, а большое значение может привести к уменьшению числа точек, выводимых на график. Рекомендуемое значение – T /100.
Шаг вывода результатов dT
Приблизительное значение шага вывода результатов в графические окна, при этом реальный шаг вывода не может быть меньше шага интегрирования. Точное значение заданного шага вывода будет соблюдаться, если установить (на закладке Дополнительные) табуляцию результатов расчета. Рекомендуется задавать dT = (10-3…10-2T. При dT = 0 результаты выводятся на каждом шаге интегрирования. Аналогичного эффекта для отдельного графического окна можно добиться, установив в окне свойств соответствующего блока: Вывод на каждом шаге ­– Да; Прореживание точек – Нет.
Точность Rtol
Допустимая относительная ошибка интегрирования (по умолчанию 0.001). Не рекомендуется задавать это значение больше 0.01 или меньше 10-10.

На закладке Дополнительные устанавливаются следующие параметры:

Значения выходов при инициализации
Начальные значения алгебраических переменных при наличии в системе алгебраических контуров (по умолчанию 10-10). Для каждой переменной можно установить свое значение, если включить в контур блок Y = F(Y).
Абсолютная ошибка Atol
Допустимая абсолютная ошибка интегрирования (по умолчанию 10-10). Ненулевое значение Atol предотвращает неоправданное уменьшение шага в тех случаях, когда значения некоторых переменных приближаются к нулю.
Максимальное значение производных (по умолчанию 10300)
Если производная некоторой переменной превышает по модулю это значение или произошла фатальная ошибка (переполнение, деление на 0 и т.д.), то шаг считается неудачным. После этого интегрирование возобновляется с уменьшенным шагом либо прекращается с выдачей сообщения об ошибке, если размер шага достиг hmin.
Метод итерирования для петель
Метод решения алгебраических уравнений при наличии в системе алгебраических контуров (Простая итерация – по умолчанию, Ньютона-Рафсона, Бройдена (секущих), Без итераций). Выбранный метод используется для расчета начального состояния алгебраических переменных (независимо от метода интегрирования), а также для расчета алгебраических переменных в процессе интегрирования явным методом. В процессе интегрирования неявным методом дифференциальные и алгебраические переменные решаются совместно, поэтому выбор метода итерирования не имеет значения. Наиболее надежным является метод Ньютона-Рафсона, но в некоторых случаях и другие методы могут иметь преимущество.
Число итераций
Максимальное число итераций при решении алгебраических уравнений (по умолчанию 20). Этот параметр, как и предыдущий, влияет на решение только в тех случаях, когда в системе есть алгебраические контуры или блоки Y = F(Y), F(Y) = 0.
Табуляция результатов расчета
Установив этот флаг, мы добьемся вывода результатов на графики с шагом, в точности равным заданному шагу вывода. Этот флаг рекомендуется также устанавливать при наличии в системе дискретных блоков, тогда шаг вывода следует устанавливать таким, чтобы периоды квантования дискретных блоков были равны или кратны ему. Если табуляция установлена, то шаг интегрирования не может превышать шага вывода результатов.

1.3. Управление шагом интегрирования

Начальный шаг устанавливается равным . В процессе интегрирования с переменным шагом необходимо, кроме решения в очередной точке, вычислять оценку ошибки, которая используется для управления величиной шага. Для этого применяют две различные формулы интегрирования, дающие на m-м шаге два решения: xm и . Пусть x = xm – вектор переменных, Δx = xm - xm-1 – приращение, δx = xm -  – оценка ошибки на последнем шаге. Для управления величиной шага используется нормированная ошибка, вычисляемая по формуле

(1.2)

,

причем шаг считается удачным, если err ≤ 1. В случае удачного шага принимаем значение xm в качестве нового вектора переменных. При неудачном шаге производится пересчет с уменьшенным размером шага. Шаг считается неудачным также и в том случае, когда одна из производных превысила максимально допустимое значение, задаваемое в окне «Параметры расчета», либо произошло прерывание, вызванное переполнением, делением на ноль, недопустимым значением аргумента и т.п. В таких случаях размер шага уменьшается сразу в 4 раза, но если он стал меньше hmin, то моделирование прекращается и выдается сообщение об ошибке.

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

В большинстве методов, реализованных в ПК «МВТУ», используется стандартная процедура управления величиной шага [6], задаваемая формулой

(1.3)

hm+1 = fachmerr,

где fac = 0.8…0.9 – множитель безопасности, a – величина, обратная порядку оценки ошибки, err – нормированная ошибка (1.2).

2. Явные методы

Явные методы представлены в ПК «МВТУ» методом Эйлера, классическими и модифицированными методами Рунге-Кутты и Мерсона, а также пятью адаптивными методами.

2.1. Классические методы

Метод Эйлера является простейшим, имеет 1-й порядок и задается расчетной формулой

xm+1 = xm + h f(tm, xm).

Его можно использовать при интегрировании с постоянным шагом и низкой точностью.

Классический метод Рунге-Кутты задается формулами

Он имеет 4‑й порядок, однако, применяя экстраполяцию по Ричардсону, можно повысить порядок до 5‑го, получив дополнительно оценку ошибки. Окончательные расчетные формулы имеют вид

,

где δxm+1 – оценка ошибки, а xm+1(h) и xm+1(h / 2) – численные решения, полученные, соответственно, путем выполнения одного шага размера h и двух шагов размера h / 2. Отметим, что такой метод требует довольно значительных вычислительных затрат, поскольку на каждом шаге функция f вычисляется 11 раз.

Метод Мерсона имеет 4-й порядок и задается формулами

По сравнению с классическим методом Рунге-Кутты, в методе Мерсона используется более экономный способ оценивания ошибки и управления величиной шага.

При решении умеренно жестких задач явными методами размер шага обычно приближается к границе устойчивости, но в этом случае процедура управления шагом может оказаться неустойчивой. В результате происходят резкие колебания величины шага, приводящие к снижению точности и увеличению вычислительных затрат. Устойчивость управления шагом (SC‑устойчивость [6]) можно обеспечить при использовании процедуры ПИ-управления [6], задаваемой формулой

(2.1)

hm+1 = fachm ∙ (errm) ∙ (errm-1)β,

где коэффициент β выбирается из условия SC-устойчивости (обычно 0 ≤ β ≤ 0.1). При β = 0 получаем стандартную процедуру (1.3). Для метода Мерсона выбраны значения α = 0.2, β = 0.08, обеспечивающие SC-устойчивость. Эксперименты показали, что использование ПИ-закона (2.1) позволяет снизить затраты и повысить точность решения умеренно жестких задач по сравнению с использованием стандартной процедуры (1.3).

Классические методы рекомендуется использовать для решения нежестких задач.

2.2. Модифицированные классические методы

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

Модифицированный метод Рунге-Кутты отличается от классического способом вычисления k4

Здесь z – вектор покомпонентных оценок наибольшего по модулю собственного значения матрицы hJ, где J – матрица Якоби. В приведенных формулах покомпонентно выполняются не только арифметические, но и логические операции. При проверке выполнения неравенств компоненты вектора z не вычисляются, а вместо этого вычисляются и сравниваются компоненты векторов k2 - k1 и k3 - k2. Для нежестких компонент, т.е. при выполнении условия zi ≥ -2 применяется классическая формула, а для жестких компонент – формула, полученная из условия стабилизации метода в полученной точке жесткого спектра, причем вычисляются и используются только компоненты вектора z-1, удовлетворяющие условию -0.5 < zi-1 < 0.

Модифицированный метод Мерсона также отличается от классического способом вычисления k4:

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

2.3. Адаптивные методы

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

Одношаговые адаптивные методы строятся на основе стадий Рунге-Кутты, которые выполняются по формулам

(2.2)

где s – число стадий, b и a – параметры метода (в общем случае самонастраиваемые [7, 8]). Далее вычисляются векторы

(2.3)

,

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

Вектор покомпонентных оценок наибольшего собственного значения получим в виде

(2.4)

z = ys / ys-1.

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

(2.5)

.

Из (2.5) видно, что порядок метода для жестких задач не может превышать s - 2. На основе формул (2.2)–(2.5) построены одношаговые методы Адаптивный 1, 2, 3, 5.

Адаптивный 1 – трехстадийный метод (s = 3), имеет 1-й порядок для жестких и 2-й для нежестких задач. В нем β = 1, тогда вектор (2.4) получим в виде z = (k3 - k2) / [α (k2 - k1)]. Компоненты вектора настраиваемых параметров вычисляются по формуле

,

а формула (2.5) принимает вид xm+1 = xm + h k1 + h c (k2 - k1). Ошибка оценивается с помощью двухшаговой формулы

где ω = hm+1 / hm – отношение текущего шага к предыдущему, di = -zi ci, i = 1, …, n. Для очень жестких компонент выполняются одна или две дополнительные стадии, расширяющие область устойчивости в окрестности точек жесткого спектра.

Адаптивный 2 – расчетная схема аналогична методу Адаптивный 1, но дополнительные стадии не используются, а значение b изменяется от 2/3 до 1 в зависимости от жесткости задачи, благодаря чему порядок метода для нежестких задач повышается до 3-го. Ошибка оценивается по правилу Рунге, т.е. с использованием одного шага размера h и двух шагов размера h / 2.

Адаптивный 3 – четырехстадийный метод (s = 4) 2-го порядка для жестких и 3-го для нежестких задач. Вектор настраиваемых параметров вычисляется по формуле

.

Для оценивания ошибки используется вложенная формула.

Адаптивный 5 – пятистадийный метод 2-го порядка для жестких и 3-го для нежестких задач, шаг которого выполняется по формуле

xm+1 = xm + h (y1 + d2 y2 + d3 y3).

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

Одношаговые адаптивные методы имеют невысокую точность, поэтому в ПК «МВТУ» реализован также многошаговый метод. На многих жестких задачах он оказался наиболее эффективным среди всех адаптивных методов. Поскольку этот метод нигде не описан, рассмотрим его более подробно.

Адаптивный 4 – многошаговый адаптивный метод переменного порядка (1…5 для жестких и 2…6 для нежестких задач), основанный на использовании разделенных разностей, определяемых рекуррентно согласно формулам

.

(аналогично определяются ). Введем также коэффициенты aijgi, зависящие от размеров текущего и предыдущих шагов и задаваемые формулами

(2.6)

В этом случае k-шаговые формулы прогноза и коррекции при интегрировании с переменным шагом запишутся в виде

(2.7а)

,

(2.7б)

,

где – разделенная разность с использованием вместо fm+1. Формула прогноза имеет порядок k, а порядок формулы коррекции на 1 больше.

На основе формул вида (2.7) обычно строят многошаговые методы Адамса переменного порядка и шага [4, 5], при этом , а оценка ошибки находится как разность между прогнозом и скорректированным решением. Для адаптивного метода эти формулы являются предварительными, на их основе получаем вектор оценок наибольшего собственного значения . Вместо него вычисляем вектор

(2.8)

,

причем только для жестких и неустойчивых компонент, удовлетворяющих условию , где m - константа, определяемая областью устойчивости и зависящая от k. Формула интегрирования для таких компонент имеет вид

(2.9)

,

где векторные коэффициенты di вычисляются по формулам

(2.10)

(en-мерный вектор с единичными компонентами). Формула (2.9) построена таким образом, что она имеет порядок k и обеспечивает затухание за один шаг при решении уравнения x' = λ x, λ < 0.

Формулы (2.6)–(2.10) задают многошаговый адаптивный метод с переменным шагом, при этом для нежестких компонент ограничиваемся формулами (2.7), а дальнейшие вычисления не производим. Оценку ошибки получаем как приращение в формуле коррекции (2.7б) для нежестких компонент или последний член в формуле (2.9) для жестких компонент. Приведенные формулы позволяют также легко изменять порядок метода, который определяется значением k. Если на предыдущем шаге k = kold, то на очередном шаге значение k может быть от 1 до kold + 1. Новое значение k определяется на этапе прогноза исходя из условия, чтобы прогноз был наиболее точным. Формулу прогноза (2.7а) можно записать в виде

.

Новое значение k определяется как максимальное число, для которого выполняется условие

|| δ1 || > || δ2 || > … > || δk||,

где норма ошибки вычисляется в соответствии с (1.2).

3. Неявные методы

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

3.1. Метод Гира

Метод Гира – многошаговый метод переменного порядка (от 1-го до 6-го), построенный на основе формул дифференцирования назад (ФДН) [4–6] вида

,

где коэффициенты a1, …, akb зависят от порядка формулы k. В методе Гира эти формулы реализованы с помощью вектора Нордсика [4, 5], позволяющего использовать эффективный алгоритм автоматического изменения порядка и величины шага. Последние компоненты вектора Нордсика позволяют получить оценку ошибки. При k = 1 (тогда a1 = b = 1) получим неявный Эйлера – метод 1-го порядка, реализованный в ПК «МВТУ».

3.2. Адаптивный метод

Адаптивный неявный метод имеет порядок от 2-го для жестких до 4-го для нежестких задач; он построен на основе метода трапеций, формула которого имеет вид

Выполняя один шаг величиной h, получим xm+1(h), а выполняя два шага величиной h / 2, получим xm+1(h / 2). Далее вычисляем вектор покомпонентных оценок наибольшего собственного значения

,

который используется в заключительной расчетной формуле

3.3. Диагонально неявные методы Рунге-Кутты

В общем случае метод Рунге-Кутты задается формулами

(3.1)

и может быть представлен в виде таблицы Бутчера

(3.2)

.

Часто приводят также формулу

,

которая используется для получения оценки погрешности численного решения . В этом случае в таблицу (3.2) добавляется строка коэффициентов .

В явном методе aij = 0 при j ≥ i, тогда формулы (3.1) задают расчетный алгоритм, который может быть непосредственно реализован. В противном случае метод является неявным и требует при своей реализации решения системы алгебраических уравнений. Среди неявных методов Рунге‑Кутты наиболее просто реализуются диагонально неявные (DIRK – Diagonally Implicit Runge‑Kutta), у которых матрица A имеет нижнюю треугольную форму.

В ПК «МВТУ» реализованы три метода DIRK, имеющие таблицу Бутчера вида

(3.3)

Эти методы имеют явную первую стадию и r = s -1 неявных стадий с одинаковыми диагональными элементами матрицы A, поэтому их обычно называют ESDIRK (Explicit first stage Singly DIRK). Явная стадия не требует вычислений, поскольку результат ее выполнения совпадает с последней стадией предыдущего шага. Такие методы называют также FSAL (First Same As Last). Преимущество методов вида (3.3) по сравнению с другими методами DIRK состоит в том, что они являются жестко точными и могут иметь второй стадийный порядок [9, 10].

Диагонально неявный – метод 2-го порядка, задаваемый таблицей Бутчера

.

Его можно интерпретировать как последовательное применение правила трапеций и формулы дифференцирования назад 2-го порядка, поэтому он получил название TR-BDF2. Этот метод реализован также и в системе MATLAB+Simulink под названием ode23tb.

DIRK33 – метод 3-го порядка с 3-мя неявными стадиями и коэффициентами

DIRK44 – метод 4-го порядка с 4-мя неявными стадиями и коэффициентами

Методы DIRK33, DIRK44 и еще три аналогичных метода 4-го и 5-го порядков предложены в [9], где изложены также детали их реализации.

4. Примеры решения тестовых задач

Решатели ОДУ ПК «МВТУ» были испытаны на множестве разнообразных задач, среди которых – задачи из тестовых наборов, приведенных в [5, 6, 11]. Результаты тестовых испытаний в сравнении с решателями системы MATLAB представлены в [12] и в Приложении. Среди задач были нежесткие, жесткие, локально-неустойчивые, осциллирующие, а также задачи с разрывами. Оценивались точность решения и вычислительные затраты. Нежесткие задачи оказались самыми легкими и были успешно решены всеми методами с автоматическим выбором шага. Поэтому остановимся на жестких задачах, численное решение которых часто встречает затруднения.

4.1. Жесткие задачи

Посмотрим, как влияет жесткость задачи на точность ее решения и вычислительные затраты. Будем решать задачу Капса

x'1 = -(μ + 2) x1 + μ x22,   x1(0) = 1,
x'2 = x1 - x2 - x22,   x2(0) = 1,   0 ≤ t ≤ 1,

которая имеет плавное решение x1(t) = exp(-2t), x2(t) = exp(-t), не зависящее от параметра жесткости m (собственные значения якобиана при больших m примерно равны – m, –1). В качестве меры вычислительной работы будем использовать число вычислений правой части Nf (включая также и вычисления, выполняемые при расчете матрицы Якоби), а точность численного решения будем оценивать по формуле

(4.1)

scd = – log10 (error),

где error – максимальная относительная ошибка на всем интервале интегрирования. Таким образом, scd – минимальное число правильных значащих цифр численного решения (significant correct digits). Результаты для явных методов и трех лучших неявных методов при hmin = 10-6, hmax = 1, Rtol = 10-3, Atol = 10-10 приведены в табл.1.

Таблица 1. Точность и вычислительные затраты при решении задачи Капса

Метод m = 101
scd Nf
m = 102
scd Nf
m = 103
scd Nf
m = 104
scd Nf
m = 105
scd Nf
Рунге-Кутты кл.
Рунге-Кутты мод.
Мерсона
Мерсона мод.
3.7   110
3.7   110
4.9   100
4.9   100
4.1   264
3.6   198
4.5   250
4.5   250
0.4   2639
3.3   528
3.8   1419
4.4   835
-0.7   25770
3.1   605
3.9   14098
3.9   4133
-1.9   255720
3.4   1210
4.1   140902
3.1   19715
Адаптивный 1
Адаптивный 2
Адаптивный 3
Адаптивный 4
Адаптивный 5
4.0   87
3.8   120
3.7   75
5.0   42
4.3   120
3.6   143
3.7   208
3.5   125
5.6   48
3.6   217
2.7   237
2.2   248
3.4   150
5.6   51
3.9   260
2.1   304
2.0   256
3.4   165
5.5   60
3.9   270
2.0   319
1.8   240
3.6   195
5.5   60
3.9   270
Гира
DIRK33
DIRK44
2.9   39
2.6   36
4.0   55
3.1   40
3.0   49
4.0   53
3.1   40
3.3   34
4.2   47
3.0   40
3.3   34
3.6   41
3.0   40
3.3   34
3.6   41

Вычислительных затраты классических явных методов возрастают практически линейно при увеличении жесткости, что вызвано необходимостью уменьшать шаг интегрирования для обеспечения устойчивости численного решения. Значительное снижение точности классического метода Рунге‑Кутты объясняется отсутствием SC‑устойчивости. Модифицированные методы эффективнее классических, но заметно уступают адаптивным и неявным методам. Адаптивные методы демонстрируют небольшое увеличение объема вычислений при повышении жесткости (а первые два из них также и снижение точности), что вызвано снижением порядка при решении жестких задач. Точность и вычислительные затраты неявных методов практически не зависят от жесткости. Результаты решения задачи Капса подтвердили эффективность метода Гира, который традиционно считается одним из лучших для жестких задач, а также показали перспективность новых методов (адаптивных, DIRK33, DIRK44).

Для более полного тестирования жестких решателей ПК «МВТУ» и системы MATLAB были выбраны первые шесть задач из тестового набора, приведенного в [6] (VDPOL, ROBER, OREGO, HIRES, E5, PLATE). Результаты приведены в [12] и в Приложении. С приемлемыми затратами и точностью решили все шесть задач четыре решателя ПК «МВТУ» (Диагонально неявный, Гира, DIRK33, DIRK44) и два решателя системы MATLAB (ode15s и ode23tb). В [12] приведен график зависимости суммарных вычислительных затрат от усредненной по всем задачам точности для этих решателей. Лучшие результаты – у методов Гира, DIRK33, DIRK44 и ode15s, при этом методы Гира и ode15s имеют небольшое преимущество при расчетах с высокой точностью, а решатели DIRK33 и DIRK44 при той же задаваемой точности обеспечивают, как правило, более высокую точность. Отметим, что явный решатель Адаптивный 4 оказался лучшим среди всех методов при решении жесткого уравнения Ван-дер-Поля (задача VDPOL).

Приведем также результаты сравнения решателей ПК «МВТУ» с решателем RADAU, реализующим методы Radau IIa порядков 5, 9 и 13 [6, 11]. По результатам тестирования наиболее известных программ этот решатель был лучшим для большинства задач из [11]. Мы выбрали задачи VDPOL и OREGO, которые решали при Rtol = Atol = 10-4 согласно условиям, приведенным в [11]. Из этой же работы заимствованы результаты для RADAU (tables II.8.4, II.9.3 в [11]). Полученные результаты представлены в табл. 2. Вычислительные затраты оценивались следующими показателями: Nf ­– число вычислений функции (сюда не включены вычисления, выполняемые при расчете якобиана); NJac – число вычислений якобиана; NLU – число LU‑разложений. При вычислении scd по формуле (4.1) использовалась максимальная относительная ошибка в конце интервала интегрирования. Отметим высокую эффективность явного адаптивного метода, который показывает превосходные результаты также и при решении других жестких задач (например, задачи CUSP из [6], имеющей 96‑й порядок и 32 жестких собственных значения). Таким образом, явные адаптивные методы могут успешно конкурировать с лучшими неявными методами при решении многих жестких задач.

Таблица 2. Точность и вычислительные затраты при решении двух жестких задач

Задача Метод scd Nf NJac NLU
VDPOL Адаптивный 4
Гира
DIRK33
DIRK44
RADAU
4.11
3.39
3.99
4.08
4.44
1756
1361
2221
2711
2214
0
151
75
94
165
0
265
278
242
231
OREGO Адаптивный 4
Гира
DIRK33
DIRK44
RADAU
3.50
1.78
2.36
3.66
3.12
2931
1410
2396
3160
3416
0
236
149
157
200
0
356
317
350
267

4.2. Локально-неустойчивые задачи

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

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

x'1 = x2,   x1(0) = 2,
x'2 = μ (1 - x12) (x1 + x2),   x2(0) = 0,   0 ≤ t ≤ 3,

построенную аналогично осциллятору Ван-дер-Поля, но имеющую более резкий переход от устойчивого состояния к неустойчивому. При больших значениях μ задача является жесткой. Правильное решение этой задачи при μ = 108 показано на рисунке синей линией. Однако неявные методы обычно дают неправильное решение, показанное на рисунке красной линией, которое практически не изменяется при изменении задаваемой точности в широких пределах. Этот факт может привести к ошибочному мнению, что данное решение является правильным. На приведенном примере видно, что при моделировании процессов, имеющих быстро нарастающий, катастрофический характер, неявные методы могут давать совершенно неверный результат, соответствующий устойчивому процессу. Явные адаптивные методы позволяют качественно правильно решать подобные задачи. При решении рассмотренной задачи все адаптивные методы давали правильное решение при умеренных требованиях к точности, при этом лучшими были методы Адаптивный 2, 3, 5.

5. Решение дифференциально-алгебраических уравнений (ДАУ)

Часто уравнения математической модели представлены не в нормальной форме Коши (1.1), а в виде системы ДАУ

(5.1а)

x' = f(x, y),   x(t0) = x0,

(5.1б)

0 = g(x, y),   y(t0) = y0.

Начальные условия x0y0 называются согласованными, если они удовлетворяют алгебраическому уравнению. Реализовать уравнения вида (5.1) в ПК «МВТУ» можно с помощью типового блока F(Y) = 0.

5.1. Метод пространства состояний

Для численного решения уравнений (5.1) можно использовать два способа [6]. Первый из них (метод пространства состояний) основан на приведении уравнений к нормальной форме (1.1) путем численного решения алгебраической подсистемы (5.1б) при заданном векторе x. Подставляя затем полученное значение вектора y в (5.1а), получим искомые значения производных. Этот способ применяется в ПК «МВТУ» для расчета начального значения вектора y, а также в процессе интегрирования явными методами. Начальное значение вектора y можно задать в диалоговом окне блока F(Y) = 0, но если оно не согласованно, то будет уточнено на этапе инициализации модели. Для решения алгебраической подсистемы можно использовать один из трех методов: простых итераций, Ньютона-Рафсона, Бройдена. Метод пространства состояний позволяет разделить задачи интегрирования ОДУ и решения алгебраических уравнений, поэтому его можно применять в сочетании с любым методом интегрирования. Но его нельзя использовать при решении задач высших индексов, когда алгебраическая подсистема вырождена.

5.2. Метод ε-вложения

Второй способ (метод ε-вложения) основан на совместном решении дифференциальной и алгебраической подсистем и может быть интерпретирован как решение сингулярно возмущенной задачи

x' = f(x, y),   x(t0) = x0,
ε y' = g(x, y),   y(t0) = y0.

при ε = 0. В ПК «МВТУ» этот способ используется в процессе интегрирования неявными методами. Применяя метод Рунге-Кутты (3.1), получим формулы одного шага интегрирования системы (5.1) в виде

(5.2а)

;

(5.2б)

,

(5.2в)

.

Формулы (5.2б), (5.2в) задают систему нелинейных алгебраических уравнений относительно векторов ki, li, i = 1, …, s (векторы стадийных значений Xi, Yi, i = 1, …, s нетрудно исключить). При реализации приведенных формул из уравнений (5.2б), (5.2в) находим векторы ki, li, которые подставляем в (5.2а). Метод ε-вложения позволяет решать задачи высших индексов, но его можно использовать только в сочетании с неявным методом интегрирования, поскольку для явных методов система алгебраических уравнений (5.2б), (5.2в) будет вырожденной.

5.3. Системы ДАУ высших индексов

Согласно определению Гира и др. (см. [6]), индекс дифференцирования системы (5.1) есть наименьшее число аналитических дифференцирований, требующихся для того, чтобы из уравнений (5.1) можно было бы получить систему ОДУ в форме Коши. При этом каждое дифференцирование понижает индекс на 1. Продифференцировав алгебраическую подсистему (5.1б) и обозначив gx = ∂g / ∂x, gy = ∂g / ∂y, получим

(5.3)

0 = gx x' + gy y'.

Если матрица gy обратима, то из (5.1), (5.3) можно получить систему в нормальной форме Коши

x' = f(x, y),   x(t0) = x0,
y' = -gy-1 gx f(x, y),   y(t0) = y0.

Итак, если матрица gy обратима в любой точки на траектории решения, то система (5.1) имеет индекс 1, а в противном случае (если матрица gy вырождена) индекс системы больше 1. Системы высших индексов (2 и выше) возникают при решении многих прикладных задач. Например, уравнения механической системы со связями, сформированные методом Лагранжа, имеют индекс 3. Такие системы наиболее трудны для численного решения и могут быть решены только неявными методами.

В качестве примера рассмотрим систему ДАУ

(5.4а)

x' = u,   u' = -x z,
y' = v,   v' = -1 - y z,

(5.4б)

0 = x2 + y2 - 1,

описывающую колебания маятника. Продифференцировав алгебраическое уравнение (5.4б), получим

(5.5)

0 = x u + y v,

а продифференцировав (5.5), получим уравнение

(5.6)

0 = -z (x2 + y2) - y + u2 + v2

из которого можно выразить алгебраическую переменную z через дифференциальные переменные x, y, u, v. Таким образом, система (5.4а), (5.6) имеет индекс 1, система (5.4а), (5.5) – индекс 2, а исходная система (5.4) – индекс 3.

Трудность решения систем ДАУ высших индексов обусловлена эффектом, который получил название «феномен снижения порядка» [6] и заключается в том, что реальный порядок метода оказывается меньше, чем полученный исходя из классических представлений. Этот эффект проявляется при решении жестких и дифференциально-алгебраических систем (особенно заметно – высших индексов). В Приложении приведены результаты решения уравнений маятника, сформированных в виде систем индексов 1, 2, 3. Видно, что снижение порядка при увеличении индекса приводит к существенному увеличению вычислительных затрат, причем только решатель DIRK44 смог решить задачу индекса 3. Для эффективного решения систем ДАУ высших индексов разрабатывают специальные методы, удовлетворяющие дополнительным условиям порядка [6]. Среди решателей ПК «МВТУ» пока нет таких методов, однако решатель DIRK44 позволяет решать задачи индексов 2 и 3. Для более эффективного решения таких задач нами разработан метод 4-го порядка с таблицей Бутчера

который предполагается реализовать в одной из последующих версий ПК «МВТУ». По сравнению с другими методами ESDIRK этот метод имеет более высокий реальный порядок при решении систем ДАУ индексов 2 и 3.

Приложение

Результаты тестирования решателей ОДУ ПК «МВТУ» и системы MATLAB

Литература и Интернет

  1. ПК «МВТУ». Учебная версия, документация, методические материалы.
  2. Козлов О.С., Кондаков Д.Е., Скворцов Л.М. и др. Программный комплекс «Моделирование в технических устройствах». Статья на сайте model.exponenta/ru.
  3. Козлов О.С., Кондаков Д.Е., Скворцов Л.М. и др. Программный комплекс для исследования динамики и проектирования технических систем // Информационные технологии. 2005. № 9. С. 20–25.
  4. Современные численные методы решения обыкновенных дифференциальных уравнений / Под ред. Дж. Холла и Дж. Уатта. М.: Мир, 1979.
  5. Хайрер Э., Нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М.: Мир, 1990.
  6. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М.: Мир, 1999.
  7. Скворцов Л. М. Адаптивные методы численного интегрирования в задачах моделирования динамических систем // Изв. РАН. Теория и системы управления. 1999. № 4. С. 72–78.
  8. Скворцов Л. М. Явные адаптивные методы численного решения жестких систем // Математическое моделирование. 2000. № 12. С. 97-107.
  9. Скворцов Л. М. Диагонально неявные FSAL-методы Рунге-Кутты для жестких и дифференциально-алгебраических систем // Математическое моделирование. 2002. Т. 14. № 2. С. 3–17.
  10. Скворцов Л. М. Точность методов Рунге Кутты при решении жестких задач // Ж. вычисл. матем. и матем. физ. 2003. Т. 43. № 9. С. 1374–1384.
  11. Test Set for Initial Value Problem Solvers. Release 2.2. August 2003.
  12. Козлов О. С., Скворцов Л. М. Тестовое сравнение решателей ОДУ системы MATLAB // Всероссийская научная конференция «Проектирование научных и инженерных приложений в среде MATLAB». М.: Изд-во ИПУ РАН, 2002. С. 53–60.

21.11.2005