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

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

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

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

Лекция 1. Примеры моделей, приводящих к дифференциальным уравнениям

Рост населения. Мальтузианская модель

Пусть скорость роста популяции какого-нибудь вида (например, рыб в пруду или бактерий в чашке Петри) в любой момент времени пропорциональна количеству особей в популяции в этот момент времени. Это предположение кажется разумным (какая-то часть популяции за единицу времени воспроизводится), если есть достаточное количество ресурсов. Обозначим размер популяции в момент времени $t$ через $x(t)$. Тогда мгновенная скорость роста равна $\frac{dx(t)}{dt}$. Обычно производная по переменной $t$ обозначается точкой $\dot x(t)$, а не штрихом. Таким образом, наш закон роста размера популяции можно записать таким образом: \begin{equation} \dot x(t)=kx(t), \end{equation} где $k>0$ — коэффициент пропорциональности (константа).

Зависимость от $t$ обычно опускают, и пишут просто \begin{equation} \dot x=kx. \end{equation}

Это — одно из простейших (и важнейших) дифференциальных уравнений. Неизвестной величиной в ней является не число (как в обычных алгебраических уравнениях) и не вектор (как в линейной алгебре), а функция $x(t)$.

Рост экономики. Модель Солоу

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

$$\dot k=sf(k)-\delta k,$$

где $k=k(t)$ — капиталовооруженность экономики в момент времени $t$, $s$ — норма сбережения, $\delta$ — норма выбытия капитала.

Механическая система. Падающий шарик

Если я возьму в руку маленький тяжелый шарик, что с ним произойдёт, когда я его отпущу? Не нужно проводить этот эксперимент на практике и даже решать дифференциальное уравение, чтобы ответить: он станет падать вниз с ускорением. Это подскажет нам наша физическая интуиция. Использование интуиции и ранее накопленного опыта очень важно при решении задач, поэтому мы время от времени будем обращаться к механическим примерам.

Пусть вертикальная координата шарика (высота) в момент времени $t$ есть $y(t)$. Известно, что на тело, находящееся в поле тяготения земли (на не слишком большой высоте) действует сила тяжести, равная

$$F=-mg,$$

где $m$ — масса тела, $g$ — ускорение свободного падения (примерно равно 10 м/с^2), знак «-» выбран, поскольку сила тяжести действует в направлении «вниз» (против направления роста $y$).

С другой стороны, второй закон Ньютона гласит, что ускорение тела пропорционально действующей на него силе и обратно пропорционально массе:

$$a=F/m\quad\Leftrightarrow\quad F=ma$$

Ускорение — это вторая производная от координаты по времени, она обозначается двумя точками. Таким образом, для шарика, имеем дифференциальное уравнение:

$$\ddot y=-g$$

Простейшие дифференциальные уравнения

Дифференциальное уравнение общего вида

Дифференциальным уравнением называется соотношение вида $\dot x=f(x,t),$ где $x=x(t)$ — неизвестная функция, $f(x,t)$ — известная функция двух переменных. Мы пока что будем рассматривать уравнения, в которых областью значений неизвестной функции являются вещественные числа $\mathbb R$, но чуть позже обсудим и более сложные случаи, когда $x$ принимает значение в многомерных пространствах.

Решением дифференциального уравнения называется дифференцируемая функция $\newcommand{\ph}{\varphi}x=\ph(t)$, такая, что при подстановке её в уравнение получается верное равенство:

$$\dot \ph(t)=f(\ph(t),t)\quad \forall t$$

Нулевая правая часть

Простейшее дифференциальное уравнение, которое только можно придумать, имеет вид

$$\dot x=0$$

Его решениями являются функции $x(t)=C$, где $C$ — любая константа. Заметим, что даже в таком простейшем случае мы имеем не одно, а сразу целое семейство решений. Аналогичная ситуация будет и в более сложных примерах.

Постоянная правая часть

Чуть более сложное уравнение:

$$\dot x=k,$$

где $k$ — константа. Это уравнение движения с постоянной скоростью, его решениями являются всевозможные функции

$$x(t)=kt+C,$$

Заметим, что в этом случае $C$ задаёт значение функции в начальный момент времени $t=0$.

Правая часть, зависящая только от времени

$$\dot x=f(t)$$

Задачу отыскания решения такого дифференциального уравнения можно сформулировать следующим образом: для каждого значения независимой переменной $t$ известна проиводная некоторой функции; найти эту функцию. Нетрудно видеть, что это в точности задача интегрирования (отыскания первообразной). Решение такого уравнения задаётся таким образом неопределенным интегралом, который можно записать в виде

$$x(t)=\int f(t) dt=\int_{t_0}^t f(\tau)d\tau+C$$

(Неопределенный интеграл по определению является семейством функций, а при записи его в виде определенного интеграла с переменным верхним пределом нужно указывать константу интегрирования явным образом.)

Начальные условия. Задача Коши

Чтобы выделить среди семейства решений дифференциального уравнения одно, обычно вместе с самим дифференциальным уравнением рассматривают дополнительное соотношение, называемое начальным условием — значение решение в какой-то момент времени (не обязательно $t=0$). Например, можно рассмотреть такую задачу:

$$\dot x=f(t),\quad x(5)=0$$

В этом случае решением такого уравнения будет уже только одна функция:

\begin{equation}\tag{1} x(t)=\int_5^t f(\tau)d\tau \end{equation} Когда задано дифференциальное уравнение и начальное условие, говорят, что поставлена задача Коши.

Заметим, что функция, заданная в (1), является решением данного уравнения. При подстановке $t=5$ решение $x(5)=\int_5^5 f(\tau)d\tau = 0$, а также, продифференцировав его, имеем: $\dot x(t)=\frac{d}{dt}\int_5^t f(\tau)d\tau = f(t)$

Контрольный вопрос: какое будет решение уравнения при начальном условии $x(5)=1$?

Простейшее линейное уравнение

Положим в уравнении роста населения $k=1$. Получим следующее уравнение:

$$\dot x=x$$

Какие функции будут его решениями? Словами можно сказать, что условие, накладываемое этим уравнением, звучит так: «Производная функции равна самой этой функции». Одна известная функция обладает таким свойством — это экспонента $x(t)=e^t$. Нетрудно видеть, что если домножить экспоненту на любое число, получающаяся функция $x(t)=Ce^t$ также будет решением этого уравнения. В частности, очевидно, что решением будет функция $x(t)\equiv 0$.

Контрольный вопрос: является ли решением этого уравнения функция $x(t)=e^t+C$?

Замечание. Мы пока не доказали, что других решений, кроме перечисленных, у этого уравнения нет. Мы вернемся к этому вопросу через одну лекцию.

Геометрические объекты

В рассмотренных выше примерах неизвстная функция $x(t)$ принимала значения во множестве вещественных чисел. В общем случае функция $x(t)$ может принимать значения в других множествах — например, в многомерных пространствах. Множество, в котором принимает значение неизвестная функция (или, иными словами, множество всевозможных значений $x(t)$ при каком-нибудь фиксированном $t$) называется фазовым пространством дифференциального уравнения. Множество точек вида $(x,t)$ (декартово произведение фазового пространсва на ось времени) называется расширенным фазовым пространством. График решения называется интегральной кривой. (Интегральные кривые живут в расширенном фазовом пространстве.)

Построим некоторые интегральные кривые для уравнения $\dot x=x$. Как мы уже знаем, ими будут графики экспонент.

In [2]:
rcParams['figure.figsize']=(10,8)
axes4x4()
for C in range(-5,6,2)+[0,0.2,-0.2]:
    mplot(linspace(-4,4),lambda t: C*exp(t),color='blue') 

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

Вот что это такое. Возьмём произвольную точку $P=(x_0,t_0)$ расширенного фазового пространства. Например, $x_0=3$, $t_0=1$. Мы можем провести в точке $P$ касательную к интегральной кривой, проходящей через эту точку. Действительно, чтобы провести прямую через фиксированную точку, нужно знать только её угловой коэффициент, но угловой коэффициент касательной к графику некоторой функции равняется производной этой функции. А производную решения мы знаем, по определению решения она равна правой части уравнения. В нашем случае, $\dot x=x$ и значит касательная, проходящая через точку $P$, имеет угловой коэффициент, равный $x_0=3$. Можно взять ещё несколько точек на прямой $t=1$ и провести соответствующие касательные через них. Получится такая картинка.

In [3]:
axes4x4()
xs=range(-4,5)
scatter([1]*len(xs),xs)
for x0 in xs:
    mplot(linspace(0,5),lambda t: x0*(t-1)+x0,linestyle='dashed',linewidth=0.5,color='black')
    mplot(linspace(1-0.2,1+0.2),lambda t: x0*(t-1)+x0,color='red',linewidth=1.5)

Контрольный вопрос: почему прямые пересекаются в начале координат?

Понятно, что можно, действуя аналогично, построить касательные к решениям не только в выбранных точках, но и вообще в любой точке расширенного фазового пространства. В данном случае правая часть не зависит от $t$ явно, поэтому через любые две точки, лежащие на одной горизонтальной прямой, будут проходить параллельные касательные. Мы будем рисовать только маленькие кусочки этих касательных.

In [4]:
axes4x4()
normdirfield(arange(-4,4,0.5),arange(-4,4,0.5),lambda x,y: y,color='red',linewidth=1,length=0.6)

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

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

In [5]:
axes4x4()
normdirfield(arange(-4,4,0.5),arange(-4,4,0.5),lambda x,y: y,color='red',linewidth=1,length=0.6)
for C in range(-5,6,2)+[0,0.2,-0.2]:
    mplot(linspace(-4,4),lambda t: C*exp(t),color='blue')