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

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

  1. Пограничный слой A'A. На этом участке за малое время O(\varepsilon ) траектория из точки {u0, v0} переходит в \varepsilon - окрестность кривой f = 0. Здесь траектория почти горизонтальна и приближенно определяется дифференциальным уравнением

    \begin{gather*}
\dot {u} = \frac{1}{\varepsilon } f(u, v), \\  
v(t)  \approx  v_0, \\  
u(0) = u_0 .
\end{gather*}

    ( f/\varepsilon  \gg g, так как f, g  \sim  O(1), 0 < \varepsilon \ll 1 ).

    В окрестности кривой f(u, v) = 0 имеем f'u < 0, поэтому допустима оценка

    $  \frac{d f}{d t} = f^{\prime}_u  \cdot {\dot {u}} = \frac{1}{\varepsilon } f \cdot {f}^{\prime}_u , $

    откуда видно, что f(u, v0) стремится к нулю как экспонента с показателем

    $ \frac{1}{\varepsilon } f^{\prime}_u   $,
    т.е. f становится малой величиной ( \approx  0 ) за время t_{A'A} = O(\varepsilon ).

  2. Квазистационарный режим. Движение точки {u(t), v(t)} продолжается уже по участку кривой AB, f(u, v) = 0 и описывается системой

    \begin{gather*}
f(u, v) = 0, \\  
\dot {v} = g(u, v).
\end{gather*}

    На этом участке за время t_{AB} \sim  O(1), что видно из второго уравнения, точка подвигается от A к B, пока система ОДУ устойчива. В случае невозмущенной системы точка могла бы и далее продвигаться по участку BD, но для полной системы эта ветвь оказывается неустойчивой (f'u > 0) и траектория "срывается" на устойчивую ветвь CD в точке B, в которой f'u = 0.

  3. Пограничный слой. На участке BC точка {u(t), v(t)} "перескакивает" из B в C за малое время O(\varepsilon ) ; движение здесь, как и на ветви AB, приближенно описывается уравнениями

    \begin{gather*}
\dot {u} = \frac{1}{\varepsilon } f(u, v), \\  
v(t)  \approx  const .
\end{gather*}

  4. Квазистационарный режим. Движение по ветви CD, как и по ветви AB, описывается уравнениями

    \begin{gather*}
f(u, v) = 0, \\  
\dot {v} = g(u, v), 
\end{gather*}

    и длится t_{CD} \sim  O(1).

  5. Пограничный слой. На неустойчивой ветви DA также, как и на ветви BC, происходит скачок из точки D за время t_{DA} \sim  O(\varepsilon ) в устойчивую точку A, и т.д.

    Такое поведение траектории (замкнутая кривая) называется предельным циклом. Для жестких систем периодические решения называют иногда релаксационными колебаниями [9.6], [9.7].

    Таким образом, это характерно для жестких систем, траектория состоит из чередующихся участков быстрого (за время t \sim  O(\varepsilon ) ) и медленного (за время t \sim  O(1) ) изменения решения.

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

$  {\tau} \cdot \frac{1}{\varepsilon } |{f^{\prime}_u}| \ll 1.  $

В зоне квазистационарного режима часто оказывается, что интегрирование с таким шагом слишком дорого. Можно, правда, разрешить уравнение f(u, v) = 0 (u = \varphi (v)) относительно u и далее интегрировать его, предварительно реализовав алгоритм перехода на другой шаг:

$ \dot {v} = g(v, \varphi (v)),  $

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

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

\begin{gather*}
\frac{{u_{k + 1} - u_k}}{{\tau}} = \frac{1}{\varepsilon } f(u_{k + 1}, v_{k + 1}), \\  
\frac{{v_{k + 1} - v_k}}{{\tau}} = g(u_{k + 1}, v_{k + 1}).
  \end{gather*}

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

Первый пограничный слой будет пройден за один шаг:

\begin{gather*}
u_1 = u_0 + \frac{{\tau}}{\varepsilon } f(u_1, v_1 ), \\  
v_1 = v_0 + {\tau}g(u_1, v_1 ), 
 \end{gather*}

так как \varepsilon \ll {\tau}, {\tau}g мало.

Далее следует процесс численного интегрирования на устойчивой ветви AB. В окрестности точки B поведение численного решения по неявной схеме осложняется. Это связано с тем, что рассматриваемая система в данной окрестности может иметь более одного решения. При этом, по крайней мере, одно из решений нелинейной алгебраической системы может лежать на неустойчивой ветви CD кривой f(u, v) = 0. Возможно, что при выборе большего шага интегрирования получится именно это нефизическое решение, к которому сойдутся итерации.

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

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