In [2]:
%pylab inline
from nesode import *
Populating the interactive namespace from numpy and matplotlib

Дифференциальные уравнения

Совместный бакалавриат ВШЭ-РЭШ, 2013-14 учебный год

И. В. Щуров, П. Ф. Соломатин, И. А. Хованская, А. Петрин, Н. Солодовников

Предлагаемое решение домашнего задания 4

благодарим Михаила Заварзина за предоставленный конспект

Задача 1

Найти уравнение фазовой кривой $$ \begin{cases} \dot x = x-4y+21\\ \dot y = \frac{1}{ x-4y+21} \end{cases} $$ удовлетворяющей начальному условию $x(0)=1$, $y(0)=5$.

Решение

Перейдём от системы автономных линейных дифференциальных уравнений к неавтономному уравнению, поделив $\dot x$ на $\dot y$: $\frac{dx}{dy} = \left(x - 4\,y + 21\right)^2,\,x-4\,y+21\neq0$.

Далее вводим замену $z = x - 4\,y + 21;\,z'_y = x'_y - 4 = (x - 4\,y + 21)^2 - 4$.

Замечаем, что при начальном условии $x(0) = 1,\,y(0) = 5$ переменная $z$ принимает знаечение $z = 1 - 4\cdot5 + 21 = 2$, а тогда дифференциальное уравнение в новых координатах $z,y$ принимает вид $z'_y = 2^2 - 4 = 0$, а значит, прямая $z = 2$ (равно как и $z = -2$) является решением неавтономного дифференциального уравнения и прямой, на которой лежит фазовая кривая (в новых координатах) системы автономных уравнений (стационарное решение в новых координатах: $z'_y = 2'_y = 0$). В старых координатах при данных начальных условиях фазовой кривой системы будет прямая $z = x - 4\,y + 21 = 2;\,y = \frac{x+19}{4}$.

Задача 2а

Рассмотрим систему уравнений

$$ \begin{cases} \dot x = (x-3)^2\\ \dot y = (y-2)^2 \end{cases} $$

Записать уравнение всех фазовых кривых в виде зависимости $y$ от $x$ или $x$ от $y$. Допускается неявное задание кривых с указанием области определения.

Решение

$$ \frac{dx}{dy} = \frac{(x-3)^2}{(y-2)^2} $$ $$ \int \frac{1}{(x-3)^2} dx = \int \frac{1}{(y-2)^2} dy $$ $$ \frac{1}{x-3}-\frac{1}{y-2}=C $$ Заметим, что здесь не зватает постоянных решений: $x=3$ и $y=2$. На них находятся три фазовых кривых (включая саму точку $(3,2)$). На общем виде решений при каждом $C$ находятся две интегральные кривые. Область определения, соответственно, ограничена прямыми $x=3$ и $y=2$.

Задача 2b

Найти все начальные условия $(x_0, y_0)$ для которых решение $(x(t),y(t))$ соответствующей задачи Коши определено для сколь угодно больших значений $t$ и стремится к точке $(3,2)$, то есть $$\lim_{t\to\infty} (x(t),y(t)) = (3,2)$$

Решение

Для данной системы дифференциальных уравнений можно найти явно найти зависимости $x(t), y(t)$.

$\dot x = (x - 3)^2;\,\int\frac{dx}{(x-3)^2} = \int\,dt;\,-\frac{1}{x-3} = t + c_1; x = - \frac{1}{t + c_1} + 3$.

Отсюда видно, что интегральные кривые этого уравнения — гиперболы, которые в зависимости от начальных условий «передвигаются» вдоль оси $t$. Причём, если $x(0)$ находится ниже 3, то в прямом времени мы будем стремиться к этому значению, а если $x(0) > 3$, то решение в прямом времени, подходя к $t = -c_1$, будет уходить на бесконечность. Не забудем про стационарное решение данного уравнения: $x = 3, \dot 3 = 0$, которое так же в прямом времени стремится к 3 (собственно, «к самому себе»).

Аналогично $\dot y = (y - 2)^2;\,\int\frac{dx}{(y-2)^2} = \int\,dt;\,-\frac{1}{y-2} = t + c_2; y = - \frac{1}{t + c_2} + 2$.

Отсюда видно, что интегральные кривые этого уравнения — гиперболы, которые в зависимости от начальных условий «передвигаются» вдоль оси $t$. Причём, если $y(0)$ находится ниже 2, то в прямом времени мы будем стремиться к этому значению, а если $x(0) > 2$, то решение в прямом времени, подходя к $t = -c_2$, будет уходить на бесконечность. Не забудем про стационарное решение данного уравнения: $y = 2, \dot 2 = 0$, которое так же в прямом времени стремится к 2 (собственно, «к самому себе»).

В итоге получается, что стремление к точке $(3,2)$ решения соответствующей задачи Коши будет достигаться только при $\left\{\begin{aligned}x(0) &\leqslant 3\\ y(0) &\leqslant 2\\ \end{aligned} \right.$.

В этом также можно убедиться, нарисовав векторное поле для данной системы уравнений:

In [3]:
axis([-1,6,-1,6])
rcParams["figure.figsize"]=(8,8)
plot([-1,3],[2,2],[3,3],[-1,2],[3,6],[2,2],'--k',[3,3],[2,6],'--k',3,2,'bo',color='blue',linewidth=2)
vectorfield(arange(-1,6,0.5),arange(-1,6,0.5),lambda x,y:(float(1)*(x-3)**2,float(1)*(y-2)**2))
eulersplot(lambda x,y:float(1)*(y-2)**2/((x-3)**2), -1, 2.96, -1, n=3000, color='turquoise',linewidth=2)
eulersplot(lambda x,y:float(1)*(y-2)**2/((x-3)**2), -1, 2.96, 1, n=3000, color='turquoise',linewidth=2)
eulersplot(lambda x,y:float(1)*(y-2)**2/((x-3)**2), -1, 2.5, 2.5, n=3000, color='orange',linewidth=2)
eulersplot(lambda x,y:float(1)*(y-2)**2/((x-3)**2), 3.05, 6, 2.05, n=3000, color='orange',linewidth=2)
eulersplot(lambda x,y:float(1)*(y-2)**2/((x-3)**2), 3.5, 6, -1, n=3000, color='orange',linewidth=2)

Задача 3

Найди такую функцию $f(x,y)$, чтобы функция $$H(x,y)=e^{-3x^3y^3-4x^2y^2+4}$$ была первым интегралом системы \begin{align} \dot x = &f(x,y) & \dot y = & -9x^3y^6-8x^2y^5 \end{align}

Решение

Сначала упростим себе немного жизнь, изменив вид первого интеграла. Как мы обсуждали на лекции, любое монотонное преобразование первого интеграла также является первым интегралом. Предлагается сначала прологарифмировать (убрать экспонту), а потом вычесть четверку. Новый первый интеграл $$\overline H(x,y)={-3x^3y^3-4x^2y^2}$$ Теперь продифференцируем функцию $H$ вдоль вектроного поля и приравняем её к нулю. $$\begin{multline} L_v H = H'_x \dot x+H'_y \dot y = (-9x^2y^3-8xy^2)f(x,y)+(-9x^3y^6-8x^2y^5)(-9x^3y^2-8x^2y)=\\ =(-9x^2y^3-8xy^2)(f(x,y)+(xy^4)(-9x^3y^2-8x^2y))=0 \end{multline} $$ $$ f(x,y)=9x^4y^6+8x^3y^5 $$

Задача 4

Найдите все значения параметра $\alpha$, при которых у следующей системы существует хотя бы один непрерывный глобально определенный непостоянный первый интеграл. Найте этот первый интеграл для всех значений $\alpha$, при которых он существует.

Решение

Сразу отметим триваальные случаи, когда при $\alpha = 5~H(x,y,z) = x$ (непрерывный и глобальный), т.к. $\langle (1,0,0), (\dot x, \dot y, \dot z)\rangle=0$. Аналогично при $\alpha = 2~H(x,y,z) = y$, a при $\alpha = 3~H(x,y,z) = z$. В этом легко убедиться, взяв производную по направлению.

Далее, исходя из того, что мы не хотим, чтобы все (кроме, тех, что сами по себе лежат на линиях уровня первого интеграла со значением 0) фазовые кривые сходились в особую точку $(0,0,0)$ (других особых точек нет) в прямом или обратном времени (ведь тогда глобального непрерывного первого интеграла не существует, т.к. фазовые кривые, соответствующие произвольным значениям первого интеграла, будут сходиться в точку с каким-то конкретным значением $H(x,y,z)$ — противоречие), можно сказать, что нам подходят все значения $\alpha \in \mathbb{R}$, кроме тех, что задают коэффициенты с одинаковыми знаками перед $x,y,$ и $z$.

$$1.~\alpha < 2: (-4\alpha + 20) > 0;\,(-3\alpha + 6) > 0;\,(-2\alpha + 6) > 0$$ $$2.~2 < \alpha < 3: (-4\alpha + 20) > 0;\,(-3\alpha + 6) < 0;\,(-2\alpha + 6) > 0$$ $$3.~3 < \alpha < 5: (-4\alpha + 20) > 0;\,(-3\alpha + 6) < 0;\,(-2\alpha + 6) < 0$$ $$4.~\alpha > 5: (-4\alpha + 20) < 0;\,(-3\alpha + 6) < 0;\,(-2\alpha + 6) < 0$$

Таким образом, нам подходят $\alpha \in [2,5]$.

Теперь попробуем найти такой $H(x,y,z)$ для промежутков $2.,3.$ Заметим, что $\forall \alpha \in (3,5)$ (включает $2.,3.$) коэффициенты, стоящие перед $x$ и $y$, имеют разные знаки. Это значит, что на ненулевых «поверхностях» уровня соответствующего первого интеграла будут лежать фазовые кривые, не попадающие в особую точку с нулевым значением $H(x,y,z)$ (ведь движение к $(0,0,0)$ подразумевает то, что все переменные имеют «одинаково направленное» направление изменения абсолютного значения). Найдём этот первый интеграл, пользуясь тем, что нам уже не нужно знать его поведение относительно $z$, т.к. знаки перед $x$ и $y$ уже различны:

$$\frac{dy}{dx} = \frac{(-3\alpha+6)\,y}{(-4\alpha+20)\,x};\,\int\frac{(-4\alpha+20)\,dy}{y} = \int\frac{(-3\alpha+6)\,dx}{x};\,(-4\alpha+20)\ln|y| = (-3\alpha+6)\ln|x| + c;$$

$$H(x,y,z) = |y|^{-4\alpha+20}|x|^{3\alpha-6}, \forall z, \left\{\begin{aligned}\alpha &\leq 5\\ \alpha &\geq 2\\ \end{aligned} \right. .$$

Задача 5а

Рассмотрим уравнение $$\ddot x = (-\alpha+x)(x-3)(x+2)$$ Найдите первый интеграл системы при $\alpha=0$

Решение

Введем стандартную замену $y=\dot x$.

Теперь найдем формулу потенциальной энергии: $$U(x)=-\int F \, dx = -\int x(x-3)(x+2) \, dx = - x^4/4+x^3/3+x^2$$ Константу интегрирования приравняли у нулю, чтобы в нуле значение потенциальной энергии было равно нулю. Кинетическая энергия равна $V(y)=y^2/2$, и находим полный интеграл системы (уравнение полной энергии): $$ H(x,y) = E(x,y)=U(x)+V(y)=y^2/2 - x^4/4+x^3/3+x^2 $$

Задача 5b

Построить фазовый портрет уравнения в координатах $(x,\dot x)$ при $\alpha =0 $. Отметить на нём:

  1. Все особые точки
  2. Как минимум одну траекторию, соответствующую непостоянному периодическому решению
  3. Все траектории, стремящиеся к какой-либо особой точке в прямом или обратном времени
  4. Как минимум одну траекторию, соответствующую решению, которое не является постоянным, периодическим и не имеет конечного предела ни в прямом, ни в обратном времени.

Решение

In [5]:
rcParams['figure.figsize']=(6,6)
axis([-5,5,-16,16])
text(4.4,0.15,"$x$",fontsize=20)
text(0.1,15.1,"$y = \dot x$",fontsize=20)
center_spines()
f=lambda x,y:-float(1)*(x**4)/4 + float(1)*(x**3)/3 + float(3)*x**2 + y**2/2
mcontour(linspace(-5,2,1000),linspace(-16,16,1000),f, [16.0/3], colors='blue' )
mcontour(linspace(-5,5,1000),linspace(-16,16,1000),f, [15.75],  colors='blue' )
mcontour(linspace(-2,2,1000),linspace(-16,16,1000),f, [3.0],    colors='red'  )
mcontour(linspace(-5,5,1000),linspace(0,16,1000),  f, [30],     colors='black')
normvectorfield(arange(-5,5,1),arange(-15,16,1),lambda x,y:(y,float(1)*x*(x-3)*(x+2)), color='turquoise')

Особые точки $(-2,0)$, $(0,0)$ и $(3,0)$.

Синии траектории (их девять) соответствуют траекториям, стремящихся к особым точкам в прямом или обратном времени.

Красная кривая соотвествует периодическому решению.

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

Задача 5c.

$$\max U = 15.75, x = 3;\,\max U = \frac{16}{3}, x = -2;\,\min U = 0, x = 0$$ Исходя из того, что мы знаем о функции полной энергии для консервативной системы с одной степенью свободы, можно сделать вывод, что точки $\left(-2,\frac{16}{3}\right), (3, 15.75)$ будут особыми точками типа «седло» , а точка $(0,0)$ — особой точкой типа «центр». Cм. левую часть графика ниже.

Различные конечные пределы у траекторий соответствуют ситуации, когда они выходят из одной особой точки (в данном случае точки типа «седло») и стремятся к другой в прямом и обратном времени. Такой ситуации при $\alpha = 0$ быть не может, ведь максимумы функции потенциальной энергии находятся на разных уровнях ($15.75$ и $\frac{16}{3}$), и поэтому линии уровня, проходящие через одну из двух таких особых точек, не могут содержать никакую точку, лежащую на более высоком или низком «уровне» (ведь $H(x,y)$ — непрерывная функция и она не может при одних и тех же $(x, y)$ принимать разные значения соответствующие разным линиям уровня). А значит, и фазовые кривые, лежащие на линиях уровня $H(x,y)$, выйдя из одной особой точки, не смогут устремиться к другой.

Поэтому, чтобы такое могло произойти, небходимо, чтобы максимумы функции были на одном уровне, тогда через них сможет пройти одна линия уровня, на которой будут лежать фазовые кривые, стремящиеся к одной особой точке в прямом врмени и к другой — в обратном. Такая линия уровня точно будет существовать в силу непрервыности $H(x,y)$ для всех $x, y \in \mathbb{R}$.

Такая ситуация будет возможна только, если поместить корень производной $U(x)$, cоответсвующий точке минимума этой функции, посередине между двумя точками максимума $U(x)$. Это можно объяснить формой функции потенциальной энергии. Функция обязана быть многочленом четвертой степени, при этом максимумы должны достигаться на одном и том же уровне. Это возможно только тогда когда функция симметрична, а значит минимум будет посередине. Более формально это можно доказать найдя уравнение потенциальной энергии (в зависимости от $x$ и $\alpha$ ), и приравнять значение этой функции $x=-2$, $x=3$.

Таким образом, данной ситуации соответствует $\alpha = \frac{3 + (-2)}{2} = 0.5$

Тогда $U(x) = -\int(x - 0.5)(x+2)(x-3) \, dx = -3x+2.75\,x^2+0.5\,x^3-0.25\,x^4$

$H(x,y) = E(x,y) = U(x) + V(y) = -3x+2.75\,x^2+0.5\,x^3-0.25\,x^4 + \frac{y^2}{2}$. $\max U = 9, x = 3, x = -2;\,\min U = -\frac{49}{64}, x = 0$. См. графики справа на рисунке выше. Теперь действительно видно, что на линии уровня, соответствующей значению 9, лежат 2 фазовые кривые, имеющие различные конечные пределы в прямом и обратном времени.

In [9]:
rcParams['figure.figsize']=(18,18)

fig = plt.figure()
#FIGURE 1
ax=fig.add_subplot(221) 
plot([-2,-2],[16,-16],'--k',[3,3],[16,-16],'--k',[0,0],[16,-16],'--k',linewidth=2)
axis([-5,5,-16,16])
text(4.4,0.15,"$x$",fontsize=27)
text(0.7,7.8,"$U(x)$",fontsize=27)
center_spines()
mplot(linspace(-5,5,1000), lambda x: -float(1)*(x**4)/4 + float(1)*(x**3)/3 + float(3)*x**2)
#FIGURE 2
ax=fig.add_subplot(222) 
plot([-2,-2],[16,-16],'--k',[3,3],[16,-16],'--k',[-5,5],[9,9],'--k',linewidth=2)
plot([0.5,0.5],[16,-16],'--k',linewidth=2,color='red')
axis([-4,5,-16,16])
text(4.4,0.15,"$x$",fontsize=27)
text(1,7.8,"$U(x)$",fontsize=27)
center_spines()
mplot(linspace(-5,5,1000), lambda x: -float(3)*x+2.75*x**2+0.5*x**3-0.25*x**4)
#FIGURE 3
ax=fig.add_subplot(223) 
plot([-2,-2],[16,-16],'--k',[3,3],[16,-16],'--k',[0,0],[16,-16],'--k',linewidth=2)
axis([-5,5,-16,16])
text(4.4,0.15,"$x$",fontsize=27)
text(0.1,15.1,"$y = \dot x$",fontsize=27)
center_spines()
mcontour(linspace(-5,5,1000),linspace(-16,16,1000),lambda x,y:-float(1)*(x**4)/4 + float(1)*(x**3)/3 + float(3)*x**2 + y**2/2, list(arange(-4,20,2))+
         list(arange(float(16)/3,1+(16)/3,1))+list(arange(15.75,16.75,1)))
normvectorfield(arange(-5,5,1),arange(-15,16,1),lambda x,y:(y,float(1)*x*(x-3)*(x+2)), color='turquoise')
#FIGURE 4
ax=fig.add_subplot(224) 
plot([-2,-2],[16,-16],'--k',[3,3],[16,-16],'--k',linewidth=2)
plot([0.5,0.5],[16,-16],'--k',linewidth=2,color='red')
plot([0.5,0.5],[1.8,-1.8],linewidth=6,color='black')
axis([-4,5,-16,16])
text(4.4,0.15,"$x$",fontsize=27)
text(0.1,15.1,"$y = \dot x$",fontsize=27)
center_spines()
mcontour(linspace(-5,5,1000),linspace(-16,16,1000),lambda x,y:-float(3)*x+2.75*x**2+0.5*x**3-0.25*x**4 + y**2/2, list(arange(-5,20,2)))
mcontour(linspace(-5,5,1000),linspace(-16,16,1000),lambda x,y:-float(3)*x+2.75*x**2+0.5*x**3-0.25*x**4 + y**2/2, list(arange(9,10,1)))
normvectorfield(arange(-5,5,1),arange(-15,16,1),lambda x,y:(y,-x**2 - 5*x - 6),color='turquoise')