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

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

8.6. Задачи

  1. Семейство явных методов Рунге - Кутты

    Пусть дана следующая система обыкновенных дифференциальных уравнений:

    \begin{gather*} \dot {y} = f(t;y), \\ 
y(0) = y_0 .
\end{gather*}

    Для численного решения задачи проведем следующие вычисления:

    \begin{gather*}
k_1 = f(t_n, y_n), \\ 
k_2 = f(t_n + c_2 {\tau}, y_n + {\tau}b_{21}k_1 ), \\ 
k_{s} = f(t_n + c_{s} {\tau}, y_n + {\tau}\sum\limits_{r = 1}^{s - 1}{b_{sr}k_r}), \\ 
y_{n + 1} = y_n + {\tau}\sum\limits_{p = 1}^s{a_p k_p}, 
\end{gather*}

    где a, b, C — действительные числа (вектор yn считается известным, т.к. y0 — начальные условия, по этим значениям можно найти y1 и т.д.).

    Численный метод называется явным методом Рунге - Кутты порядка S или S - стадийным.

    Будем также записывать метод Рунге - Кутты в виде таблицы Бутчера — таблицы коэффициентов метода:

    0 0 0 \ldots 0 0
    C2 b21 0 \ldots 0 0
    \ldots \ldots \ldots \ldots \ldots \ldots
    cs bs1 bs2 \ldots bss - 1 0
    a1 a2 \ldots as - 1 as

    Для метода порядка S вывести условия, связывающие коэффициенты таблицы Бутчера, необходимые для обеспечения p - го порядка аппроксимации метода ( p \le s ).

    Указание 1. Разложить проекцию на сетку точного решения задачи Коши в ряд Тейлора в окрестности точки tn. Использовать следствия исходного дифференциального уравнения, например:

    $ 
\frac{d^2 y}{d t^2 } = \frac{d}{dt}f(t, y) = \frac{\partial f}{\partial t} + 
\frac{\partial f}{\partial y} \frac{\partial y}{\partial t} = \frac{\partial f}{\partial t} + \frac{\partial f}{\partial y}f. $

    Указание 2. Для простоты все вычисления провести для случая одного нелинейного уравнения. Требуемые условия приведены в [8.3, с. 153, 162]. В качестве упражнения найти все явные методы Рунге - Кутты с S = p = 3, S = p = 4.

  2. Семейство неявных методов Рунге - Кутты

    Метод вида

    \begin{gather*}
k_1 = f(t_n + c_1 {\tau}, y_n + {\tau}\sum\limits_{l = 1}^{s}{b_{1l}k_l} ), \\ 
\ldots \\ 
k_{s} = f(t_n + c_{s} {\tau}, y_n + {\tau}\sum\limits_{l = 1}^{s}{b_{sl}k_l}), \\ 
y_{n + 1} = y_n + {\tau}\sum\limits_{l = 1}^{s}{a_l k_l}, 
\end{gather*}

    где k1, ..., ks определяются как решение нелинейной системы уравнений, называется неявным методом Рунге - Кутты порядка S ( S - стадийным ). Как будет выглядеть для него таблица Бутчера?

    Вывести условия аппроксимации порядка p на решении ( p = 1, 2, 3, 4 ; S = 2 ).

    Следует обратить внимание на то, что для определения k1, k2, , ks необходимо решать систему нелинейных уравнений. Какова ее размерность?

    В чем состоит особенность методов с bij = 0 при j > i?

    (Это так называемые полуявные или диагонально - неявные методы).

  3. Управление длиной шага
    • a) Система ОДУ решается с помощью метода Рунге - Кутты порядка аппроксимации p > 1. Описать алгоритм выбора шагов интегрирования, таких, чтобы достигнутая погрешность на каждом из них не превосходила заданную величину \varepsilon.

      Решение. Пусть переход от значения yn к yn + 1 осуществляется с шагом h. Тогда в результате работы метода получим

      {y^{n + 1} = \tilde {y}^{n + 1} + Ch^{p + 1} + O(h^{{p} + 2}), } ( 8.14)

      где $ \tilde {y}^{n + 1} $ — точное решение задачи, взятое в точке n + 1, Chp + 1 — главный член погрешности, постоянная C зависит как от коэффициентов конкретного метода, так и от решения задачи. (Из решения задачи 8.1 следует, что разложения yn + 1 и \tilde {y}^{n + 1} в ряд Тейлора совпадают до членов с номерами p + 1 включительно).

      Осуществим теперь переход от yn к yn + 1 за два этапа, используя тот же метод Рунге - Кутты с шагом h/2. Тогда имеем

      $ \tilde {\tilde {y}}^{n + 1} = \tilde {y}^{n + 1} + 2C{(h/2)}^{p + 1} + O(h^{p + 2}).} $ ( 8.15)

      Так как C оценивается через максимальное по абсолютной величине значение (p + 1) - й производной на отрезке [tn, tn + 1], можно считать, что в последнем выражении константа совпадает с константой в предыдущем равенстве. Используя два равенства, возможно либо явно вычислить Chp и исключить его из (8.15), повысив точность аппроксимации на порядок, либо поступить следующим образом. Вычитая из (8.15) равенство (8.14), имеем

      y^{n + 1} - \tilde {\tilde {y}}^{n + 1} = Ch^{p + 1} (1 - 2^{- p}) + O(h^{p + 2}),

      отсюда главный член погрешности (8.15) будет равен

      $ 
{Ch}^{p + 1} = \frac{{y^{n + 1} - \tilde {\tilde {y}}^{n + 1}}}{{1 - 2^{- p}}}.} $ ( 8.16)

      Необходимо выбрать шаг интегрирования hnew таким образом, чтобы главный член погрешности не превосходил заданное значение (Ch}_{new}^{p + 1} < \varepsilon ):

      $ \left({\frac{{h_{new}}}{h}}\right)^{p + 1} \frac{{y^{{n} + 
1} - \tilde {\tilde {y}}^{n + 1}}}{{1 - 2^{- {p}}}} < \varepsilon , $

      откуда имеем

      $ {h_{new} = \left({\frac{{y^{n + 1} - \tilde {\tilde 
{y}}^{n + 1}}}{{1 - 2^{- {p}}}}}\right)^{1/(p + 1)}h. } $ ( 8.17)

      Если же на текущем шаге погрешность (8.16) превышает \varepsilon, то шаг считается отклоненным и расчет выполняется снова со значением, найденным по формуле (8.17).

    • б) Система ОДУ решается с помощью явных методов p и p + 1 порядка аппроксимации. На каждом шаге погрешность расчетов не должна превышать \varepsilon. Как выбрать длину шага интегрирования?

      Рассуждаем аналогично тому, как это сделано в пункте 1. При расчете методом p порядка аппроксимации имеем

      {y_p^{n + 1} = \tilde {y}^{n + 1} + Ch}^{p + 1} + K_1 h^{{p} + 2} + O(h^{p + 3}), } ( 8.18)

      а p + 1 порядок дает

      {y_{p + 1}^{n + 1} = \tilde {y}^{n + 1} + K_2 h^{{p} + 2} + O(h^{{p} + 3}).} ( 8.19)

      Отсюда сразу получаем {Ch}^{p + 1} \sim y_p^{n + 1} - y_{p + 1}^{n + 1} . Эти величины сравниваются с заданной точностью.

      Существует семейство методов Рунге - Кутты, таких, что члены погрешности для результата старшего ( p + 1 ) порядка минимизируются, а вычисление погрешности p - го порядка используется лишь для управления длиной шага. При этом для порядка p и p + 1 коэффициенты b, c в таблице Бутчера (2.4.1) совпадают, а различаются лишь коэффициенты a. Кроме того, bsi = ai(yp). Такие методы построены Дорманом и Принсом и носят их имя.

  4. Рассмотрим уравнение Ферхюльста (логистическое уравнение) , описывающее динамику численности популяции:

    \begin{gather*}
{y^{\prime} = ay(1 - y), } \\ 
y(0) = y_0; 0 < y_0 < 1, a > 0
\end{gather*} ( 8.20)

    • а) Решить уравнение (8.20) точно, учитывая, что это уравнение в разделяющихся переменных. Получившаяся кривая — график точного решения — носит название логистической кривой.
    • б) Рассмотреть для (8.20) явный метод Эйлера:

      {y_{n + 1} = y_n + {\tau}ay_n \left({1 - y_n}\right).} ( 8.21)

      При каких шагах \tau метод является устойчивым?

      Решение. При t \longrightarrow \infty точное решение (8.20) асимптотически стремится к единице. Очевидно, что последовательность y_n \equiv 1 является и решением (8.20). Непосредственно проверяется, что (8.21) аппроксимирует (8.20) с точностью O(\tau ). В случае, если метод (8.21) устойчив , то решение (8.21) сходится к решению (8.20) и должно выполняться неравенство

      | yn + 1 - 1 | < | yn - 1 | .

      Если рассмотреть (8.21) как запись метода простых итераций для (8.20), то он должен сходиться к корню y = 1 при любом начальном приближении y0, лежащем между 0 и 1.

      Перепишем (8.21) в виде

      $ 
{z_{n + 1} = bz_n (1 - z_n), z_n = \frac{{{\tau}a}}{{{\tau}a + 1}}y_n; b= 1 + {\tau}a} $ ( 8.22)

      и для (8.22) воспользуемся результатами, полученными в "Численное решение нелинейных алгебраических уравнений и систем" . Таким образом, условия устойчивости для (8.21) будут следующими:

      b = 1 + \tau  a < 3, откуда \ \tau  < 2/a.
    • с) Описать сценарий развития вычислительной неустойчивости для задачи (8.21), можно воспользоваться материалом "Численное решение нелинейных алгебраических уравнений и систем" . Следует обратить внимание на то, что когда метод неустойчив, в отличие от линейных задач, модуль разности численного и точного решений остается ограниченным при любом сколь угодно большом t при {2/a} \le {\tau}< {3/a}.