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

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

9.4. Численные методы решения ЖС ОДУ. Семейства неявных методов Рунге - Кутты и Розенброка

Рассмотрим другие методы для численного решения как линейных жестких систем ОДУ вида (9.1), так и нелинейных систем общего вида

$ \dot {u} = {F}(t;{u}), {u}(0) = u_0 $

и автономных ЖС ОДУ

$ \dot {u} = {F}(u), \quad {u}(0) = u_0, $

см. также [9.1], [9.4], [9.9], [9.13], [9.14], [9.15], [9.16]. Из соображений устойчивости метода предпочтение естественно отдать неявным методам. Простейшими из них являются:

  1. неявный метод Эйлера (приведем его вид для случая автономной ЖС ОДУ, для неавтономной системы формулы очевидны):

    $  \frac{u_{n + 1} - u_n}{{\tau}} = f(u_{n + 1}), u_0 = a; 
  $

  2. метод трапеций

    $  \frac{u_{n + 1} - u_n}{{\tau}} = \frac{f(u_{n + 1}, t_{n + 1}) + f(u_n, t_n)}{2}, u_0 = a;

  3. метод прямоугольников (правило средней точки)

    $  \frac{u_{n + 1} - u_n}{{\tau}} = f(u_{n + 1/2}, t_{n + 1/2}), 
u_0 = a, $

    где

    $  t_{n + 1/2} = \frac{t_{n + 1} + t_n}{2}, u_{n + 1/2} = \frac{u_{n + 1} + u_n}{2}.  $

Среди одношаговых методов для решения жестких систем наиболее известны методы Рунге-Кутты. Все приведенные в данном параграфе формулы можно рассматривать как частные случаи неявных одно - и двухстадийных методов из этого семейства. Не останавливаясь на получении приводимых коэффициентов, выпишем наиболее известные из формул Рунге - Кутты, используя таблицу Бутчера. О представлении методов Рунге - Кутты в виде таблиц Бутчера — в "Численные методы решения задачи Коши для систем обыкновенных дифференциальных уравнений" .

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

  1. Методы Гаусса соответственно 2 - го, 4 - го, и 6 - го порядков представлены в табл. 9.1, 9.2, 9.3. Основаны эти методы на соответствующих квадратурных формулах Гаусса, рассмотренных в "Численное интегрирование" . Первый метод совпадает с правилом средней точки. Второй метод (4 - го порядка) носит название метода Хаммера - Холлинсворта.
    Таблица 9.1.
    1/2 1/2
    1
    Таблица 9.2.
    $ \frac{1}{2} - \frac{\sqrt{3}}{6} $
    $ \frac{1}{4} $
    $ \frac{1}{4} - \frac{\sqrt{3}}{6} $
    $ \frac{1}{2} + \frac{\sqrt{3}}{6} $
    $ \frac{1}{4} + \frac{\sqrt{3}}{6} $
    $ \frac{1}{4} $
    1/2 1/2
    Таблица 9.3.
    $ \frac{1}{2} - \frac{\sqrt{15}}{10} $
    $ \frac{5}{36} $
    $ \frac{2}{9} - \frac{\sqrt{15}}{15} $
    $ \frac{5}{36} - \frac{\sqrt{15}}{30} $
    $ \frac{1}{2} $
    $ \frac{5}{36} + \frac{\sqrt{15}}{24} $
    $ \frac{2}{9} $
    $ \frac{5}{36} - \frac{\sqrt{15}}{24} $
    $ \frac{1}{2} + \frac{\sqrt{15}}{10} $
    $ \frac{5}{36} + \frac{\sqrt{15}}{30} $
    $ \frac{2}{9} + \frac{\sqrt{15}}{15} $
    $ \frac{5}{36} $
    $ \frac{5}{18} $
    $ \frac{4}{9} $
    $ \frac{5}{18} $
  2. Методы Радо IIA порядков 1, 3, 5 представлены соответственно в табл. 9.4, 9.5, 9.6. Основаны на квадратурных формулах Радо, принадлежащих к семейству формул Гаусса. Метод первого порядка является неявным методом Эйлера.
    Таблица 9.4.
    1 1
    1
    Таблица 9.5.
    1/3 5/12 - 1/12
    1 3/4 1/4
    3/4 1/4
    Таблица 9.6.
    $ \frac{4 - \sqrt{6}}{10} $
    88 - 7\sqrt{6} $
    $ \frac{296 - 169\sqrt{6}}{1 800} $
    $ \frac{- 2 + 3\sqrt{6}}{225} $
    $ \frac{4 + \sqrt{6}}{10} $
    $ \frac{296 + 169\sqrt{6}}{1 800} $
    $ \frac{88 + 7\sqrt{6}}{360} $
    $ \frac{- 2 - 3\sqrt{6}}{225} $
    1
    $ \frac{16 - \sqrt{6}}{36} $
    $ \frac{16 + \sqrt{6}}{36} $
    $ \frac{1}{9} $
    $ \frac{16 - \sqrt{6}}{36} $
    $ \frac{16 + \sqrt{6}}{36} $
    $ \frac{1}{9} $
  3. Методы Лобатто IIIА 2-го, 4-го и 6-го порядков точности см. в табл. 9.7, 9.8, 9.9 и 9.10 соответственно. Видно, что метод второго порядка точности является неявным методом трапеций.
    Таблица 9.7.
    0 0 0
    1 1/2 1/2
    1/2 1/2
    Таблица 9.8.
    0 0 0 0
    1/2 5/24 1/3 - 1/24
    1 1/6 2/3 1/6
    1/6 2/3 1/6
    Таблица 9.9.
    0 0 0 0 0
    $ \frac {5 - \sqrt{5}}{10} $
    $ \frac{11 + \sqrt{5}}{120} $
    $ \frac{25 - \sqrt{5}}{120} $
    $ \frac{25 - 13\sqrt{5}}{120} $
    $ \frac{- 1 + \sqrt{5}}{120} $
    $ \frac {5 + \sqrt{5}}{10} $
    $ \frac{11 - \sqrt{5}}{120} $
    $ \frac{25 + 13\sqrt{5}}{120} $
    $ \frac{25 + \sqrt{5}}{120} $
    $ \frac{- 1 - \sqrt{5}}{120} $
    1 1/12 5/12 5/12 1/12
    1/12 5/12 5/12 1/12

    Рассмотрим теперь, как строится функция устойчивости для методов Рунге - Кутты [9.4]. Уже отмечалось выше, что при построении функции устойчивости рассматривается модельное уравнение вида

    $ \dot {v} = \lambda v, v(0) = v_0. $

    Запишем теперь метод Рунге - Кутты для решения приведенного выше уравнения в общей форме с использованием новых переменных Y:

    \begin{gather*}
v_{{n} + 1} = v_n + {\tau}\sum {\gamma_p f(t_n + \alpha_p {\tau}, Y_p) \equiv v_n + {\tau}\lambda \sum {\gamma_p Y_p} }, \\  
Y_p = v_n + {\tau}\sum\limits_{k = 1}^{s}{\beta_{{kp}} f(t_n + \gamma_k {\tau}, Y_k) \equiv v_n + {\tau}\lambda \sum {\beta_{{kp}} Y_k} }.
\end{gather*}

    Приведенные выше формулы можно рассматривать как систему линейных уравнений относительно новых переменных Y1, ..., Yk, vn + 1 следующего вида:

    \left( \begin{array}{ccccc}
 {1 - z\beta_{11}} & {- z\beta_{12}} & \ldots & {- z\beta_{1s}} & 0  \\ 
 {- z\beta_{21}} & {1 - z\beta_{22}} & \ldots & {- z\beta_{2s}} & 0  \\ 
  \ldots & \ldots & \ldots & \ldots & \ldots   \\ 
  \ldots & \ldots & \ldots & \ldots & \ldots   \\ 
 {- z\beta_{s1}} & {- z\beta_{s2}} & \ldots & {1 - z\beta_{ss}} & 0 \\ 
 {- z\gamma_1 } & {- z\gamma_2 } & \ldots & {- z\gamma_{s}} & 1  \\ 
\end{array} \right) \left( \begin{array}{l}
   {Y_1 }  \\ 
   {Y_2 }  \\ 
    \ldots   \\ 
    \ldots   \\ 
   {Y_s}  \\ 
   {v_{n + 1}}  \\ 
\end{array} \right) = \left( \begin{array}{l}
   {v_n}  \\ 
   {v_n}  \\ 
    \ldots   \\ 
    \ldots   \\ 
   {v_n}  \\ 
   {v_n}  \\ 
\end{array} \right) .

    Необходимо выразить vn + 1 через vn. Для этого воспользуемся правилом Крамера. Ответ можно записать в следующей форме:

    $  x_{n + 1} = \frac{\det (E - zB + zea^T)}{\det (E - zB)}x_n,   $

    где E — единичная матрица размера S x S, B - матрица коэффициентов \beta_{ij}, входящая в таблицу Бутчера, E — единичный вектор размерности S (столбец), aT — строка коэффициентов \gamma _{1}, входящая в таблицу Бутчера, S — число стадий метода Рунге - Кутты.

    Условием устойчивости метода будет

    $  \left|{\frac{\det (E - z B + z ea^T)}{\det (E - zB)}}\right| < 1. $

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

$  (E - a {\tau}B  - b {\tau}^2 B^2 ) \frac{u_{n + 1} - u_n}{{\tau}} = 
f[u_n + c{\tau}f(u_n)] . $

Здесь

$  B = \frac{\partial f}{\partial u}(u_n)  $
матрица, постоянная на данном шаге по времени. Параметры a, b, c подбираются таким образом, чтобы обеспечить максимально возможный порядок точности.

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

a = 1, 077;   b = - 0, 372;   c = - 0, 577 .

Такую схему иногда называют методом с одной итерацией, имея в виду, что вычисление обратной матрицы сравнимо по количеству арифметических операций с одной итерацией метода Ньютона. Преимущество методов типа Розенброка перед прочими классами численных методов ЖС ОДУ заключается в том, что для определения решения на верхнем временном слое необходимо решать уже линейную систему алгебраических уравнений.

Рассмотрим связь между методами Розенброка и Рунге - Кутты.

Определение. ([9.9]) S - стадийный метод Розенброка для решения автономной системы ЖС ОДУ имеет вид

\begin{gather*}
k_i = {\tau}f \left({u_n + \sum\limits_{j = 1}^{i - 1}{\beta_{i, j}k_j}}\right) + {\tau}B \sum\limits_{j = 1}^{i}{\mu_{i, j}k_j}, \\  
i = 1, \ldots , S, \\  
u_{n + 1} = u_n + \sum\limits_{k = 1}^{S}{\gamma_k k_k}, 
\end{gather*}

где \beta, \gamma, \muуправляющие коэффициенты метода. Как и выше, матрица Якоби правой части системы B вычисляется по данным в точке tn. Связь последних формул с выражениями для явных методов Рунге - Кутты очевидна.

Несколько сложнее представляются методы Розенброка в случае неавтономной ЖС ОДУ [9.9]:

\begin{gather*}
k_i = {\tau}f \left({t_n + \alpha_i {\tau}, u_n + \sum\limits_{j = 1}^{i - 1}{\beta_{i, j}k_j} }\right) + {\tau}B\sum\limits_{j = 1}^{i}{\mu_{i, j}k_j} + \mu_i \tau^2 B, \\  
i = 1, \ldots , S, \\  
u_{n + 1} = u_n + \sum\limits_{k = 1}^{S}{\gamma_k k_k}, \\  
\alpha_i = \sum\limits_j \beta_{i, j}, \mu_i = \sum\limits_j{\mu_{i, j}}. 
\end{gather*}

Вывод условий порядка методов Розенброка — достаточно громоздкий, и всем заинтересованным читателям можно рекомендовать разобраться в выкладках в книге [9.9]. Отметим, что широкое распространение получают в последнее время вложенные методы Розенброка высокого порядка аппроксимации, имеющие очень хорошие вычислительные качества. Возможно, что новые методы типа Розенброка способны будут вытеснить из вычислительной практики наиболее распространенные до этого в вычислительной практике методы Гира.