Московский физико-технический институт
Опубликован: 25.10.2007 | Доступ: свободный | Студентов: 3229 / 659 | Оценка: 4.50 / 4.33 | Длительность: 24:00:00
ISBN: 978-5-9556-0065-9
Специальности: Программист, Математик

Лекция 10: Численные методы решения жестких систем обыкновенных дифференциальных уравнений

9.5. Формулы дифференцирования назад и методы Гира. Представление Нордсика

В "Численные методы решения задачи Коши для систем обыкновенных дифференциальных уравнений" вкратце был изложен альтернативный методам Рунге - Кутты подход к построению численных методов для решения ОДУ, основанный на расширении шаблона разностной схемы. При переходе от точки tn к tn + 1 использовались значения решения (или функций от него) в предыдущих точках. Полученные численные методы (первые схемы такого рода были получены Адамсом) носят название линейных многошаговых. Однако применение методов Адамса для решения ЖС ОДУ приводит к неутешительному результату. Отчасти он обоснован результатом, полученным Далквистом.

Теорема (Далквиста (или второй барьер Далквиста) [9.8], [9.9].)

Не существует A - устойчивых линейных многошаговых схем с порядком аппроксимации выше второго.

Попробуем ввести другой класс линейных многошаговых методов — это так называемые формулы дифференцирования назад (или ФДН - методы) [9.8].

Общий вид ФДН - метода таков:

{u_{n + 1} + \sum\limits_{j = 1}^{k}{\alpha_j u_{n + 1 - j}} = {\tau}\beta f_{n + 1}, } ( 9.5)

где коэффициенты \alpha_j выбираются из условий аппроксимации метода. В отличие от методов типа Рунге - Кутты, при использовании ФДН - методов нелинейная система алгебраических уравнений для определения un + 1 имеет меньшую размерность, следовательно, требуется меньшее число операций для нахождения решения.

Семейство A(\alpha) - устойчивых ФДН - методов с достаточно большим значением угла полураствора \alpha носит общее название методов Гира. Коэффициенты методов Гира, представленных в виде (9.5), приведены в таблице 9.10 (см. также книгу [9.13]).

Таблица 9.10. Коэффициенты методов Гира
k B \alpha_0 \alpha_1 \alpha_2 \alpha_3 \alpha_4 $ \alpha_5 $
1 1 - 1
2 2/3 1/3 - 4/3
3(86) 6/11 - 2/11 9/11 - 18/11
4(73, 35) 12/25 - 3/25 16/25 - 36/25 48/25
5(51, 84)
$ \frac{60}{137} $
$ \frac{- 12}{137} $
$ \frac{75}{137} $
$ \frac{- 200}{137} $
$ \frac{300}{137} $
$ \frac{ - 300}{137} $
6(17, 84)
$ \frac{600}{147} $
$ \frac{10}{147} $
$ \frac{- 72}{147} $
$ \frac{225}{147} $
$ - \frac{400}{147} $
$ \frac{450}{147} $
$ \frac{- 360}{147} $

Величина k характеризует количество точек в шаблоне и порядок аппроксимации ФДН - метода. В первом столбце таблицы также приведен угол полураствора \alpha для A(\alpha) - устойчивых методов. ФДН - методы с порядком аппроксимации 7 и выше безусловно неустойчивы.

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

Рассмотрим в качестве примера построение метода Нордсика [9.2], [9.8] для модельной задачи

\dot {u} = f(t;  u), u(0) = u_0 ( 9.6)

и введем также в рассмотрение вектор Нордсика

$  z_n = {(u_n, {\tau}u^{\prime}_n, \frac{\tau^2}{2}u^{\prime\prime}_n, 
\ldots , \frac{\tau^k}{k!}u_n^{(k)})}^T, $

в который, кроме искомого значения функции, включены приближенные значения первых k производных в узлах сетки. Для того чтобы осуществить переход к следующему значению zn + 1 необходимо задать правила вычисления компонентов вектора z по данным в текущий момент времени. Вслед за [9.2], [9.8] ограничимся рассмотрением случая k = 3. Для этого разложим (9.6) в ряд Тейлора в окрестности tn:

$  {u_{n + 1} = u_n + {\tau}u^{\prime}_n + \frac{\tau^2}{2}u^{\prime\prime}_n + \frac{\tau^3 }{3!}u_n^{(3)} + \frac{\tau^4}{4!}u^4 (\xi), }  $ ( 9.7)

для записи ряда использован остаточный член в форме Лагранжа. Кроме того, используем следующие следствия приведенного выше разложения:

$  {{\tau}u^{\prime}_{n + 1} = {\tau}u^{\prime}_n + 2\frac{\tau^2}{2}u^{\prime\prime}_n + 3\frac{\tau^3}{3!}u_n^{(3)} +  4\frac{\tau^4}{4!}u^{(4)}(\xi), }  $ ( 9.8)

$  {\frac{\tau^2}{2}u^{\prime\prime}_{n + 1} = 2\frac{\tau^2}{2}u^{\prime\prime}_n + 3\frac{\tau^3}{3!}u_n^{(3)} + 6\frac{\tau^4}
{4!}u^{(4)}(\xi), }  $ ( 9.9)

$  {\frac{\tau^3}{3!}u^{(3)}_{n + 1} = 3\frac{\tau^3}{3!}u_n^{(3)} +  4\frac{\tau^4}{4!}u^{(4)}(\xi)}  $ ( 9.10)

а конкретное значение \xi выберем из условий аппроксимации уравнения (9.6), т.е. u'n + 1 = f(tn + 1, un + 1).

Подставляя последнее равенство в (9.8), получим:

$   {4\frac{{{\tau}^4 }}{{4!}}u^{(4)} (\xi ) = {\tau}(f_{n + 1} - 
f_n^{p}), }   $ ( 9.11)

а f_n^{p} выражается через компоненты вектора Нордсика:

$  f_n^{p} = u^{\prime}_n + {\tau}^2 u^{\prime\prime}_n + 
\frac{{{\tau}^3 }}{2}u_n^{(3)} . $

Заманчивая идея подставить выражение (9.11) во все соотношения (9.7), (9.8), (9.9), (9.10) приводит к неустойчивому методу (см. [9.2], [9.8]). При этом (неустойчивый) метод имеет четвертый порядок аппроксимации. Основная идея метода Нордсика заключается в том, чтобы, несколько "испортив" порядок аппроксимации, добавляя в выражения (9.7), (9.8), (9.9), (9.10) представление (9.11), добиться максимальной устойчивости метода при приемлемом порядке аппроксимации.

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

z_{n + 1} = Pz_n + l(\tau f_{n + 1} - e_1^{T}Pz_n), ( 9.12)

где матрица P — треугольная матрица Паскаля, определяемая соотношением (для произвольного числа компонентов вектора Нордсика k ): P_{ij} = C_i^{j}, (0 \le i \le j \le k) здесь C^{j}_iбиномиальный коэффициент. В противном случае Pij = 0. Вектор E1 определяется как (0, 1, 0, 0, ..., 0), а вектор l — как (l_0, 1, l_2, ..., lk), значение 1 выбирается из условия нормировки. Относительно первого компонента вектора Норсдика un + 1 система уравнений нелинейна, все остальные переменные входят в систему линейным образом.

Для рассматриваемого случая k = 3 Норсдик получил набор коэффициентов l = (3/8, 1, 3/4, 1/6) из условий минимума погрешности и обращения в нуль собственных чисел матрицы M= P - le_1^{T}P. Можно получать набор коэффициентов для метода Гира в представлении Норсдика, например, из условий жесткой устойчивости (например, ослабленного требования L - устойчивости, когда требуется лишь |R(z)|  \to 0 при {\mathop{\rm Re}\ \nolimits}  {\tau}\lambda  \to - \infty ). В этом случае для рассмотренного выше примера l = (6/11, 1, 6/11, 1/11).

Покажем, что последний из приведенных методов Нордсика эквивалентен методу ФДН Гира при k = 3. Для этого просто для трех последовательных шагов по времени исключаем из формул типа (9.7), (9.8), (9.9), (9.10) (точнее, из (9.12)) значения производных решения. После нетрудных преобразований получим формулу Гира. Другой рассмотренный вариант метода Нордсика приведет к одной из неявных формул Адамса. Подробнее в [9.8].

Отметим, что метод Гира в представлении Нордсика оказывается самостоятельно стартующим. При старте можно положить вектор Нордсика для данной задачи равным, например, z0 = (u0, 0, 0, ..., 0), что позволит начать вычисления, но приведет к уменьшению порядка аппроксимации. Таким образом, многозначный (по введенной в [9.2] терминологии) вариант метода Гира обладает переменным порядком аппроксимации: стартуя как метод первого порядка, по завершении разгонного участка метод стремится к максимально возможному для данной формулы порядку.

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

Евгений Бурцев
Евгений Бурцев
Россия
Максим Виноградов
Максим Виноградов
Россия, Москва