Рассмотрим в качестве примера две задачи Коши для систем обыкновенных дифференциальных уравнений (ОДУ) [9.1], [9.2]:
$$$ \dot {u} = au + \frac{1}{\varepsilon } v, \dot {v} = - \frac{1}{\varepsilon } v, $$$с начальными данными u(0) = u0, v(0) = v0 ; здесь $$a \sim O(1), \varepsilon \ll 1$$ ; и линейную систему с постоянными коэффициентами
Решением первой задачи Коши являются функции
$$\begin{gather*} u(t) = u_0 e^{at} + \frac{{v_0 }}{{1 + a\varepsilon }} (e^{at} - e^{- t/\varepsilon }), \\ v(t) = v_0 e^{- t/\varepsilon }, \end{gather*}$$а второй -
u(t) = 4e- t - 3e- 1000t, v(t) = - 2e- t + 3e- 1000t.
В обоих случаях решение состоит из двух экспонент: быстро убывающей и относительно медленно изменяющейся. Отметим, что абсолютные величины собственных значений матриц рассматриваемых линейных систем ОДУ при их представлении в виде
$$$ \dot {u} = Au, $$$( u — вектор - столбец, A — матрица с
постоянными коэффициентами) существенно различаются. Так, в первом случае $$\lambda _{1} \approx (4\varepsilon )^{- 1},$$ $$\lambda _{2} \sim O(1)$$ ; во втором: $$\lambda _{1} \approx - 1,$$ $$\lambda _{2} = 10^{- 3}.$$ В обоих случаях имеем:
При моделировании физических процессов причина такой разницы в собственных числах заключена в существенно различных характерных временах процессов, описываемых системами ОДУ. Наиболее часто подобные системы встречаются при моделировании процессов в ядерных реакторах, при решении задач радиофизики, астрофизики, физики плазмы, биофизики, химической кинетики. Последние задачи часто могут быть записаны в виде [9.3]:
$$$ \frac{d u_k}{d t} = \sum\limits_{i = 1}^{N}{\sum\limits_{j = 1}^{N}{a_{_{ij}}^{k} u_i u_j}, } k = 1 \div N; $$$где uk — концентрации веществ, участвующих в химических реакциях, скорости протекания которых характеризуются коэффициентами $$a_{ij}^{k}.$$ В качестве примера приведем одну из систем химической кинетики, описывающую изменение концентрации трех веществ, участвующих в реакции для случая полного перемешивания [9.1].
Пример 1. Обозначим концентрации трех веществ, участвующих в
реакции, через u1, u2 и u3, тогда
Участки решения, характеризующиеся быстрым и медленным его изменением, называются пограничным слоем и квазистационарным режимом, соответственно.
Трудности численного решения подобных 1012 раз. Следовательно, если при численном решении системы
выбирать шаг из условия
$$\tau\|f^{\prime}_u (u)\|\ll 1,$$то он будет соответствовать самому быстрому процессу. В данном случае затраты машинного времени для исследования самых медленных процессов будут неоправданно велики. По этой причине имеются следующие альтернативы в выборе подхода к численному решению рассматриваемых задач.
$${\tau}\ll \|f^{\prime}_u (u)\|^{- 1},$$
т.е. с учетом характерных времен всех процессов, описываемых данной системой.
$${\tau}\gg \|f^{\prime}_u (u)\|^{- 1}.$$
Определение. ([9.3])
называется
J ={f'u(u)}
разделяется на две части.
$${\mathop{\mathrm{Re}}\nolimits}\Lambda_i (u) \le - \Lambda_0, \left|{\mathop{\mathrm{Im}}\nolimits}\Lambda_k\right| < \left|{\mathop{\mathrm{Re}}\nolimits}\Lambda_k\right|, k = 1 \div N_1 $$
( $$\Lambda _{i}$$ — собственные значения матрицы Якоби ); $$\Lambda _{0} > 0$$
$$|{\lambda_j}| \le \lambda_0, j = 1 \div N_2, \lambda_0 > 0.$$
При этом $$\lambda_0 \ll \Lambda_0 .$$
Отношение $$\Lambda _{0}/\lambda _{0}$$ называется показателем жесткости системы. В дальнейшем будем полагать $$\lambda _{0} \sim O(1).$$
Проблему численного решения
Ее точное решение задается формулой
$${{u}(t) = \sum\limits_{j = 1}^{N_1}{b_j e^{\Lambda_j t}{\Omega }_j} + \sum\limits_{k = 1}^{N_2 }{\bar {b}_k e^{\lambda_k t}{\omega }_k}, }$$где константы интегрирования $$b_j, \bar b_k$$ соответствуют
жесткой и мягкой частям
В этом решении видны две части: первая (жесткая) убывает как $$e^{- \Lambda_0 t}$$ на временном интервале $$[t_0, O(\Lambda_0^{- 1})]$$ (пограничный слой), вторая заметно изменяется на интервале $$[t_0, O(\Lambda_0^{- 1})]$$ (квазистационарный режим).
Если провести аппроксимацию линейной системы ОДУ с помощью явного метода Эйлера
$$$ \frac{u_{n + 1} - u_n}{{\tau}} = Bu_n, $$$или
$$u_{n + 1} = (E + \tau B)u_{n},$$то общее решение такой системы разностных уравнений будет иметь вид
$$u_n (t) = \sum\limits_{j = 1}^{N_1}{c_j(1 + {\tau}\Lambda_j)^n \Omega_j} + \sum\limits_{k = 1}^{N_2 }{c_k (1 + {\tau}\lambda_k)^{n}{\Omega_k}.$$Второе слагаемое в этом решении аппроксимирует второе слагаемое в точном решении (9.2), а первое быстро растет и приводит к абсурдному результату.
Теперь проведем аппроксимацию линейной системы ОДУ (9.1) с помощью
или
$$u_{n + 1} = (E + \tau B)^{ - 1}u_{n}.$$Общее решение такого разностного уравнения имеет следующий вид:
$$u_n (t) = \sum\limits_{j = 1}^{N_1 }{c_j (1 - {\tau}\Lambda_j)^{- n}{\Omega }_j} + \sum\limits_{k = 1}^{N_2 }{\bar c_k (1 - {\tau}\lambda_k)^{- n}{\omega }_k}.$$В этом случае второе слагаемое ведет себя так же, как и точное решение, а первое стремится к нулю как $$(\tau \Lambda _{0})^{- n},$$ т.е. его поведение качественно совпадает с точным в области пограничного слоя.
В практике численных исследований жестких задач часто не нужно изучать
поведение решения в пограничном слое, и можно воспользоваться
Устойчивость методов численного интегрирования
Положим, что численный метод, применяемый к решению этого уравнения, может быть записан в виде
$$u_{n + 1} = R(z)u_{n}, \ z = \tau \lambda ,$$где R(z) называется функцией устойчивости [9.1], [9.4]. О построении функции
устойчивости речь пойдет ниже.
Определение. Численный метод для решения уравнения (9.4) является абсолютно устойчивым, если выполнено условие
$$|R(z)| \le 1.$$Из определения следует, что $$|{u_{n + 1}}| \le |u_n|.$$
Это требование является естественным при $${\mathop{\mathrm{Re}}\nolimits} z \le 0,$$ поскольку в таком случае модуль точного решения есть невозрастающая функция.
Множество всех точек z, для которых $$|R(z)| \le 1,$$ называется областью абсолютной устойчивости.
(рис 9.2a) (рис 9.1) 
(рис 9.3) (рис 9.2b) Определение. Если область абсолютной устойчивости $$(|R(z)| \le
1)$$ занимает левую полуплоскость комплексной плоскости $$({Re} z \le
0),$$ то метод является
В случае, когда область абсолютной устойчивости включает в себя угол (в левой полуплоскости комплексной плоскости) с вершиной в нуле и углом полураствора $$\alpha,$$ то метод называется $$А(\alpha)$$ - устойчивым.
Области $$A(\alpha)$$ - и A(0) - устойчивости изображены на рис. 9.2а и 9.2b, соответственно.
В случае, когда вся область абсолютной устойчивости включает в себя часть левой полуплоскости (граница ее лежит вне заштрихованной части на рис. 9.3), то метод называется
Определение. Численный метод называется L - устойчивым, если он
В частности, рассмотренный выше L - устойчивым.
Решения, полученные такими методами, будут затухающими.
Рассмотрим простейшую нелинейную
или
$$\begin{gather*} \dot {u} = \frac{1}{\varepsilon } f(u, v), \\ \dot {v} = g(u, v), \end{gather*}$$с начальными условиями
u(0) = u0, v(0) = v0.
Из
т.е. $$\lambda _{1}$$ — жесткая, а $$\lambda _{2}$$
— мягкая часть
Первое, что хотелось бы сделать, — упростить задачу, положив $$\varepsilon \approx 0.$$ В этом случае система приобретает вид:
$$\begin{gather*} f(u, v) = 0, \\ \dot {v} = g(u, v), \\ u(0) = u_0, v(0) = v_0 \end{gather*} $$и описывает траекторию, проходящую по кривой f(u, v) = 0. Эта
кривая играет существенную роль в решении исходной системы ОДУ. Однако полностью поведение решения упрощенная (невозмущенная) система описывать не будет, поскольку начальные данные не обязательно должны лежать на этой кривой.
Для того чтобы понять эту ситуацию, рассмотрим фазовый портрет конкретной системы ОДУ вида:
$$\begin{gather*} \varepsilon \dot {u} = v - \frac{1}{3} u^3 + u, \\ \dot {v} = - u, \\ u(0) = u_0, v(0) = v_0 . \end{gather*}$$Эта система эквивалентна нелинейному уравнению второго порядка — уравнению Ван - дер - Поля. Приведенный выше вид иногда называется представлением Льенара. Подробнее о свойствах рассматриваемой задачи можно прочитать в [9.6], [9.7], [9.8], [9.9].
(рис 9.4) Кривая f(u, v) = 0 для этого случая изображена на рис. 9.4. Она делит плоскость {u, v} на две части: f > 0 и f < 0 ; вдали от кривой, поле скоростей {du, dv} направлено почти горизонтально влево (вправо) в зависимости от знака f. На самой кривой выделяются две устойчивые ветви AB и CD, где f'u < 0, и неустойчивая ветвь BD, на которой f'u > 0. Опишем теперь качественно поведение траектории рассматриваемой системы ОДУ, состоящей из следующих участков.
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, поэтому допустима оценка
откуда видно, что f(u, v0) стремится к нулю как экспонента с показателем $$$ \frac{1}{\varepsilon } f^{\prime}_u $,$$ т.е. f становится малой величиной ( $$\approx 0$$ ) за время $$t_{A'A} = O(\varepsilon ).$$
{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.
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*}$$
CD, как и по ветви AB, описывается уравнениями$$\begin{gather*} f(u, v) = 0, \\ \dot {v} = g(u, v), \end{gather*} $$
и длится $$t_{CD} \sim O(1).$$
DA также, как и на
ветви BC, происходит скачок из точки D за время $$t_{DA} \sim O(\varepsilon )$$ в устойчивую точку A, и т.д.Такое поведение траектории (замкнутая кривая) называется предельным циклом.
Для
Таким образом, это характерно для
Рассмотрим проблемы, которые могут возникнуть при численном интегрировании
подобных
В зоне квазистационарного режима часто оказывается, что интегрирование с
таким шагом слишком дорого. Можно, правда, разрешить уравнение f(u, v) =
0 $$(u = \varphi (v))$$ относительно u и далее
интегрировать его, предварительно реализовав алгоритм перехода на другой шаг:
однако в случаях более сложных, построение подобных численных методов может оказаться отдельной сложной задачей.
Чаще всего в практике численных расчетов целесообразно использовать
Эта
Первый пограничный слой будет пройден за один шаг:
$$\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. Возможно, что при
выборе большего шага интегрирования получится именно это нефизическое решение, к которому сойдутся итерации.
Такую опасность необходимо всегда учитывать при проведении численного интегрирования по
Рассмотрим один из наиболее простых численных методов решения жестких линейных систем ОДУ [9.10], основанный на представлении решения в явном виде
$$u_{n + 1} = e^{B\tau }u_{n}$$для линейных
Его численная реализация связана с вычислением матричной экспоненты. О
свойствах матричной экспоненты (и в целом функций от матриц) можно прочитать в [9.11]. Использование для этого
представляется непригодным, так как простые оценки показывают, что по
вычислительным затратам этот алгоритм сопоставим с
Следовательно, для того, чтобы k - й член разложения ряда
Тейлора был, по крайней мере, порядка O(1), необходимо выполнение условия
что соответствует условию численного интегрирования с шагом $$\tau || B || < 1\}.$$ Количество членов ряда при этом $$N \sim \tau\|{B}\| \gg 1,$$ что неприемлемо для решения.
Представим экспоненту в следующем эквивалентном виде:
$$$ e^{{B}{\tau}} = (e^{{\tau}A/2^p})^{2^p} \approx {\left[{E + \frac{{\tau}}{2^p}B + \ldots + \frac{1}{k!}{\left({\frac{{\tau}}{2^p}}\right)}^k B^k}\right]}^{2^p}.$$$При этом значение параметра p выбирают таким, чтобы
и можно было использовать ряд Тейлора с небольшим количеством членов. Действительно, в этом случае
$$$ \frac{{\tau\|{B}\|}}{{k \cdot 2^{p}}} \sim O(1) \mbox{ и } N \sim 2^{- p} \tau\|{B}\|, $$$что вполне приемлемо при соответствующем выборе параметра p.
В этом алгоритме сначала вычисляют матрицу
$$$ D = E + \sum\limits_{k = 1}^N{\left({\frac{\tau}{2^p}}\right)^k B^k}, $$$затем $$D^{2^p}$$ путем последовательных перемножений. При таком способе вычисления матрицы также имеется опасность, связанная с тем, что при некоторых p, $$\lambda _{0}, \Lambda _{0}$$ слагаемые, соответствующие мягкой части
а последовательность вычислений имеет вид
$$$ D = E + \frac{{\tau}}{2^p}B; C = D^{- 1}; C^{\prime} = {(C^2 )}^p . $$$При этом на мягкой части
на жесткой части, поскольку теперь p небольшие и $${\tau}|{\Lambda_j}|/2^p \gg 1,$$ получим оценку
что уже приемлемо для вычислений.
Рассмотрим другие методы для численного решения как линейных
и автономных
см. также [9.1], [9.4], [9.9], [9.13], [9.14], [9.15], [9.16]. Из соображений устойчивости метода предпочтение естественно отдать
$$$ \frac{u_{n + 1} - u_n}{{\tau}} = f(u_{n + 1}), u_0 = a; $$$
$$$ \frac{u_{n + 1} - u_n}{{\tau}} = \frac{f(u_{n + 1}, t_{n + 1}) + f(u_n, t_n)}{2}, u_0 = a;$$
$$$ \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 |
1/2 |
1 |
| $$$ \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 |
| $$$ \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} $$$ |
1 |
1 |
1 |
1/3 |
5/12 |
- 1/12 |
1 |
3/4 |
1/4 |
3/4 |
1/4 |
| $$$ \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} $$$ |
0 |
0 |
0 |
1 |
1/2 |
1/2 |
1/2 |
1/2 |
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 |
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 |
Рассмотрим теперь, как строится функция устойчивости для
Запишем теперь Y:
Приведенные выше формулы можно рассматривать как систему линейных уравнений
относительно новых переменных Y1, ..., Yk, vn + 1 следующего вида:
Необходимо выразить vn + 1 через vn. Для этого воспользуемся правилом Крамера. Ответ можно записать в следующей форме:
где 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. $$$Одноитерационные
Здесь $$$ B = \frac{\partial f}{\partial u}(u_n) $$$ матрица, постоянная на данном шаге по времени. Параметры a, b, c подбираются таким образом, чтобы обеспечить максимально возможный порядок точности.
Например, для схемы третьего порядка точности получим
a = 1, 077; b = - 0, 372; c = - 0, 577 .
Такую схему иногда называют методом с одной итерацией, имея в виду, что
вычисление обратной матрицы сравнимо по количеству арифметических операций с одной итерацией метода Ньютона. Преимущество
Рассмотрим связь между
Определение. ([9.9]) S - стадийный
где $$\beta, \gamma, \mu$$ — управляющие коэффициенты метода.
Как и выше, B вычисляется по
данным в точке tn. Связь последних формул с выражениями для явных
Несколько сложнее представляются
Вывод условий порядка
В предыдущей лекции вкратце был изложен альтернативный tn к tn + 1 использовались значения решения (или функций от него) в предыдущих точках. Полученные численные методы (первые схемы такого рода были получены Адамсом) носят название линейных многошаговых. Однако
применение
Теорема (Далквиста (или второй барьер Далквиста) [9.8], [9.9].)
Не существует
Попробуем ввести другой класс линейных многошаговых методов — это так называемые формулы дифференцирования назад (или ФДН - методы) [9.8].
Общий вид ФДН - метода таков:
$${u_{n + 1} + \sum\limits_{j = 1}^{k}{\alpha_j u_{n + 1 - j}} = {\tau}\beta f_{n + 1}, }$$где коэффициенты $$\alpha_j$$ выбираются из условий аппроксимации
метода. В отличие от un + 1 имеет меньшую размерность, следовательно, требуется меньшее число операций для нахождения решения.
Семейство $$A(\alpha)$$ - устойчивых ФДН - методов с достаточно
большим значением угла полураствора $$\alpha$$ носит общее название
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$$и введем также в рассмотрение вектор Нордсика
$$$ 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:
для записи ряда использован остаточный член в форме Лагранжа. Кроме того, используем следующие следствия приведенного выше разложения:
$$$ {{\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), } $$$$$$ {\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), } $$$
$$$ {\frac{\tau^3}{3!}u^{(3)}_{n + 1} = 3\frac{\tau^3}{3!}u_n^{(3)} + 4\frac{\tau^4}{4!}u^{(4)}(\xi)} $$$а конкретное значение $$\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}), } $$$а $$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), добиться максимальной устойчивости метода при приемлемом порядке аппроксимации.
Окончательный ответ для
где матрица 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)) значения производных решения. После нетрудных преобразований получим формулу Гира. Другой рассмотренный вариант метода Нордсика приведет к одной из
Отметим, что z0 = (u0, 0, 0, ..., 0), что позволит начать вычисления, но приведет к уменьшению порядка аппроксимации. Таким образом, многозначный (по введенной в [9.2] терминологии) вариант
Отметим также, что для системы с релаксационными колебаниями лучшие результаты могут давать многозначные
Простейшая математическая модель периодической химической реакции Белоусова - Жаботинского состоит из трех уравнений:
$$\begin{gather*} \dot {y}_1 = 77, 27(y_2 + y_1 (1 - 8, 375 \cdot 10^{- 6} y_1 - y_2 )), \\ \dot {y}_2 = \frac{1}{{77, 27}}(y_3 - (1 + y_1 )y_2), \\ \dot {y}_3 = 0, 161(y_1 - y_3 ) \end{gather*}$$На то, что
Так как переменные системы — концентрации ( HBrO2, Br - и Ce(IV) соответственно) то начальные условия для системы следует выбирать положительными, как правило, близкими к 0. Конечное время интегрирования системы Tk = 800.
Типичным примером жесткой задачи малой размерности является уравнение Ван - дер - Поля [9.8], [9.9], [9.18], [9.19]. Его возможно записать в виде системы
$$\begin{gather*} y^{\prime}_1 = y_2, \\ {y^{\prime}_2 = - a(y_2 (y_1^2 - 1) + y_1 ), } \end{gather*}$$или в виде
$$\begin{gather*} y^{\prime}_1 = - a(\frac{y_1^3}{3} - y_1) + ay_2, \\ y^{\prime}_2 = - y_1, \end{gather*}$$(представление Льенара). Считаем, что параметр a — большой. В расчетах рассмотреть два случая: a = 103 и a = 106. Для тестов обычно полагают y1 = 2, y2 = 0.
Конечное время интегрирования системы, записанной в виде (9.13), Tk = 20.
Периодические решения
Дополнительный вопрос: указать преобразование, переводящее представление (9.13) в представление Льенара (9.14).
Рассмотрим неавтономную систему уравнений Ван - дер - Поля:
$$\begin{gather*} y^{\prime}_1 = a(- (\frac{y_1^3}{3} - y_1 ) + y_2), \\ y^{\prime}_2 = - y_1 + A\cos \omega t \end{gather*}$$Как и в предыдущей задаче считаем, что a = 103 и a = 106, y1 = 2, y2 = 0. Рассмотреть численно случаи 0 < A < 1 и $$$ 1 < A < \sqrt{1 + \frac{1}{64\omega ^2 }} $$$ Tk = 200.
О траекториях - утках в системе Ван - дер - Поля см. [9.19] (строгое математическое исследование) и [9.20](популярное изложение).
Рассмотрим простейшую математическую модель колебаний концентрации озона в атмосфере [9.2]. Она описывается следующей неавтономной системой ОДУ:
$$\begin{gather*} \dot {y_1} = - k_1 y_1 y_2 - k_2 y_1 y_3 + 2k_3 (t)y_2 + k_4 (t)y_3 , \\ \dot {y_2} = 0, \\ \dot {y_1} = k_1 y_1 y_2 - k_2 y_1 y_3 - k_4 (t)y_3 \end{gather*}$$В данной модели уравнения описывают изменение концентрации атомарного кислорода, молекулярного кислорода и озона соответственно. Считается, что изменения концентрации молекулярного кислорода невелики. Начальные значения для задачи таковы:
$$y_1 (0) = 10^6 (\mbox{см}^{- 3}), y_2 (0) = 3, 7 \cdot 10^{16}(\mbox{см}^{- 3}), y_3 (0) = 10^{12}(\mbox{см}^{- 3}),$$значения констант скоростей химических реакций
$$k_1 = 1, 63 \cdot 10^{- 16}, \quad k_2 = 4, 66 \cdot 10^{- 16}.$$Две другие химические реакции зависят от локальной освещенности участка земной поверхности и приближаются следующим выражением:
$$k_i (t) = \left\{ \begin{array}{cc} \exp (- c_i /\sin \omega t), \sin \omega t > 0, \\ 0, \sin \omega t < 0, \\ \end{array} \right.$$где $$\omega = \pi /43 200c^{- 1},$$ c3 = 22, 62, c4 = 7, 601. Значения констант скоростей обращаются в нуль ночью, резко возрастают на рассвете, достигают максимума в полдень и падают до нуля на закате. Конечное время интегрирования Tk = 172800 c (двое суток).
Данная
Рассмотрим еще один пример жесткой задачи малой размерности, имеющей периодическое решение [9.19], [9.21].
$$\begin{gather*} y^{\prime}_1 = a(- (\frac{y_1^3}{3} - y_1 ) + y_2), \\ y^{\prime}_2 = - y_1 - by_2 + c \end{gather*}$$Здесь a = 103 и a = 106, y1 = 2, y2= 0.
Уравнение описывает протекание тока через клеточную мембрану. Постоянная компонента тока c в безразмерной записи системы такова, что 0 < c < 1, b > 0. Tk = 20.
Система более высокой размерности, имеющая решение в виде релаксационного цикла, приведена в [9.18] (см. также [9.21]). Она имеет вид:
$$\begin{gather*} \varepsilon \dot {x}_1 = - \alpha (y_1 - y_2 ) + \varphi (x_1) - x_2, \\ \varepsilon \dot {x}_2 = \alpha (y_1 - y_2 ) + \varphi (x_2 ) - x_1, \\ \dot {y}_1 = x_1, \\ \dot {y}_2 = x_2 . \end{gather*} $$Здесь $$\alpha > 0$$ — константа порядка единицы, функция $$\varphi (u) = {- \tg(\pi u/2)},$$ x1(0) = x2(0) = 0, y1 = 2, y2 = 0, Tk = 20, $$\varepsilon = 10^{- 3}, 10^{- 6}.$$
Простейшая модель гликолиза описывается уравнениями следующего вида [9.21]:
$$\begin{gather*} \dot {y}_1 = 1 - y_1 y_2, \\ \dot {y}_2 = \alpha y_2 \left({y_1 - \frac{{1 + \beta }}{{y_2 + \beta }}}\right), \end{gather*}$$предложенными Дж. Хиггинсом. В системе $$\beta = 10, \alpha = 100,$$ 200, 400, 1000. Начальные условия
для системы: y1(0) = 1, y2(0) = 0, 001, Tk = 50. Решение этой системы — релаксационные автоколебания (жесткий предельный цикл).
Один из первых и самых популярных примеров
Система Робертсона имеет вид [9.9]
$$\begin{gather*} \dot {y}_1 = - 0.04y_1 + 10^4 y_2 y_3, \\ \dot {y}_2 = 0.04y_1 - 10^4 y_2 y_3 - 3 \cdot 10^7 y_2^2, \\ \dot {y}_3 = 3 \cdot 10^7 y_2^2 . \end{gather*} $$Начальные условия для системы таковы: y1(0) = 1, y2(0) = 0, y3(0) = 0. Рассматриваются следующие величины отрезка интегрирования: Tk = 40 (в работе Робертсона рассматривался именно такой отрезок интегрирования), Tk = 100, 1000, ..., 1011. О свойствах
задачи см. в [9.9].
Данный пример из [9.9] — типичный случай биохимической модели "умеренной" размерности (современные модели, например, фотосинтеза включают сотни уравнений подобного типа). Хотя данная
модель является умеренно жесткой, тем не менее, ее лучше решать с помощью методов, предназначенных для решения
Начальные значения всех переменных системы равны 0, кроме y1(0) = 1 и y8(0) = 0.0057. Длина отрезка интегрирования Tk = 421, 8122.
Еще одна модель химической реакции из [9.9], получившая свое название Е5 в более ранних публикациях.
$$\begin{gather*} \dot {y}_1 = - {Ay}_1 - {By}_1 y_3, \\ \dot {y}_2 = {Ay}_1 - {MCy}_2 y_3, \\ \dot {y}_3 = {Ay}_1 - {By}_1 y_3 - {MCy}_2 y_3 + {Cy}_4, \\ \dot {y}_4 = {By}_1 y_3 - {Cy}_4 . \end{gather*} $$Начальные условия: $$y_1(0) = 1, 76\cdot 10^{- 3},$$ а все остальные
переменные равны 0. Значения коэффициентов модели следующие: $$A = 7, 89 \cdot 10^{- 10}, B = 1, 1\cdot 10^7, C = 1, 13\cdot 10^3, M = 10^6.$$ Первоначально задача ставилась на отрезке Tk = 1000, но впоследствии было обнаружено, что она обладает нетривиальными свойствами вплоть до времени Tk = 1013 (подробнее см. [9.9]).
Обратить особое внимание, что в процессе расчетов приходится иметь дело с очень малыми концентрациями реагентов (малы значения y2, y3 и y4 ). Как "подправить" постановку задачи E5?
Уравнение Релея во многом похоже на уравнение Ван - дер - Поля [9.21]. Рассматривается задача вида
$$$ \ddot {x} - \mu (1 - \dot {x}^2 )\dot {x} + x = 0 . $$$Решить задачу, записав уравнение Релея в виде системы ОДУ. Начальные условия: $$$ x(0) = 0, \dot {x}(0) = 0, 001 $,$$ $$\mu = 1000,$$ Tk = 1000.
Рассмотрим пример системы уравнений, которая описывает изменения численности популяций двух видов и эволюцию некого генетического признака $$\alpha.$$ Система ОДУ имеет вид
$$\begin{gather*} \dot {x} = x(1 - 0, 5x - \frac{2}{{7\alpha ^2}}y), \\ \dot {y} = y(2\alpha - 3, 5\alpha ^2 x - 0, 5y), \\ \dot {\alpha} = \varepsilon (2 - 7\alpha x) . \end{gather*}$$Параметры задачи таковы: $$\varepsilon \le 0, 01, 0 \le x_0 \le 3, 0 \le y_0 \le 15,
\alpha_0 = 0,$$ Tk = 1500. Наличие малого параметра в третьем уравнении системы показывает, что генетический признак меняется медленнее, чем численность популяций. Решение системы — релаксационные колебания.
Задача описана в статье [9.22].
Еще один пример
Параметры задачи таковы: $$\varepsilon \le 0, 01, 0 \le x_0 \le 40, 0 \le y_0 \le 40,
\alpha_{10} = 0, \alpha_{20} = 10,$$ Tk = 2000.
Рассмотреть также модификацию предыдущей системы [9.22]:
$$\begin{gather*} \dot {x} = x(2\alpha_1 - 0, 5x - \alpha_1^3 \alpha_2^{- 3}y), \\ \dot {y} = y(2\alpha_2 - \alpha_1^{- 3} \alpha_{2 }^3 x - 0, 5y), \\ \dot {\alpha}_1 = \varepsilon (2 - 3\alpha_1^2 \alpha_2^{- 3}y), \\ \dot {\alpha}_2 = \varepsilon (2 - 3\alpha_1^{- 3} \alpha_2^2x) . \end{gather*} $$Параметры задачи: $$\varepsilon \le 0, 01, 0 \le x_0 \le 40, 0 \le y_0 \le 40,
\alpha_{10} = 0, \alpha_{20} = 10,$$ Tk = 2000.
Рассмотрим в качестве примера две задачи Коши для систем обыкновенных дифференциальных уравнений (ОДУ) [9.1], [9.2]:
$$$ \dot {u} = au + \frac{1}{\varepsilon } v, \dot {v} = - \frac{1}{\varepsilon } v, $$$с начальными данными u(0) = u0, v(0) = v0 ; здесь $$a \sim O(1), \varepsilon \ll 1$$ ; и линейную систему с постоянными коэффициентами
Решением первой задачи Коши являются функции
$$\begin{gather*} u(t) = u_0 e^{at} + \frac{{v_0 }}{{1 + a\varepsilon }} (e^{at} - e^{- t/\varepsilon }), \\ v(t) = v_0 e^{- t/\varepsilon }, \end{gather*}$$а второй -
u(t) = 4e- t - 3e- 1000t, v(t) = - 2e- t + 3e- 1000t.
В обоих случаях решение состоит из двух экспонент: быстро убывающей и относительно медленно изменяющейся. Отметим, что абсолютные величины собственных значений матриц рассматриваемых линейных систем ОДУ при их представлении в виде
$$$ \dot {u} = Au, $$$( u — вектор - столбец, A — матрица с
постоянными коэффициентами) существенно различаются. Так, в первом случае $$\lambda _{1} \approx (4\varepsilon )^{- 1},$$ $$\lambda _{2} \sim O(1)$$ ; во втором: $$\lambda _{1} \approx - 1,$$ $$\lambda _{2} = 10^{- 3}.$$ В обоих случаях имеем:
При моделировании физических процессов причина такой разницы в собственных числах заключена в существенно различных характерных временах процессов, описываемых системами ОДУ. Наиболее часто подобные системы встречаются при моделировании процессов в ядерных реакторах, при решении задач радиофизики, астрофизики, физики плазмы, биофизики, химической кинетики. Последние задачи часто могут быть записаны в виде [9.3]:
$$$ \frac{d u_k}{d t} = \sum\limits_{i = 1}^{N}{\sum\limits_{j = 1}^{N}{a_{_{ij}}^{k} u_i u_j}, } k = 1 \div N; $$$где uk — концентрации веществ, участвующих в химических реакциях, скорости протекания которых характеризуются коэффициентами $$a_{ij}^{k}.$$ В качестве примера приведем одну из систем химической кинетики, описывающую изменение концентрации трех веществ, участвующих в реакции для случая полного перемешивания [9.1].
Пример 1. Обозначим концентрации трех веществ, участвующих в
реакции, через u1, u2 и u3, тогда
Участки решения, характеризующиеся быстрым и медленным его изменением, называются пограничным слоем и квазистационарным режимом, соответственно.
Трудности численного решения подобных 1012 раз. Следовательно, если при численном решении системы
выбирать шаг из условия
$$\tau\|f^{\prime}_u (u)\|\ll 1,$$то он будет соответствовать самому быстрому процессу. В данном случае затраты машинного времени для исследования самых медленных процессов будут неоправданно велики. По этой причине имеются следующие альтернативы в выборе подхода к численному решению рассматриваемых задач.
$${\tau}\ll \|f^{\prime}_u (u)\|^{- 1},$$
т.е. с учетом характерных времен всех процессов, описываемых данной системой.
$${\tau}\gg \|f^{\prime}_u (u)\|^{- 1}.$$
Определение. ([9.3])
называется
J ={f'u(u)}
разделяется на две части.
$${\mathop{\mathrm{Re}}\nolimits}\Lambda_i (u) \le - \Lambda_0, \left|{\mathop{\mathrm{Im}}\nolimits}\Lambda_k\right| < \left|{\mathop{\mathrm{Re}}\nolimits}\Lambda_k\right|, k = 1 \div N_1 $$
( $$\Lambda _{i}$$ — собственные значения матрицы Якоби ); $$\Lambda _{0} > 0$$
$$|{\lambda_j}| \le \lambda_0, j = 1 \div N_2, \lambda_0 > 0.$$
При этом $$\lambda_0 \ll \Lambda_0 .$$
Отношение $$\Lambda _{0}/\lambda _{0}$$ называется показателем жесткости системы. В дальнейшем будем полагать $$\lambda _{0} \sim O(1).$$
Проблему численного решения
Ее точное решение задается формулой
$${{u}(t) = \sum\limits_{j = 1}^{N_1}{b_j e^{\Lambda_j t}{\Omega }_j} + \sum\limits_{k = 1}^{N_2 }{\bar {b}_k e^{\lambda_k t}{\omega }_k}, }$$где константы интегрирования $$b_j, \bar b_k$$ соответствуют
жесткой и мягкой частям
В этом решении видны две части: первая (жесткая) убывает как $$e^{- \Lambda_0 t}$$ на временном интервале $$[t_0, O(\Lambda_0^{- 1})]$$ (пограничный слой), вторая заметно изменяется на интервале $$[t_0, O(\Lambda_0^{- 1})]$$ (квазистационарный режим).
Если провести аппроксимацию линейной системы ОДУ с помощью явного метода Эйлера
$$$ \frac{u_{n + 1} - u_n}{{\tau}} = Bu_n, $$$или
$$u_{n + 1} = (E + \tau B)u_{n},$$то общее решение такой системы разностных уравнений будет иметь вид
$$u_n (t) = \sum\limits_{j = 1}^{N_1}{c_j(1 + {\tau}\Lambda_j)^n \Omega_j} + \sum\limits_{k = 1}^{N_2 }{c_k (1 + {\tau}\lambda_k)^{n}{\Omega_k}.$$Второе слагаемое в этом решении аппроксимирует второе слагаемое в точном решении (9.2), а первое быстро растет и приводит к абсурдному результату.
Теперь проведем аппроксимацию линейной системы ОДУ (9.1) с помощью
или
$$u_{n + 1} = (E + \tau B)^{ - 1}u_{n}.$$Общее решение такого разностного уравнения имеет следующий вид:
$$u_n (t) = \sum\limits_{j = 1}^{N_1 }{c_j (1 - {\tau}\Lambda_j)^{- n}{\Omega }_j} + \sum\limits_{k = 1}^{N_2 }{\bar c_k (1 - {\tau}\lambda_k)^{- n}{\omega }_k}.$$В этом случае второе слагаемое ведет себя так же, как и точное решение, а первое стремится к нулю как $$(\tau \Lambda _{0})^{- n},$$ т.е. его поведение качественно совпадает с точным в области пограничного слоя.
В практике численных исследований жестких задач часто не нужно изучать
поведение решения в пограничном слое, и можно воспользоваться
Устойчивость методов численного интегрирования
Положим, что численный метод, применяемый к решению этого уравнения, может быть записан в виде
$$u_{n + 1} = R(z)u_{n}, \ z = \tau \lambda ,$$где R(z) называется функцией устойчивости [9.1], [9.4]. О построении функции
устойчивости речь пойдет ниже.
Определение. Численный метод для решения уравнения (9.4) является абсолютно устойчивым, если выполнено условие
$$|R(z)| \le 1.$$Из определения следует, что $$|{u_{n + 1}}| \le |u_n|.$$
Это требование является естественным при $${\mathop{\mathrm{Re}}\nolimits} z \le 0,$$ поскольку в таком случае модуль точного решения есть невозрастающая функция.
Множество всех точек z, для которых $$|R(z)| \le 1,$$ называется областью абсолютной устойчивости.
(рис 9.2a) (рис 9.1) (рис 9.3) (рис 9.2b) Определение. Если область абсолютной устойчивости $$(|R(z)| \le
1)$$ занимает левую полуплоскость комплексной плоскости $$({Re} z \le
0),$$ то метод является
В случае, когда область абсолютной устойчивости включает в себя угол (в левой полуплоскости комплексной плоскости) с вершиной в нуле и углом полураствора $$\alpha,$$ то метод называется $$А(\alpha)$$ - устойчивым.
Области $$A(\alpha)$$ - и A(0) - устойчивости изображены на рис. 9.2а и 9.2b, соответственно.
В случае, когда вся область абсолютной устойчивости включает в себя часть левой полуплоскости (граница ее лежит вне заштрихованной части на рис. 9.3), то метод называется
Определение. Численный метод называется L - устойчивым, если он
В частности, рассмотренный выше L - устойчивым.
Решения, полученные такими методами, будут затухающими.
Рассмотрим простейшую нелинейную
или
$$\begin{gather*} \dot {u} = \frac{1}{\varepsilon } f(u, v), \\ \dot {v} = g(u, v), \end{gather*}$$с начальными условиями
u(0) = u0, v(0) = v0.
Из
т.е. $$\lambda _{1}$$ — жесткая, а $$\lambda _{2}$$
— мягкая часть
Первое, что хотелось бы сделать, — упростить задачу, положив $$\varepsilon \approx 0.$$ В этом случае система приобретает вид:
$$\begin{gather*} f(u, v) = 0, \\ \dot {v} = g(u, v), \\ u(0) = u_0, v(0) = v_0 \end{gather*} $$и описывает траекторию, проходящую по кривой f(u, v) = 0. Эта
кривая играет существенную роль в решении исходной системы ОДУ. Однако полностью поведение решения упрощенная (невозмущенная) система описывать не будет, поскольку начальные данные не обязательно должны лежать на этой кривой.
Для того чтобы понять эту ситуацию, рассмотрим фазовый портрет конкретной системы ОДУ вида:
$$\begin{gather*} \varepsilon \dot {u} = v - \frac{1}{3} u^3 + u, \\ \dot {v} = - u, \\ u(0) = u_0, v(0) = v_0 . \end{gather*}$$Эта система эквивалентна нелинейному уравнению второго порядка — уравнению Ван - дер - Поля. Приведенный выше вид иногда называется представлением Льенара. Подробнее о свойствах рассматриваемой задачи можно прочитать в [9.6], [9.7], [9.8], [9.9].
(рис 9.4) Кривая f(u, v) = 0 для этого случая изображена на рис. 9.4. Она делит плоскость {u, v} на две части: f > 0 и f < 0 ; вдали от кривой, поле скоростей {du, dv} направлено почти горизонтально влево (вправо) в зависимости от знака f. На самой кривой выделяются две устойчивые ветви AB и CD, где f'u < 0, и неустойчивая ветвь BD, на которой f'u > 0. Опишем теперь качественно поведение траектории рассматриваемой системы ОДУ, состоящей из следующих участков.
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, поэтому допустима оценка
откуда видно, что f(u, v0) стремится к нулю как экспонента с показателем $$$ \frac{1}{\varepsilon } f^{\prime}_u $,$$ т.е. f становится малой величиной ( $$\approx 0$$ ) за время $$t_{A'A} = O(\varepsilon ).$$
{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.
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*}$$
CD, как и по ветви AB, описывается уравнениями$$\begin{gather*} f(u, v) = 0, \\ \dot {v} = g(u, v), \end{gather*} $$
и длится $$t_{CD} \sim O(1).$$
DA также, как и на
ветви BC, происходит скачок из точки D за время $$t_{DA} \sim O(\varepsilon )$$ в устойчивую точку A, и т.д.Такое поведение траектории (замкнутая кривая) называется предельным циклом.
Для
Таким образом, это характерно для
Рассмотрим проблемы, которые могут возникнуть при численном интегрировании
подобных
В зоне квазистационарного режима часто оказывается, что интегрирование с
таким шагом слишком дорого. Можно, правда, разрешить уравнение f(u, v) =
0 $$(u = \varphi (v))$$ относительно u и далее
интегрировать его, предварительно реализовав алгоритм перехода на другой шаг:
однако в случаях более сложных, построение подобных численных методов может оказаться отдельной сложной задачей.
Чаще всего в практике численных расчетов целесообразно использовать
Эта
Первый пограничный слой будет пройден за один шаг:
$$\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. Возможно, что при
выборе большего шага интегрирования получится именно это нефизическое решение, к которому сойдутся итерации.
Такую опасность необходимо всегда учитывать при проведении численного интегрирования по
Рассмотрим один из наиболее простых численных методов решения жестких линейных систем ОДУ [9.10], основанный на представлении решения в явном виде
$$u_{n + 1} = e^{B\tau }u_{n}$$для линейных
Его численная реализация связана с вычислением матричной экспоненты. О
свойствах матричной экспоненты (и в целом функций от матриц) можно прочитать в [9.11]. Использование для этого
представляется непригодным, так как простые оценки показывают, что по
вычислительным затратам этот алгоритм сопоставим с
Следовательно, для того, чтобы k - й член разложения ряда
Тейлора был, по крайней мере, порядка O(1), необходимо выполнение условия
что соответствует условию численного интегрирования с шагом $$\tau || B || < 1\}.$$ Количество членов ряда при этом $$N \sim \tau\|{B}\| \gg 1,$$ что неприемлемо для решения.
Представим экспоненту в следующем эквивалентном виде:
$$$ e^{{B}{\tau}} = (e^{{\tau}A/2^p})^{2^p} \approx {\left[{E + \frac{{\tau}}{2^p}B + \ldots + \frac{1}{k!}{\left({\frac{{\tau}}{2^p}}\right)}^k B^k}\right]}^{2^p}.$$$При этом значение параметра p выбирают таким, чтобы
и можно было использовать ряд Тейлора с небольшим количеством членов. Действительно, в этом случае
$$$ \frac{{\tau\|{B}\|}}{{k \cdot 2^{p}}} \sim O(1) \mbox{ и } N \sim 2^{- p} \tau\|{B}\|, $$$что вполне приемлемо при соответствующем выборе параметра p.
В этом алгоритме сначала вычисляют матрицу
$$$ D = E + \sum\limits_{k = 1}^N{\left({\frac{\tau}{2^p}}\right)^k B^k}, $$$затем $$D^{2^p}$$ путем последовательных перемножений. При таком способе вычисления матрицы также имеется опасность, связанная с тем, что при некоторых p, $$\lambda _{0}, \Lambda _{0}$$ слагаемые, соответствующие мягкой части
а последовательность вычислений имеет вид
$$$ D = E + \frac{{\tau}}{2^p}B; C = D^{- 1}; C^{\prime} = {(C^2 )}^p . $$$При этом на мягкой части
на жесткой части, поскольку теперь p небольшие и $${\tau}|{\Lambda_j}|/2^p \gg 1,$$ получим оценку
что уже приемлемо для вычислений.
Рассмотрим другие методы для численного решения как линейных
и автономных
см. также [9.1], [9.4], [9.9], [9.13], [9.14], [9.15], [9.16]. Из соображений устойчивости метода предпочтение естественно отдать
$$$ \frac{u_{n + 1} - u_n}{{\tau}} = f(u_{n + 1}), u_0 = a; $$$
$$$ \frac{u_{n + 1} - u_n}{{\tau}} = \frac{f(u_{n + 1}, t_{n + 1}) + f(u_n, t_n)}{2}, u_0 = a;$$
$$$ \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 |
1/2 |
1 |
| $$$ \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 |
| $$$ \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} $$$ |
1 |
1 |
1 |
1/3 |
5/12 |
- 1/12 |
1 |
3/4 |
1/4 |
3/4 |
1/4 |
| $$$ \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} $$$ |
0 |
0 |
0 |
1 |
1/2 |
1/2 |
1/2 |
1/2 |
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 |
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 |
Рассмотрим теперь, как строится функция устойчивости для
Запишем теперь Y:
Приведенные выше формулы можно рассматривать как систему линейных уравнений
относительно новых переменных Y1, ..., Yk, vn + 1 следующего вида:
Необходимо выразить vn + 1 через vn. Для этого воспользуемся правилом Крамера. Ответ можно записать в следующей форме:
где 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. $$$Одноитерационные
Здесь $$$ B = \frac{\partial f}{\partial u}(u_n) $$$ матрица, постоянная на данном шаге по времени. Параметры a, b, c подбираются таким образом, чтобы обеспечить максимально возможный порядок точности.
Например, для схемы третьего порядка точности получим
a = 1, 077; b = - 0, 372; c = - 0, 577 .
Такую схему иногда называют методом с одной итерацией, имея в виду, что
вычисление обратной матрицы сравнимо по количеству арифметических операций с одной итерацией метода Ньютона. Преимущество
Рассмотрим связь между
Определение. ([9.9]) S - стадийный
где $$\beta, \gamma, \mu$$ — управляющие коэффициенты метода.
Как и выше, B вычисляется по
данным в точке tn. Связь последних формул с выражениями для явных
Несколько сложнее представляются
Вывод условий порядка
В предыдущей лекции вкратце был изложен альтернативный tn к tn + 1 использовались значения решения (или функций от него) в предыдущих точках. Полученные численные методы (первые схемы такого рода были получены Адамсом) носят название линейных многошаговых. Однако
применение
Теорема (Далквиста (или второй барьер Далквиста) [9.8], [9.9].)
Не существует
Попробуем ввести другой класс линейных многошаговых методов — это так называемые формулы дифференцирования назад (или ФДН - методы) [9.8].
Общий вид ФДН - метода таков:
$${u_{n + 1} + \sum\limits_{j = 1}^{k}{\alpha_j u_{n + 1 - j}} = {\tau}\beta f_{n + 1}, }$$где коэффициенты $$\alpha_j$$ выбираются из условий аппроксимации
метода. В отличие от un + 1 имеет меньшую размерность, следовательно, требуется меньшее число операций для нахождения решения.
Семейство $$A(\alpha)$$ - устойчивых ФДН - методов с достаточно
большим значением угла полураствора $$\alpha$$ носит общее название
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$$и введем также в рассмотрение вектор Нордсика
$$$ 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:
для записи ряда использован остаточный член в форме Лагранжа. Кроме того, используем следующие следствия приведенного выше разложения:
$$$ {{\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), } $$$$$$ {\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), } $$$
$$$ {\frac{\tau^3}{3!}u^{(3)}_{n + 1} = 3\frac{\tau^3}{3!}u_n^{(3)} + 4\frac{\tau^4}{4!}u^{(4)}(\xi)} $$$а конкретное значение $$\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}), } $$$а $$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), добиться максимальной устойчивости метода при приемлемом порядке аппроксимации.
Окончательный ответ для
где матрица 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)) значения производных решения. После нетрудных преобразований получим формулу Гира. Другой рассмотренный вариант метода Нордсика приведет к одной из
Отметим, что z0 = (u0, 0, 0, ..., 0), что позволит начать вычисления, но приведет к уменьшению порядка аппроксимации. Таким образом, многозначный (по введенной в [9.2] терминологии) вариант
Отметим также, что для системы с релаксационными колебаниями лучшие результаты могут давать многозначные
Простейшая математическая модель периодической химической реакции Белоусова - Жаботинского состоит из трех уравнений:
$$\begin{gather*} \dot {y}_1 = 77, 27(y_2 + y_1 (1 - 8, 375 \cdot 10^{- 6} y_1 - y_2 )), \\ \dot {y}_2 = \frac{1}{{77, 27}}(y_3 - (1 + y_1 )y_2), \\ \dot {y}_3 = 0, 161(y_1 - y_3 ) \end{gather*}$$На то, что
Так как переменные системы — концентрации ( HBrO2, Br - и Ce(IV) соответственно) то начальные условия для системы следует выбирать положительными, как правило, близкими к 0. Конечное время интегрирования системы Tk = 800.
Типичным примером жесткой задачи малой размерности является уравнение Ван - дер - Поля [9.8], [9.9], [9.18], [9.19]. Его возможно записать в виде системы
$$\begin{gather*} y^{\prime}_1 = y_2, \\ {y^{\prime}_2 = - a(y_2 (y_1^2 - 1) + y_1 ), } \end{gather*}$$или в виде
$$\begin{gather*} y^{\prime}_1 = - a(\frac{y_1^3}{3} - y_1) + ay_2, \\ y^{\prime}_2 = - y_1, \end{gather*}$$(представление Льенара). Считаем, что параметр a — большой. В расчетах рассмотреть два случая: a = 103 и a = 106. Для тестов обычно полагают y1 = 2, y2 = 0.
Конечное время интегрирования системы, записанной в виде (9.13), Tk = 20.
Периодические решения
Дополнительный вопрос: указать преобразование, переводящее представление (9.13) в представление Льенара (9.14).
Рассмотрим неавтономную систему уравнений Ван - дер - Поля:
$$\begin{gather*} y^{\prime}_1 = a(- (\frac{y_1^3}{3} - y_1 ) + y_2), \\ y^{\prime}_2 = - y_1 + A\cos \omega t \end{gather*}$$Как и в предыдущей задаче считаем, что a = 103 и a = 106, y1 = 2, y2 = 0. Рассмотреть численно случаи 0 < A < 1 и $$$ 1 < A < \sqrt{1 + \frac{1}{64\omega ^2 }} $$$ Tk = 200.
О траекториях - утках в системе Ван - дер - Поля см. [9.19] (строгое математическое исследование) и [9.20](популярное изложение).
Рассмотрим простейшую математическую модель колебаний концентрации озона в атмосфере [9.2]. Она описывается следующей неавтономной системой ОДУ:
$$\begin{gather*} \dot {y_1} = - k_1 y_1 y_2 - k_2 y_1 y_3 + 2k_3 (t)y_2 + k_4 (t)y_3 , \\ \dot {y_2} = 0, \\ \dot {y_1} = k_1 y_1 y_2 - k_2 y_1 y_3 - k_4 (t)y_3 \end{gather*}$$В данной модели уравнения описывают изменение концентрации атомарного кислорода, молекулярного кислорода и озона соответственно. Считается, что изменения концентрации молекулярного кислорода невелики. Начальные значения для задачи таковы:
$$y_1 (0) = 10^6 (\mbox{см}^{- 3}), y_2 (0) = 3, 7 \cdot 10^{16}(\mbox{см}^{- 3}), y_3 (0) = 10^{12}(\mbox{см}^{- 3}),$$значения констант скоростей химических реакций
$$k_1 = 1, 63 \cdot 10^{- 16}, \quad k_2 = 4, 66 \cdot 10^{- 16}.$$Две другие химические реакции зависят от локальной освещенности участка земной поверхности и приближаются следующим выражением:
$$k_i (t) = \left\{ \begin{array}{cc} \exp (- c_i /\sin \omega t), \sin \omega t > 0, \\ 0, \sin \omega t < 0, \\ \end{array} \right.$$где $$\omega = \pi /43 200c^{- 1},$$ c3 = 22, 62, c4 = 7, 601. Значения констант скоростей обращаются в нуль ночью, резко возрастают на рассвете, достигают максимума в полдень и падают до нуля на закате. Конечное время интегрирования Tk = 172800 c (двое суток).
Данная
Рассмотрим еще один пример жесткой задачи малой размерности, имеющей периодическое решение [9.19], [9.21].
$$\begin{gather*} y^{\prime}_1 = a(- (\frac{y_1^3}{3} - y_1 ) + y_2), \\ y^{\prime}_2 = - y_1 - by_2 + c \end{gather*}$$Здесь a = 103 и a = 106, y1 = 2, y2= 0.
Уравнение описывает протекание тока через клеточную мембрану. Постоянная компонента тока c в безразмерной записи системы такова, что 0 < c < 1, b > 0. Tk = 20.
Система более высокой размерности, имеющая решение в виде релаксационного цикла, приведена в [9.18] (см. также [9.21]). Она имеет вид:
$$\begin{gather*} \varepsilon \dot {x}_1 = - \alpha (y_1 - y_2 ) + \varphi (x_1) - x_2, \\ \varepsilon \dot {x}_2 = \alpha (y_1 - y_2 ) + \varphi (x_2 ) - x_1, \\ \dot {y}_1 = x_1, \\ \dot {y}_2 = x_2 . \end{gather*} $$Здесь $$\alpha > 0$$ — константа порядка единицы, функция $$\varphi (u) = {- \tg(\pi u/2)},$$ x1(0) = x2(0) = 0, y1 = 2, y2 = 0, Tk = 20, $$\varepsilon = 10^{- 3}, 10^{- 6}.$$
Простейшая модель гликолиза описывается уравнениями следующего вида [9.21]:
$$\begin{gather*} \dot {y}_1 = 1 - y_1 y_2, \\ \dot {y}_2 = \alpha y_2 \left({y_1 - \frac{{1 + \beta }}{{y_2 + \beta }}}\right), \end{gather*}$$предложенными Дж. Хиггинсом. В системе $$\beta = 10, \alpha = 100,$$ 200, 400, 1000. Начальные условия
для системы: y1(0) = 1, y2(0) = 0, 001, Tk = 50. Решение этой системы — релаксационные автоколебания (жесткий предельный цикл).
Один из первых и самых популярных примеров
Система Робертсона имеет вид [9.9]
$$\begin{gather*} \dot {y}_1 = - 0.04y_1 + 10^4 y_2 y_3, \\ \dot {y}_2 = 0.04y_1 - 10^4 y_2 y_3 - 3 \cdot 10^7 y_2^2, \\ \dot {y}_3 = 3 \cdot 10^7 y_2^2 . \end{gather*} $$Начальные условия для системы таковы: y1(0) = 1, y2(0) = 0, y3(0) = 0. Рассматриваются следующие величины отрезка интегрирования: Tk = 40 (в работе Робертсона рассматривался именно такой отрезок интегрирования), Tk = 100, 1000, ..., 1011. О свойствах
задачи см. в [9.9].
Данный пример из [9.9] — типичный случай биохимической модели "умеренной" размерности (современные модели, например, фотосинтеза включают сотни уравнений подобного типа). Хотя данная
модель является умеренно жесткой, тем не менее, ее лучше решать с помощью методов, предназначенных для решения
Начальные значения всех переменных системы равны 0, кроме y1(0) = 1 и y8(0) = 0.0057. Длина отрезка интегрирования Tk = 421, 8122.
Еще одна модель химической реакции из [9.9], получившая свое название Е5 в более ранних публикациях.
$$\begin{gather*} \dot {y}_1 = - {Ay}_1 - {By}_1 y_3, \\ \dot {y}_2 = {Ay}_1 - {MCy}_2 y_3, \\ \dot {y}_3 = {Ay}_1 - {By}_1 y_3 - {MCy}_2 y_3 + {Cy}_4, \\ \dot {y}_4 = {By}_1 y_3 - {Cy}_4 . \end{gather*} $$Начальные условия: $$y_1(0) = 1, 76\cdot 10^{- 3},$$ а все остальные
переменные равны 0. Значения коэффициентов модели следующие: $$A = 7, 89 \cdot 10^{- 10}, B = 1, 1\cdot 10^7, C = 1, 13\cdot 10^3, M = 10^6.$$ Первоначально задача ставилась на отрезке Tk = 1000, но впоследствии было обнаружено, что она обладает нетривиальными свойствами вплоть до времени Tk = 1013 (подробнее см. [9.9]).
Обратить особое внимание, что в процессе расчетов приходится иметь дело с очень малыми концентрациями реагентов (малы значения y2, y3 и y4 ). Как "подправить" постановку задачи E5?
Уравнение Релея во многом похоже на уравнение Ван - дер - Поля [9.21]. Рассматривается задача вида
$$$ \ddot {x} - \mu (1 - \dot {x}^2 )\dot {x} + x = 0 . $$$Решить задачу, записав уравнение Релея в виде системы ОДУ. Начальные условия: $$$ x(0) = 0, \dot {x}(0) = 0, 001 $,$$ $$\mu = 1000,$$ Tk = 1000.
Рассмотрим пример системы уравнений, которая описывает изменения численности популяций двух видов и эволюцию некого генетического признака $$\alpha.$$ Система ОДУ имеет вид
$$\begin{gather*} \dot {x} = x(1 - 0, 5x - \frac{2}{{7\alpha ^2}}y), \\ \dot {y} = y(2\alpha - 3, 5\alpha ^2 x - 0, 5y), \\ \dot {\alpha} = \varepsilon (2 - 7\alpha x) . \end{gather*}$$Параметры задачи таковы: $$\varepsilon \le 0, 01, 0 \le x_0 \le 3, 0 \le y_0 \le 15,
\alpha_0 = 0,$$ Tk = 1500. Наличие малого параметра в третьем уравнении системы показывает, что генетический признак меняется медленнее, чем численность популяций. Решение системы — релаксационные колебания.
Задача описана в статье [9.22].
Еще один пример
Параметры задачи таковы: $$\varepsilon \le 0, 01, 0 \le x_0 \le 40, 0 \le y_0 \le 40,
\alpha_{10} = 0, \alpha_{20} = 10,$$ Tk = 2000.
Рассмотреть также модификацию предыдущей системы [9.22]:
$$\begin{gather*} \dot {x} = x(2\alpha_1 - 0, 5x - \alpha_1^3 \alpha_2^{- 3}y), \\ \dot {y} = y(2\alpha_2 - \alpha_1^{- 3} \alpha_{2 }^3 x - 0, 5y), \\ \dot {\alpha}_1 = \varepsilon (2 - 3\alpha_1^2 \alpha_2^{- 3}y), \\ \dot {\alpha}_2 = \varepsilon (2 - 3\alpha_1^{- 3} \alpha_2^2x) . \end{gather*} $$Параметры задачи: $$\varepsilon \le 0, 01, 0 \le x_0 \le 40, 0 \le y_0 \le 40,
\alpha_{10} = 0, \alpha_{20} = 10,$$ Tk = 2000.
Для получения официальных документов о завершении программы дополнительного профессионального образования (удостоверения о повышении квалификации, дипломов о профессиональной переподготовке и MBA) необходимо предоставить:
Внимание! Вы можете не заказывать доставку бумажной версии официального документы, а скачать его в электронном виде и распечатать самостоятельно. Информация о выданном документе в течение 1 месяца загружается в Федеральную информационную систему «Федеральный реестр сведений о документах об образовании и (или) о квалификации, документах об обучении» - ФИС ФРДО.
Доступ на новый сайт осуществляется с использованием адреса электронной почты, который был указан вами при регистрации на "старом". Мы постарались перенести все ваши данные с прежнего ресурса, однако не исключена вероятность потери части информации.
При возникновении проблемы со входом, воспользуйтесь функцией сброса пароля
Если вы обнаружите несоответствия, пожалуйста, сообщите нам.