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

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

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

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

Лекция 2. Дифференциальные уравнения на прямой

Напоминание:

  1. Дифференциальное уравнение $\dot x = f(x,t)$
  2. Решение дифференциального уравнения: $\newcommand{\ph}{\varphi}x=\ph(t)$ такая что выполнено тождество $$\dot \ph(t)=f(\ph(t),t)\quad \forall t$$

На протяжении этой лекции неизвестная функция будет принимать значения на вещественной прямой, то есть будет отображением $$x\colon \mathbb R\to \mathbb R$$

Поле направлений: Вместе с каждым дифференциальным уравнением связано поле направлений — картинка, на которой через каждую точку проведена прямая, угловой коэффициент которой равен значению функции $f(x,t)$. Понятно, что такие прямые — это всеовозможные касательные к решениям уравнения.

Приведем пример для дифференциального уравнения $\dot x = t$:

In [6]:
axes4x4()
rcParams['figure.figsize']=(8,6)
normdirfield(arange(-4,4,0.5),arange(-4,4,0.5),lambda t,x: t,color='red',linewidth=1,length=0.6)

Даже просто посмотрев на этот график, то видно, что решением этого уравнение должна быть парабола (либо функция, очень похожая на параболу). Синим изображена интегральная кривая этого дифференциального уравнения.

In [3]:
axes4x4()
normdirfield(arange(-4,4,0.5),arange(-4,4,0.5),lambda t,x: t,color='red',linewidth=1,length=0.4)
mplot(linspace(-4,4),lambda x: x**2/2, color='blue')

Напомним также, как выглядит поле направлений и интегральные кривые для уравнения $\dot x=x$:

In [4]:
axes4x4()
normdirfield(arange(-4,4,0.5),arange(-4,4,0.5),lambda t,x: x,color='red',linewidth=1,length=0.6)
mplot(linspace(-4,4),lambda t: exp(t), color='blue')

Численное решение дифференциальных уравнений. Метод Эйлера

Пусть поставлена задача Коши: $$\dot x=f(x,t),\quad x(t_0)=x_0$$

Можем приблизительно решать её таким образом. Возьмём точку $(x_0,t_0)$. Решение, проходящее через эту точку, имеет в ней касательную с угловым коэффициентом $f(x_0,t_0)$. Касательная — это прямая, которая хорошо приближает график функции. Давайте на секундочку представим, что график нашего решения в точности совпадает с касательной на некотором небольшом промежутке, содержащем точку $t_0$. Отступим от точки $t_0$ вправо на некоторую маленькую величину $\Delta t$ и рассмотрим точку $(x_1, t_0+\Delta t)$, построенную следующим образом:

$$x_1=x_0+f(x_0,t_0)\cdot \Delta t$$ $$t_1=t_0+\Delta t$$

Она лежит на касательной, проходящей через точку $(x_0, t_0)$. Если $\Delta t$ мало, эта точка должна лежать близко к графику настоящего решения. Затем мы можем взять за стартовую точку $(x_1, t_1)$, построить в ней уже новую касательную, и пройти по ней ещё на $\Delta t$ вправо. Действуя таким образом, получим набор точек, связанных соотношением:

$$x_{n+1}=x_n+f(x_n,t_0)+\Delta t$$ $$t_n=t_{n-1}+\Delta t=t_0+n\Delta t$$

Если соединить эти точки отрезками прямых, они будут проходить вблизи касательных к решению, и сама получающаяся ломаная будет приближаться к настоящему решению. Естественно, с уменьшением шага точность приближения увеличивается.

На графике ниже синим изображено истинное решение уравнения $\dot x = t, \quad x(-3)=4$, а красным, розовым, фиолетовым и зеленым изображены численные решения уравнения методом Эйлера 5, 10, 20 и 100 шагами соответственно.

In [5]:
#import math
def eulersplot(f, xa, xb, ya, n, **kw):
    """plots numerical solution y'=f
    # xa: initial value of independent variable
    # xb: final value of independent variable
    # ya: initial value of dependent variable
    # n : number of steps (higher the better)
    """
    h = (xb - xa) / float(n)
    x = [xa] 
    y = [ya]
    for i in range(1,n+1):
        y.append(y[-1] + h * f(x[-1], y[-1]))
        x.append(x[-1] + h)
    plot(x,y, **kw)

axes4x4()
eulersplot(lambda t,x: t, -3, 4, 4, 5, color='red')
eulersplot(lambda t,x: t, -3, 4, 4, 10, color='pink')
eulersplot(lambda t,x: t, -3, 4, 4, 20, color='purple')
eulersplot(lambda t,x: t, -3, 4, 4, 100, color='green')
mplot(linspace(-4,4),lambda x: x**2/2-0.5, color='blue')

Заметим, что уже сто шагов дает достаточно хорошее приближение решения.

Задача

Найти число $e$, зная, что функция $y=e^x$ удовлетворяет дифференциальному уравнению $y'=y$ и $y(0)=1$.

Аналитическое решение автономных дифференциальных уравнений на прямой

Рассмотрим задачу Коши для автономного дифференциального уравнения

\begin{equation}\tag{*} \dot x=f(x),\quad x(t_0)=x_0 \end{equation}

Пусть $x=\phi(t)$ — решение этой задачи, и $(x_1, t_1)$ — некоторая точка, лежащая на графике. В этом случае $x_1=\phi(t_1)$. Пусть $f(x_1)\ne 0$.

Рассмотрим функцию $t=\psi(x)$, обратную к функции $x=\phi(t)$. Тогда по теореме о производной сложной функции,

$$\psi'(x_1)=\frac{d\psi(x_1)}{dx}=\frac{1}{\dot \phi(t_1)}=\frac{1}{f(x_1)}$$

Это равенство выполняется в любой точке $x_1$. Значит, функция $\psi(x_1)$ является решением дифференциального уравнения

$$\psi'=\frac{1}{f(x)},$$

где $x$ выступает в роли независимой переменной. Такое уравнение мы умеем решать:

$$\psi(x)=\int_{x_0}^x \frac{dx}{f(x)}+t_0$$

Вспоминая, что $t=\psi(x)$ — обратная функция к решению $x=\phi(t)$, имеем:

$$t-t_0=\int_{x_0}^{\phi(t)} \frac{dx}{f(x)}$$

Эта формула называется формулой Барроу. Её можно понимать так: если $f(x)$ — это скорость движения в точке $x$, то время, необходимое для попадания из точки $x_0$ в точку $x$, вычисляется как интеграл от величины, обратной скорости. Кажется, довольно убедительно.

Это соотношение можно понимать как неявное выражение функции $\phi(t)$, то есть решения.

Часто используют такую символическую запись:

$\begin{align} \dot x&=f(x)\\ \frac{dx}{dt}&=f(x)\\ \frac{dx}{f(x)}&=dt\\ \int_{x_0}^x \frac{dx}{f(x)}&=\int_{t_0}^t dt \end{align}$

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

Пример

Решим уравнение $\dot x=x$.

$$x = \varphi(t), \quad x(t_0)=x_0$$

$\varphi^{-1}(x)$ — обратная

$\begin{align} &(\varphi^{-1}(x))'=\frac{1}{x}\\ &\varphi^{-1}(x) = \psi (x)\\ &\psi' (x) = \frac 1 x\\ &\psi (x) = \int \frac{dx} x = \ln |x| + C\\ &t=\ln|x|+C\\ &x = \pm e^{-C}e^t = C_1 e^t \end{align}$

Заметим, что если бы мы забыли модуль под логарифмом при интегрировании, то константа $C_1=e^{-C}$ принимала бы только положительные значения. Но из-за модуля она может принимать и отрицательные значения. Заметим также, что в ходе преобразований (деления на $x$) мы «потеряли» решение $x=0$. Если в ответ подставить значение $C_1=0$, получим как раз его. Таким образом, формула $x(t)=Ce^t$, $c\in \mathbb R$ даёт все известные нам решениям. (Мы пока не доказали, что других нет — на следующей лекции мы обсудим этот вопрос.

In []: