%pylab inline
from nesode import *
Совместный бакалавриат ВШЭ-РЭШ, 2013-14 учебный год
И. В. Щуров, П. Ф. Соломатин, И. А. Хованская, А. Петрин, Н. Солодовников
Напоминание:
На протяжении этой лекции неизвестная функция будет принимать значения на вещественной прямой, то есть будет отображением \[x\colon \mathbb R\to \mathbb R\]
Поле направлений: Вместе с каждым дифференциальным уравнением связано поле направлений — картинка, на которой через каждую точку проведена прямая, угловой коэффициент которой равен значению функции \(f(x,t)\). Понятно, что такие прямые — это всеовозможные касательные к решениям уравнения.
Приведем пример для дифференциального уравнения \(\dot x = t\):
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)
Даже просто посмотрев на этот график, то видно, что решением этого уравнение должна быть парабола (либо функция, очень похожая на параболу). Синим изображена интегральная кривая этого дифференциального уравнения.
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\):
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 шагами соответственно.
#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=\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\) даёт все известные нам решениям. (Мы пока не доказали, что других нет — на следующей лекции мы обсудим этот вопрос.