Skip to content

Instantly share code, notes, and snippets.

@mitya57
Last active December 24, 2015 17:59
Show Gist options
  • Save mitya57/6839254 to your computer and use it in GitHub Desktop.
Save mitya57/6839254 to your computer and use it in GitHub Desktop.

Отчёт о проведённом исследовании системы дифференциальных уравнений

Автор:Дмитрий Шачнев
Группа:411
Дата:2013-10-10

Постановка задачи

Дана система дифференциальных уравнений в \mathbb{R}^2:

\left\{
\eqalign{
  \dot{x} &= P(x, y) &= y \\
  \dot{y} &= Q(x, y) &= ky - x - y^3
}
\right.

При значениях параметра k = 1, -1, 3 требуется:

  1. Найти собственные значения и собственные вектора матрицы Якоби.
  2. Найти особую точку и определить её тип.
  3. Найти длину периода, если существует периодическое решение.
  4. Нарисовать фазовый портрет системы.

Исследование системы

Из уравнения y = ky - x - y^3 = 0 видно, что единственная особая точка — это (0, 0). Чтобы найти её тип, надо сначала вычислить собственные значения матрицы Якоби:

\det\left(
\begin{matrix}
 -\lambda & 1 \\
 -1 & k - \lambda
\end{matrix}
\right) = 0
\lambda^2 - k\lambda + 1 = 0
\lambda = \frac{k \pm \sqrt{k^2 - 4}}{2}
  • При k = 1 имеем \lambda = \frac{1 \pm \sqrt{3}i}{2}, численные значения 0.5 + 0.866i и 0.5 - 0.866i, собственные вектора — (1, \frac{1 \pm \sqrt{3}i}{2}).

    Корни комплексные с положительной вещественной частью, то есть особая точка — неустойчивый фокус.

    Периодическое решение существует, длина периода равна 7.7.

  • При k = -1 имеем \lambda = \frac{-1 \pm \sqrt{3}i}{2}, численные значения -0.5 + 0.866i и -0.5 - 0.866i, собственные вектора — (1, \frac{-1 \pm \sqrt{3}i}{2}).

    Корни комплексные с отрицательной вещественной частью, то есть особая точка — устойчивый фокус.

    Периодического решения не существует.

  • При k = 3 имеем \lambda = \frac{3 \pm \sqrt{5}}{2}, численные значения 2.618 и 0.382, собственные вектора — (1, \frac{3 \pm \sqrt{5}}{2}).

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

    Периодическое решение существует, длина периода равна 17.1.

В случаях k = 1, 3 имеет место предельный цикл, на который «накручиваются» траектории.

Граничные значения: k = 0 — центр, k = 2 — нестабильный вырожденный узел.

Алгоритм работы

Построение траектории движения системы

Для построения траектории из точки (x, y) сначала выбирается некоторый луч \mathcal{L}, выходящий из особой точки (например, ось Ox).

Далее, для построения каждого шага траектории используется метод Рунге–Кутты. Ниже описан алгоритм метода Рунге–Кутты четвёртого порядка точности:

\left\{
\eqalign{
  P_1 &= P \left( x + \frac{\tau}{2} \cdot P(x, y),\;
                 y + \frac{\tau}{2} \cdot Q(x, y) \right), \\
  Q_1 &= Q \left( x + \frac{\tau}{2} \cdot P(x, y),\;
                 y + \frac{\tau}{2} \cdot Q(x, y) \right)
}
\right.
\left\{
\eqalign{
  P_2 &= P \left( x + \frac{\tau}{2} \cdot P_1,\;
                 y + \frac{\tau}{2} \cdot Q_1 \right), \\
  Q_2 &= Q \left( x + \frac{\tau}{2} \cdot P_1,\;
                 y + \frac{\tau}{2} \cdot Q_1 \right)
}
\right.
\left\{
\eqalign{
  P_3 &= P ( x + \tau \cdot P_2,\; y + \tau \cdot Q_2 ), \\
  Q_3 &= Q ( x + \tau \cdot P_2,\; y + \tau \cdot Q_2 )
}
\right.
\left\{
\eqalign{
  x' &= x + \frac{\tau}{6} \left( P(x, y) + 2 \cdot P_1
          + 2 \cdot P_2 + P_3 \right), \\
  y' &= y + \frac{\tau}{6} \left( Q(x, y) + 2 \cdot P_1
          + 2 \cdot P_2 + P_3 \right)
}
\right.

Соединив (x, y) и (x', y'), получим шаг траектории. Такие шаги будем продолжать делать до тех пор, пока не будет выполнено одно из трёх условий конца цикла:

  1. Траектория пересечёт некоторый луч два раза в точках, отличающихся не более чем на \frac{\tau}{2} (циклическое решение).
  2. Тректория подойдёт к особой точке на расстояние \frac{\tau}{2} (стремление к особой точке).
  3. Траектория отойдёт от начала координат на расстояние больше \frac{10}{\tau}, либо будет достигнуто максимальное число шагов \frac{50}{\tau} (стремление к бесконечности).

Поиск неподвижной точки отображения Пуанкаре

Пусть задан луч \mathcal{L} и точка (x_0, y_0) = (x_0(\rho), y_0(\rho)) \in \mathcal{L} (где \rho — локальная координата на луче \mathcal{L}).

Пустим из этой точки траекторию движения системы в положительном времени, и пусть (x_1, y_1) = (x_1(\chi), y_1(\chi)) — первое пересечение этой траектории с лучом \mathcal{L}. Отображение \rho \to \chi = \chi(\rho) называется отображением Пуанкаре этого луча.

Обозначим \psi(\rho) = \chi(\rho) - \rho. Для нахождения неподвижной точки отображения \chi (то есть нуля \psi) используется упрощённый метод Ньютона:

  • Выбираются две точки a и b, такие, что \psi(a) \cdot \psi(b) < 0.

  • Пока расстояние \left|a - b\right| не достигнет необходимой минимальной величины, выполняется замена

    \eqalign{
      a &\to b - (b - a) \cdot \frac{\psi(b)}{\psi(b) - \psi(a)}, \\
      b &\to a - (a - b) \cdot \frac{\psi(a)}{\psi(a) - \psi(b)}
    }
    

Поиск периода

Для поиска периода из неподвижной точки отображения Пуанкаре выпускается траектория с шагом \tau. Длиной периода будет N \cdot \tau, где N — число шагов, за которое алгоритм Рунге–Кутты вернётся в исходную точку.

Определение \varepsilon — погрешности метода

Определим погрешность метода Рунге–Кутты с шагом \tau, запущенного из начальной точки (x_0, y_0). Для этого из той же точки запустим алгоритм с шагом \frac{\tau}{2}.

Пусть (x_i, y_i) и (\tilde{x}_i, \tilde{y}_i) — точки с номером i первой и второй траекторий соответственно. Тогда погрешность метода вычисляется по формуле

\varepsilon = \max_{0 \leqslant i \leqslant N}
\sqrt{ \left( x_i - \tilde{x}_{2i} \right)^2 +
       \left( y_i - \tilde{y}_{2i} \right)^2 }.

Определение \delta — погрешности поиска неподвижной точки

Неподвижная точка существует, если существует решение, пересекающее некоторый луч два раза в точках, остоящих друг от друга не более, чем на \frac{\tau}{2}. Тогда расстояние между этими точками и будет погрешностью поиска неподвижной точки \delta.

Численные значения погрешностей

k = 1
\tau \varepsilon \delta
0.0001 1.09613e-07 1.76144e-11
0.0002 7.7602e-07 3.24754e-11
0.0005 9.25363e-06 1.72758e-10
0.001 5.59294e-05 1.46514e-10
0.002 0.000129327 8.42441e-10
0.005 0.000938488 5.29843e-06
0.01 0.00334614 3.76117e-05
0.02 0.0130162 0.000242305
0.05 0.0885789 0.000137598
0.1 0.376933 0.00013097
k = -1
\tau \varepsilon \delta
0.0001 1.31033e-10 8.31218e-09
0.0002 3.36056e-10 3.32489e-08
0.0005 3.96062e-10 2.07774e-07
0.001 8.73237e-10 8.30426e-07
0.002 1.72283e-09 3.30989e-06
0.005 2.09922e-09 2.0184e-05
0.01 3.59741e-09 7.44739e-05
0.02 4.21477e-09 0.00021614
0.05 3.99769e-08 0.000112655
0.1 6.56572e-07 0.000268919
k = 3
\tau \varepsilon \delta
0.0001 2.6328e-07 1.73772e-11
0.0002 1.59201e-06 7.64651e-11
0.0005 1.43449e-05 3.25251e-11
0.001 8.42552e-05 4.14602e-12
0.002 0.000578843 1.51516e-09
0.005 0.00145192 5.17844e-11
0.01 0.010968 1.04374e-08
0.02 0.0388027 8.47973e-07
0.05 0.233205 3.07063e-06
0.1 1.00551 9.50295e-05

Реализация

Программа написана на языке C++ с использованием фреймворка Qt.

Система уравнений считывается из файла, в котором хранятся коэффициенты многочленов \dot{x} и \dot{y}.

Реализованы два варанта отрисовки:

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

Основным объектом является система типа PolynomSystem, для которой реализованы следующие методы:

/* Получить следующую точку траектории, с шагом eps */
Point getNextValue(Point point, double eps);

/* Отображение Пуанкаре в точке p, с шагом eps */
double poincareFunction(double p, double eps);

/* Найти неподвижную точку отображения Пуанкаре
 * с начальными точками a и b и шагом eps */
double findPoincareStaticPoint(Point a, Point b, double eps);

/* Найти собственные значения матрицы Якоби */
complex[2] getEigenValues();

/* Определить тип особой точки */
PointType getPointType();
scr0.png

Параметр: k = 1

scr1.png

Параметр: k = 3

scr2.png

Параметр: k = -1

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment