§ 9. Особенности интегрирования систем уравнений

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

Формальное отличие состоит в том, что в соответствующих соотношениях вместо скалярных величин участвуют некоторые матрицы или тензоры.

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

В случае использования конечно-разностной аппроксимации (5.2) соответствующая система конечно-разностных уравнений имеет вид

Для простоты предположим, что жорданова форма матрицы простая: — диагональная матрица с диагональными элементами . Положим и умножим систему (3) слева на . Получим

или

Эта система распадается на систему скалярных конечно-разностных уравнений относительно компонент векторов :

Соотношение (4) совпадает с конечно-разностной аппроксимацией для уравнения

Если в (2) перейти к новой неизвестной вектор-функции и умножить (2) слева на матрицу , то получится система скалярных уравнений (5), . В соответствии с определением вектор-функций имеем равенство

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

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

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

Шаг интегрирования должен быть существенно меньше характерного размера изменения решения, т.е. . Отсюда следует оценка снизу для числа шагов интегрирования

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

В случае, когда

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

Конечно, возникает вопрос, о каких проблемах идет речь, поскольку решение системы выписывается в явном виде? Дело в том, что мы говорим об этой задаче как о модельной; реально же метод применяется для решения какой-то, как правило, сложной задачи, и мы смотрим, как ведет себя метод в применении к простейшей задаче, где все выписывается в явном виде.

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

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

В качестве модели таких систем берется система уравнений

удовлетворяющая определенным условиям на собственные значения матрицы А. Не существует установившегося определения жестких систем. Обычно систему (6) относят к классу жестких, если величина не является большим положительным числом, а величина и

а) величина не является большим положительным числом, или

б) при умеренных значениях и .

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

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

Из проводившихся выше рассуждений для случая системы (6) видно, что численное решение задачи Коши для таких систем требует разработки специальных методов решения. Такие методы в настоящее время разработаны, и на их основе созданы соответствующие комплексы стандартных программ.

Рассмотрим простейшие варианты наиболее распространенных методов решения жестких систем.

1. Пусть уже найдено приближение к значению и ищется приближение к значению . Разложим правую часть в ряд Тейлора в точке

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

при начальном условии . Произведем следующую замену переменных и введем обозначения

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

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

Матрица является решением системы

при начальном условии .

В частном случае имеем , где матрица имеет вид

Значение можно было бы попытаться вычислить, пользуясь разложением в ряд Тейлора

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

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

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

Матрица удовлетворяет соотношению

где Е — единичная матрица. Поэтому часто бывает целесообразно пойти по следующему пути. Выбираем s такое, что ; основываясь на (9), вычисляем

а затем с помощью рекуррентной формулы (10).

В случае линейной системы алгоритм решения задачи может быть несколько упрощен. Здесь в описанном выше алгоритме на каждом шаге возникает необходимость решения системы при начальном условии . Тогда

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

поэтому зададимся некоторым s и вычисляем

а затем

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

и затем , пользуясь рекуррентной формулой

Далее находим .

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

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

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

Выражение оставим без изменения. Получим конечно-разностную аппроксимацию

Рассмотрим случай модельного уравнения , когда (12) превращается в конечно-разностное уравнение

Решение этого уравнения выписывается через корни характеристического уравнения

При корни этого уравнения удовлетворяют условию в области значений при — условию в области значений , где .

В нелинейном случае значение требуется находить из нелинейной системы (12).

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

(13)

Могло бы показаться разумным найти начальное приближение к с помощью явной формулы Эйлера

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

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

где — номер итерации.

В одном из методов решения жестких систем за принимается значение , получаемое из (14) при . Тогда имеем

Рассмотрим случай скалярного уравнения ; тогда

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

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

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

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

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

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

Получаемая в итоге аппроксимация имеет уже второй порядок.

Задача 1. На примере уравнения и получить явную формулу, выражающую через . При каких М выполняется неравенство ?

4. Как уже отмечалось ранее, явные методы численного интегрирования для жестких систем обыкновенных дифференциальных уравнений неприемлемы вследствие ограничений на шаг интегрирования. Действительно, рассмотрим модельное уравнение

и применим для его решения метод Эйлера:

Решение задачи (15) при имеет вид монотонно убывающей экспоненты , в то время как решение (16) имеет вид . Таким образом, при решение задачи методом Эйлера экспоненциально возрастает и не имеет ничего общего с реальным решением. Более того, если даже , но , то решение (15) экспоненциально убывает, но при переходе от шага к шагу меняет знак, т.е. метод Эйлера в этом случае также неприменим.

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

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

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

Изложим основные идеи метода на примере модельной задачи. В качестве такой задачи рассмотрим систему обыкновенных дифференциальных уравнений с симметричной неотрицательной матрицей А:

Будем предполагать, что спектр матрицы А может быть разделен на две части: гладкую — и жесткую — . Решение задачи (18) имеет вид

Решение задачи (18) будем искать на отрезке , используя метод Эйлера

откуда

Разлагая по собственным векторам матрицы А, получим

Тогда из (20) имеем

Вследствие (19) норма решения дифференциальной задачи не возрастает. Если потребовать, чтобы норма решения дискретной задачи также не возрастала, то из (21) можно получить, что шаг интегрирования h должен удовлетворять условию

откуда получаем оценку для :

Условие (22) является необходимым для устойчивости метода Эйлера (20). Даже в случае, когда мы находимся в области, где составляющая решения, соответствующая жесткой части спектра, близка к нулю, мы вынуждены выбирать шаг интегрирования, удовлетворяющий (22); в противном случае происходит экспоненциальное накопление вычислительной погрешности.

Предположим, что мы задали число шагов интегрирования, равное N. Тогда за N шагов, используя метод Эйлера (20) и условие устойчивости (22), мы можем получить приближенное решение до момента времени . Поставим следующую задачу. Используя переменный шаг метода Эйлера, получить за N шагов приближенное решение на как можно большем промежутке. Естественно, что при этом используемый метод должен быть устойчивым.

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

здесь — некоторые постоянные, которые наряду с надо определить.

Рассмотрение метода (23) будем проводить на примере задачи (18). Для случая задачи (18) метод (23) имеет вид

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

Таким образом, мы можем записать соотношение, связывающее в виде

где

Заметим, что

Из постановки задачи видно, что коэффициент при в многочлене как раз и является величиной, которую мы должны минимизировать. Кроме этого должно выполняться условие устойчивости , что эквивалентно выполнению неравенства при .

Таким образом, мы пришли к следующей постановке задачи. Введем класс Р многочленов степени со свободным членом 1, не превосходящих по модулю 1 на отрезке На этом классе требуется найти многочлен, производная которого в нуле является минимальной, т.е. требуется найти такой, что

Справедлива следующая

Лемма.

здесь — многочлен Чебышева степени .

Доказательство. Из определения многочленов Чебышева легко видеть, что . Предположим, что утверждение леммы неверно. Тогда существует многочлен такой, что . Рассмотрим разность . Эта разность является многочленом степени и в точке 0 обращается в нуль, так как . Пусть — точки экстремума , при этом . Так как многочлен четной степени, то . Отсюда следует, что удовлетворяет условию

Найдем количество корней многочлена на отрезке [0, М]. На отрезке всегда имеется два корня. Действительно, и, по предположению леммы, . Поскольку , то это и означает, что на имеется не менее двух корней.

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

Следует отметить, что исходная задача до сих пор не решена, так как мы лишь нашли функцию из класса Р, удовлетворяющую условию экстремума. Однако ниоткуда пока не следует, что эта функция (многочлен) представима в виде (26). Покажем, что это действительно так. Многочлен имеет на отрезке [0, М] ровно корней: , которые расположены симметрично относительно точки отрезка [0, М]. Вводя величины такие, что , искомый многочлен можно записать в виде

Приравнивая при фиксированном j выражение, стоящее под знаком произведения в (29), к выражению под знаком произведения в (26) при том же j, получаем соотношение для определения параметров :

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

Тогда

Так как

то из (32) имеем

Таким образом, за шагов численного интегрирования жесткой системы уравнений с помощью метода Эйлера мы можем получить приближенное решение на отрезке , в то время как при интегрировании той же системы уравнений с помощью явного метода с переменным шагом можно получить приближенное решение на отрезке , за то же количество шагов. Отсюда следует, что использование переменного шага интегрирования позволяет увеличить при тех же затратах процессорного времени длину отрезка интегрирования в раз, т. е. метод будет особо эффективным при больших N. В настоящее время расчеты по указанному методу проводятся до значений достигающих порядка . При этом N выбирается в каждой расчетной точке в зависимости от гладкости решения. Как правило, N выбирают среди чисел N — . Следует отметить тот факт, что при практической реализации рассмотренного выше метода, как и в случае Чебышевского ускорения итерационных методов, очень важной является проблема правильного упорядочивания параметров процесса. В настоящее время эти задачи решены и на их основе создан комплекс программ Лебедева «DUMKA» для решения жестких систем обыкновенных дифференциальных уравнений, который показал свою высокую эффективность на многочисленных задачах этого класса. Особенно эффективен данный метод по сравнению с другими при использовании многопроцессорной вычислительной техники, так как он легко распараллеливается.