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\).

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

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

Пусть \(x=\phi(t)\) — решение этой задачи, и \((x_1, t_1)\) — некоторая точка, лежащая на графике. В этом случае \(x_1=\ohi(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)\), то есть решения.

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

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

Пример

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

\[x = \varphi(t), \quad x(t_0)=x_0\]

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

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