This chapter hasn’t been translated into English yet, so here is the Russian original. Your browser can translate the page; the formulas and widgets work the same. Back to the English contents →

Часть IV · Анализ Глава 31 из 60

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

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

1–2 курс 55 минут

Опирается на: 30 · Ряды

Вы научитесь

читать дифференциальное уравнение как закон изменения и видеть его решения в поле направлений
решать уравнения численно методом Эйлера и понимать, откуда берётся ошибка
решать линейные уравнения второго порядка через характеристическое уравнение и объяснять порог эпидемии через R₀

Линия IV · Анализ

Вы здесь: 31 · Дифференциальные уравнения

Пересадки 59 · Хаос и фракталы

Прошлая глава закончилась наблюдением: законы природы говорят не «сколько», а «как меняется». Чашка кофе остывает тем быстрее, чем она горячее. Маятник ускоряется тем сильнее, чем дальше отклонён. Эпидемия растёт тем быстрее, чем больше больных и чем больше тех, кто ещё не болел. В каждом таком законе неизвестная функция стоит рядом со своей производной, и просто проинтегрировать скорость, как в главе 28, не выходит: скорость зависит от того, что мы ищем.

Эта глава — лаборатория. У каждого стенда своя модель и свои ручки: остывание и рост, пружина, волки и зайцы, эпидемия. На каждом мы зададим вопрос «а что будет, если…», посмотрим на ответ и объясним его. Правила лаборатории, без которых все опыты были бы бессмысленны, записаны в конце — вместе с опытом, который эти правила едва не ломает.

Закон вместо ответа

Начнём с кофе. В 1701 году в «Философских трудах» Королевского общества без подписи вышла заметка, которую потом признали работой Ньютона: тело остывает со скоростью, пропорциональной разности его температуры и температуры воздуха. Если в комнате $20°$, а температура кофе через $t$ минут равна $T(t)$, закон записывается так:

$$T'(t) = -k\,\bigl(T(t) - 20\bigr), \qquad k > 0.$$

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

Дифференциальное уравнение — уравнение, в котором неизвестна функция, а в запись входят её производные. Наибольший порядок производной называют порядком уравнения. Решение — функция, которая обращает уравнение в верное равенство при всех $t$ из некоторого промежутка. Уравнение вместе с начальным условием $y(t_0) = y_0$ называют задачей Коши; для уравнения второго порядка к нему добавляют ещё и начальную скорость $y'(t_0)$.

Здесь неизвестная функция зависит от одной переменной, и такие уравнения называют обыкновенными. Уравнения, где функция зависит от нескольких переменных, как температура от времени и точки в комнате, называют уравнениями в частных производных, и они сложнее. Самое простое обыкновенное уравнение мы уже встречали в главе 29: скорость роста пропорциональна самой величине.

Функция $y(t)$, определённая на промежутке, удовлетворяет уравнению $y' = ky$ тогда и только тогда, когда $y = Ce^{kt}$ при некоторой постоянной $C$. У задачи Коши $y(t_0) = y_0$ ровно одно решение: $y = y_0e^{k(t - t_0)}$.

Функции $Ce^{kt}$ подходят: $(Ce^{kt})' = kCe^{kt}$. Обратно, пусть $y' = ky$. Хитрость в том, чтобы поделить $y$ на ожидаемый ответ $e^{kt}$ и убедиться, что частное не меняется. По правилу производной произведения $\left(y e^{-kt}\right)' = y'e^{-kt} - ky\,e^{-kt} = (y' - ky)\,e^{-kt} = 0$. Функция с нулевой производной на промежутке постоянна (это следствие теоремы Лагранжа из главы 27), поэтому $y e^{-kt} = C$ и $y = Ce^{kt}$. Условие $y(t_0) = y_0$ даёт $C = y_0e^{-kt_0}$, и никакого другого $C$ оно не допускает.

Начальное значение — то, что было в момент $t_0$.

Коэффициент в законе $y' = ky$: при $k > 0$ рост с временем удвоения $\frac{\ln 2}{k}$, при $k < 0$ распад с периодом полураспада $\frac{\ln 2}{|k|}$.

Момент, с которого начали отсчёт.

Пример: для кофе обозначим $u = T - 20$. Тогда $u' = T' = -ku$, и $u = u_0e^{-kt}$. При $T(0) = 90°$ и $k = 0{,}05$ в минуту получаем $T = 20 + 70e^{-0{,}05t}$. Через $10$ минут $T = 20 + 70e^{-0{,}5} \approx 62{,}5°$, а до $60°$ кофе остынет за $\frac{\ln(70/40)}{0{,}05} \approx 11{,}2$ минуты.

Стенд первый: поле направлений

Большинство уравнений вида $y' = f(t, y)$ формулой не решаются. Но даже без формулы уравнение знает о своих решениях всё. В каждой точке $(t, y)$ оно сообщает наклон, с которым через эту точку проходит решение. Нарисуем в узлах сетки короткие чёрточки с этими наклонами — получится картина течения, а решения будут линиями тока.

Поле направлений уравнения $y' = f(t, y)$ — набор коротких отрезков: в каждой точке $(t, y)$ отрезок с наклоном $f(t, y)$. График решения называют интегральной кривой; в каждой своей точке он касается чёрточки поля.

Коснитесь плоскости — из этой точки вперёд и назад пойдёт решение. Переключайте уравнения или впишите своё. Почему у $y' = t - y$ все решения прижимаются к одной прямой, а у $y' = y^2$ улетают вверх, не дойдя до края?

У уравнения $y' = t - y$ все интегральные кривые прижимаются к прямой $y = t - 1$. Она сама решение: наклон у неё $1$, и $t - (t - 1) = 1$. Остальные решения — $y = t - 1 + Ce^{-t}$ (проверьте подстановкой), и добавка $Ce^{-t}$ тает. У логистического уравнения $y' = y(1 - y)$ видны две горизонтальные прямые, $y = 0$ и $y = 1$: на них наклон нулевой, и решение, начавшееся там, никуда не уходит. Такие решения-константы называют положениями равновесия. К прямой $y = 1$ соседние кривые стягиваются, от прямой $y = 0$ разбегаются: одно равновесие устойчиво, другое нет.

Стенд второй: шаг за шагом

Поле подсказывает, как найти решение без формулы: встать в начальную точку, пройти немного вдоль чёрточки, посмотреть на новую чёрточку и шагнуть снова. Этот рецепт Леонард Эйлер описал в «Интегральном исчислении» (1768).

Метод Эйлера — приближённое решение задачи Коши $y' = f(t, y)$, $y(t_0) = y_0$ шагами длины $h$: каждый шаг делается вдоль касательной, которую задаёт поле в текущей точке.

Где мы сейчас — приближение к $y(t_n)$.

Длина шага по времени. Чем она меньше, тем точнее и тем дольше считать.

Наклон поля в текущей точке: уравнение говорит, куда идти.

Пример: $y' = y$, $y(0) = 1$, шаг $h = \frac1n$. Каждый шаг умножает $y$ на $1 + \frac1n$, и после $n$ шагов $y(1) \approx \left(1 + \frac1n\right)^n$ — те самые сложные проценты из главы 29. При $n = 10$ выходит $2{,}594$ вместо $e = 2{,}718$.

Метод Эйлера ошибается: он идёт по касательной, а решение изгибается. Насколько ошибается, видно на том же примере. При $n = 10$ ошибка $0{,}125$, при $n = 100$ — $0{,}0135$, при $n = 1000$ — $0{,}00136$: шаг в десять раз меньше — ошибка в десять раз меньше. Причина — в ряде Тейлора из прошлой главы. Отклонение касательной от кривой за один шаг примерно равно $\frac{y''h^2}{2}$, а шагов до цели $\frac1h$, так что ошибки набегает порядка $h$. Для $y' = y$ это можно проверить точно: $\ln\left(1 + h\right)^{1/h} = \frac{1}{h}\left(h - \frac{h^2}{2} + \dots\right) = 1 - \frac h2 + \dots$, поэтому $\left(1 + h\right)^{1/h} \approx e\left(1 - \frac h2\right)$, и ошибка близка к $\frac{e}{2}h \approx 1{,}36h$.

Методы, у которых ошибка убывает быстрее, берут на каждом шаге несколько наклонов и усредняют их. Самый известный из них, классический метод Рунге — Кутты (Карл Рунге, 1895; Мартин Кутта, 1901), с десятью шагами даёт $e$ с ошибкой $2 \cdot 10^{-6}$, а при уменьшении шага вдвое ошибка падает примерно в $16$ раз. Им, в разных вариантах, решают уравнения почти во всех виджетах этой главы.

Меняйте шаг и сравнивайте ломаную Эйлера с точным решением. Включите метод Рунге — Кутты с тем же шагом. А на уравнении колебаний посмотрите, что Эйлер делает с амплитудой.

Последний пример стоит запомнить. Для пружины $y'' = -y$ метод Эйлера на каждом шаге чуть увеличивает энергию, и вместо вечных колебаний выходит раскручивающаяся спираль: амплитуда растёт в $\sqrt{1 + h^2}$ раз за шаг. Численный метод может нарушать законы сохранения, которые соблюдает точное решение, поэтому для долгих расчётов орбит планет придуманы специальные методы, которые энергию не накачивают.

Стенд третий: рост с потолком

Закон $y' = ky$ при $k > 0$ обещает бесконечный рост, а у любой популяции кончается еда и место. В 1838 году бельгийский математик Пьер Ферхюльст предложил тормозить рост множителем, который обращается в ноль, когда популяция достигает предела.

Скорость размножения, пока особей мало и места вдоволь: тогда множитель в скобках близок к $1$, и $y' \approx ry$.

Ёмкость среды — сколько особей она может прокормить. При $y = K$ скорость роста равна нулю, при $y > K$ популяция убывает.

Пример: на острове кролики с $r = 0{,}5$ в месяц, остров прокормит $K = 1000$ особей, а завезли $y_0 = 10$. Половины ёмкости популяция достигнет через $\frac{\ln 99}{0{,}5} \approx 9{,}2$ месяца (это следует из формулы ниже).

Уравнение $y' = ry\left(1 - \frac yK\right)$ с $r, K > 0$ называют логистическим, а графики его решений — логистическими кривыми.

Если $0 < y_0 < K$, то решение задачи Коши $y' = ry\left(1 - \frac yK\right)$, $y(0) = y_0$ единственно и равно

$$y(t) = \frac{K}{1 + \left(\frac{K}{y_0} - 1\right)e^{-rt}}.$$

Хитрость в том, чтобы перевернуть неизвестную: для $u = \frac1y$ уравнение становится линейным. Пусть $y$ — решение, и пока $y > 0$. Тогда $u' = -\frac{y'}{y^2} = -\frac{r}{y}\left(1 - \frac{y}{K}\right) = -ru + \frac rK$. Для $v = u - \frac1K$ отсюда $v' = u' = -r\left(u - \frac1K\right) = -rv$, и по теореме об экспоненциальном росте $v = \left(\frac1{y_0} - \frac1K\right)e^{-rt}$. Возвращаемся: $y = \frac1u = \frac{1}{\frac1K + \left(\frac1{y_0} - \frac1K\right)e^{-rt}}$; умножив числитель и знаменатель на $K$, получаем формулу из условия. Все шаги обратимы: для этой функции $u = \frac1y$ удовлетворяет уравнению $u' = -ru + \frac rK$, а значит, сама она — логистическому. Знаменатель этой формулы положителен и конечен при всех $t$, так что $y$ нигде не обращается в ноль, и оговорка «пока $y > 0$» выполнена всегда. Значит, любое решение с тем же начальным условием совпадает с найденным.

Кривая похожа на растянутую букву S. Сначала рост почти экспоненциальный, потом замедляется и прижимается к $K$. Быстрее всего популяция растёт при $y = \frac K2$: правая часть $ry\left(1 - \frac yK\right)$ — парабола по $y$ с вершиной в середине. Логистические кривые неплохо описывают и распространение новых технологий, и рост колоний бактерий, а дискретный родственник этого уравнения, логистическое отображение, в главе 59 откроет дверь в хаос. Поле этого уравнения есть на первом стенде.

Стенд четвёртый: пружина

Груз массы $m$ на пружине. Пружина тянет его к положению равновесия с силой $-\kappa x$, где $x$ — смещение, а $\kappa$ — жёсткость (закон Гука). Трение о воздух тормозит с силой $-cx'$. По второму закону Ньютона $mx'' = -\kappa x - cx'$. Поделим на $m$ и обозначим $\frac{\kappa}{m} = \omega^2$, $\frac{c}{m} = 2\gamma$:

$$x'' + 2\gamma x' + \omega^2 x = 0.$$

Это уравнение второго порядка, и для задачи Коши нужны два числа: где груз в начале и с какой скоростью движется. Без трения, при $\gamma = 0$, решения знакомы из механики: $x = A\cos\omega t + B\sin\omega t$, колебания с периодом $\frac{2\pi}{\omega}$. Груз $0{,}5$ кг на пружине жёсткостью $50$ Н/м даёт $\omega = 10$ рад/с и период около $0{,}63$ с. А с трением? Эйлер нашёл приём, который работает для всех линейных уравнений с постоянными коэффициентами: искать решение в виде $x = e^{\lambda t}$. Подставим: $\lambda^2 e^{\lambda t} + 2\gamma\lambda e^{\lambda t} + \omega^2 e^{\lambda t} = 0$, и после деления на $e^{\lambda t}$ остаётся квадратное уравнение.

Трение: коэффициент при $x'$.

Жёсткость на единицу массы: коэффициент при $x$.

Дискриминант решает, какой будет картина: колебания, критический случай или плавное сползание.

Пример: $x'' + 2x' + 5x = 0$ даёт $\lambda^2 + 2\lambda + 5 = 0$ и $\lambda = -1 \pm 2i$. Мнимая часть $2$ — частота колебаний, действительная часть $-1$ — скорость затухания: решения имеют вид $e^{-t}(A\cos 2t + B\sin 2t)$.

Для уравнения $x'' + px' + qx = 0$ с постоянными коэффициентами уравнение $\lambda^2 + p\lambda + q = 0$ называют характеристическим. Это характеристический многочлен из главы 38: если записать уравнение как систему для $x$ и $v = x'$, её матрица $\left(\begin{smallmatrix} 0 & 1 \\ -q & -p \end{smallmatrix}\right)$ имеет именно этот многочлен.

Пусть $\lambda_1$, $\lambda_2$ — корни характеристического уравнения для $x'' + px' + qx = 0$ с действительными $p$ и $q$. Все действительные решения этого уравнения таковы:

$x = Ae^{\lambda_1 t} + Be^{\lambda_2 t}$, если корни действительны и различны;
$x = (A + Bt)\,e^{\lambda t}$, если корни совпали: $\lambda_1 = \lambda_2 = \lambda$;
$x = e^{\alpha t}(A\cos\beta t + B\sin\beta t)$, если корни комплексные: $\lambda_{1,2} = \alpha \pm i\beta$, $\beta \ne 0$.

Здесь $A$ и $B$ — любые действительные числа, и при любых $x(0)$, $x'(0)$ они находятся однозначно.

Идея: разложить уравнение второго порядка на два уравнения первого, каждое из которых мы уже умеем решать.

По теореме Виета (глава 10) $p = -(\lambda_1 + \lambda_2)$ и $q = \lambda_1\lambda_2$. Поэтому $x'' + px' + qx = x'' - (\lambda_1 + \lambda_2)x' + \lambda_1\lambda_2 x$. Если обозначить $w = x' - \lambda_2 x$, то $w' - \lambda_1 w = (x'' - \lambda_2 x') - \lambda_1(x' - \lambda_2 x)$ — ровно то же выражение. Значит, для решения $x$ функция $w$ удовлетворяет уравнению $w' = \lambda_1 w$, и по теореме об экспоненциальном росте $w = Ce^{\lambda_1 t}$. Если корни комплексные, функции тоже приходится брать комплекснозначными. Доказательство той теоремы работает и для них: по формуле Эйлера из главы 29 $e^{(a + ib)t} = e^{at}(\cos bt + i\sin bt)$, и прямое дифференцирование даёт $\left(e^{\lambda t}\right)' = \lambda e^{\lambda t}$ для комплексных $\lambda$.

Теперь решаем $x' - \lambda_2 x = Ce^{\lambda_1 t}$. Умножим на $e^{-\lambda_2 t}$: слева получится производная произведения, $\left(xe^{-\lambda_2 t}\right)' = Ce^{(\lambda_1 - \lambda_2)t}$. Если $\lambda_1 \ne \lambda_2$, первообразная правой части — $\frac{C}{\lambda_1 - \lambda_2}e^{(\lambda_1 - \lambda_2)t} + D$, и $x = Ae^{\lambda_1 t} + De^{\lambda_2 t}$. Если $\lambda_1 = \lambda_2 = \lambda$, правая часть равна просто $C$, и $xe^{-\lambda t} = Ct + D$, то есть $x = (D + Ct)e^{\lambda t}$. Обратно, все такие функции — решения: достаточно пройти эти выкладки в обратную сторону.

При комплексных корнях $\lambda_2 = \bar\lambda_1$, и действительное решение $x$ совпадает со своей действительной частью. Действительная часть функции $Ae^{(\alpha + i\beta)t}$ — это $e^{\alpha t}(a\cos\beta t + b\sin\beta t)$ с какими-то действительными $a$, $b$, и такой же вид у второго слагаемого. Наконец, постоянные: при различных корнях $x(0) = A + B$, $x'(0) = \lambda_1 A + \lambda_2 B$ — система с определителем $\lambda_2 - \lambda_1 \ne 0$, она однозначно решается; при кратном корне $x(0) = A$, $x'(0) = \lambda A + B$.

Для пружины это три режима. При $\gamma < \omega$ корни комплексные, и груз колеблется с частотой $\sqrt{\omega^2 - \gamma^2}$ внутри затухающей огибающей $e^{-\gamma t}$. При $\gamma > \omega$ корни действительные и отрицательные: груз, отпущенный из отклонённого положения, сползает к равновесию, ни разу его не перейдя. Граница, $\gamma = \omega$, — критическое затухание с двойным корнем $-\omega$ и решением $(A + Bt)e^{-\omega t}$. Это тот самый множитель $t$, о котором говорила глава 38: при двойном корне матрица системы не диагонализуема. При критическом затухании груз возвращается быстрее всего, не качаясь, поэтому стрелочные приборы и дверные доводчики настраивают близко к этому режиму.

Оттяните груз и отпустите. Справа — корни характеристического уравнения на комплексной плоскости: увеличивайте трение и следите, как они идут по дуге навстречу друг другу, сталкиваются в точке $-\omega$ и разбегаются по действительной оси.

Найдите решение задачи Коши $x'' + 2x' + 5x = 0$, $x(0) = 0$, $x'(0) = 2$.

Корни характеристического уравнения $\lambda^2 + 2\lambda + 5 = 0$ — это $-1 \pm 2i$, поэтому $x = e^{-t}(A\cos 2t + B\sin 2t)$. Из $x(0) = 0$ получаем $A = 0$. Тогда $x' = e^{-t}(-B\sin 2t + 2B\cos 2t)$, и $x'(0) = 2B = 2$, $B = 1$. Ответ: $x = e^{-t}\sin 2t$ — колебания с периодом $\pi$ внутри огибающей $\pm e^{-t}$.

Набить руку на этих уравнениях — от роста и распада до задачи Коши для пружины — можно в тренажёре.

Стенд пятый: волки и зайцы

После Первой мировой войны итальянский биолог Умберто Д'Анкона разбирал статистику рыбных рынков Адриатики и заметил странность. В военные годы, когда ловили мало, доля хищных рыб в уловах выросла. Почему ослабление промысла помогло хищникам больше, чем их добыче? Д'Анкона спросил своего тестя, математика Вито Вольтерру, и тот в 1926 году построил модель. Годом раньше ту же модель независимо опубликовал американец Альфред Лотка.

Жертвы $x$ (зайцы) без хищников размножались бы экспоненциально со скоростью $\alpha$.

Встречи с хищниками: их число пропорционально произведению $xy$, и каждая уменьшает число жертв.

Те же встречи кормят хищников $y$ (волков) и дают им потомство.

Без еды хищники вымирают со скоростью $\gamma$.

Пример: $\alpha = \gamma = 0{,}6$ в год, $\beta = 0{,}03$, $\delta = 0{,}015$. Если зайцев $x = 40$, а волков $y = 20$, обе скорости равны нулю: $0{,}6 \cdot 40 - 0{,}03 \cdot 40 \cdot 20 = 0$ и $0{,}015 \cdot 40 \cdot 20 - 0{,}6 \cdot 20 = 0$. Это положение равновесия.

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

Фазовая плоскость системы двух уравнений — плоскость, точки которой — состояния $(x, y)$. Решение изображается движущейся точкой, её путь называют фазовой траекторией, а картину всех траекторий — фазовым портретом. Состояние, в котором все скорости равны нулю, называют положением равновесия.

Коснитесь фазовой плоскости — точка побежит по своей траектории, а внизу будут видны графики обеих величин. Двигайте «промысел» и следите за крестиком равновесия. Переключитесь на маятник: чем замкнутые петли отличаются от волнистых линий сверху и снизу?

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

Для любого решения системы Лотки — Вольтерры с $x, y > 0$ величина $V(x, y) = \delta x - \gamma\ln x + \beta y - \alpha\ln y$ не меняется со временем.

Продифференцируем $V$ вдоль решения по правилу для сложной функции и подставим уравнения: $V' = \left(\delta - \frac{\gamma}{x}\right)x' + \left(\beta - \frac{\alpha}{y}\right)y' = \left(\delta - \frac\gamma x\right)x(\alpha - \beta y) + \left(\beta - \frac\alpha y\right)y(\delta x - \gamma)$. Раскроем скобки: $(\delta x - \gamma)(\alpha - \beta y) + (\beta y - \alpha)(\delta x - \gamma) = (\delta x - \gamma)\bigl((\alpha - \beta y) + (\beta y - \alpha)\bigr) = 0$. Производная равна нулю на промежутке, значит, $V$ постоянна.

Функция $\delta x - \gamma\ln x$ убывает до $x = \frac\gamma\delta$ и потом растёт, так же ведёт себя $\beta y - \alpha\ln y$ около $y = \frac\alpha\beta$. Поэтому $V$ устроена как чаша с дном в положении равновесия, а её линии уровня — замкнутые кривые вокруг дна. Траектория обязана оставаться на одной такой кривой, и точка ходит по ней кругами. Около равновесия период близок к $\frac{2\pi}{\sqrt{\alpha\gamma}}$; в нашем примере это $\frac{2\pi}{0{,}6} \approx 10{,}5$ года. Похожие десятилетние циклы видны в знаменитых данных о шкурках рыси и зайца-беляка, которые скупала канадская Компания Гудзонова залива, хотя настоящая экология, конечно, сложнее двух уравнений.

Средние значения $x$ и $y$ за полный цикл равны координатам равновесия: $\bar x = \frac\gamma\delta$, $\bar y = \frac\alpha\beta$.

Разделим первое уравнение на $x$: $(\ln x)' = \frac{x'}{x} = \alpha - \beta y$. Проинтегрируем по циклу длины $T$: слева $\ln x(T) - \ln x(0) = 0$, потому что за цикл точка вернулась, а справа $\alpha T - \beta\int_0^T y\,dt$. Значит, среднее $\bar y = \frac1T\int_0^T y\,dt = \frac\alpha\beta$. Точно так же из $(\ln y)' = \delta x - \gamma$ получаем $\bar x = \frac\gamma\delta$.

Теперь загадка Д'Анконы решается в две строки. Промысел отнимает у обеих популяций одинаковую долю $f$ в единицу времени: $\alpha$ превращается в $\alpha - f$, а $\gamma$ — в $\gamma + f$. Среднее число жертв $\frac{\gamma + f}{\delta}$ от промысла растёт, а среднее число хищников $\frac{\alpha - f}{\beta}$ падает. Война сократила промысел, $f$ уменьшилось, и хищников стало больше. Этот вывод называют принципом Вольтерры.

На том же стенде есть маятник. Его уравнение $\theta'' = -\frac gL\sin\theta$ нелинейно, но при малых углах $\sin\theta \approx \theta$ (ряд из прошлой главы), и получается пружина с $\omega = \sqrt{g/L}$ и периодом $2\pi\sqrt{L/g}$. Для метрового маятника это $2{,}006$ с. При больших размахах приближение портится: при отклонении на $90°$ период на $18\,\%$ больше. На фазовом портрете маятника видны и колебания (замкнутые петли), и вращение через верх (волнистые линии), а разделяет их особая траектория, ведущая в перевёрнутое положение. Замкнутость петель объясняет тот же приём, что у Вольтерры: величина $\frac{(\theta')^2}{2} - \frac gL\cos\theta$, то есть энергия на единицу массы, делённая на $L^2$, вдоль решения не меняется.

Стенд шестой: эпидемия

В 1927 году Уильям Кермак и Андерсон Маккендрик предложили модель, которой пользуются до сих пор. Население делится на три группы: восприимчивые $S$, которые ещё не болели, заразные $I$ и выбывшие $R$, которые выздоровели и больше не заражаются. Все три величины — доли населения, $S + I + R = 1$.

Восприимчивые заражаются при встречах с заразными; число таких встреч пропорционально $SI$.

Заразных прибавляется ровно столько, сколько заразилось, и убавляется за счёт выздоровевших.

Выздоровевшие копятся.

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

Скорость выздоровления: в среднем человек заразен $\frac1\gamma$ дней.

Пример: $\beta = 0{,}25$ в день, заразность длится в среднем $10$ дней, $\gamma = 0{,}1$. Один больной за время болезни заражает $\frac{0{,}25}{0{,}1} = 2{,}5$ человека.

Базовое репродуктивное число $R_0 = \frac\beta\gamma$ — сколько человек в среднем заражает один больной за время болезни, если все вокруг восприимчивы.

Число заразных растёт в тот момент, когда $R_0 S > 1$, и убывает, когда $R_0 S < 1$. В частности, вспышка начинается, только если $R_0 S(0) > 1$.

Вынесем $I$ во втором уравнении: $I' = I(\beta S - \gamma) = \gamma I\,(R_0 S - 1)$. Множитель $\gamma I$ положителен, пока заразные есть, и знак $I'$ совпадает со знаком $R_0 S - 1$. А так как $S' = -\beta SI \le 0$, доля восприимчивых только убывает: если в начале $R_0 S < 1$, так будет и дальше, и число заразных убывает с первого дня.

Отсюда главная формула вакцинации. Если заранее привить долю $p$ населения, то $S(0) \approx 1 - p$, и вспышки не будет при $R_0(1 - p) < 1$, то есть при $p > 1 - \frac1{R_0}$. Для $R_0 = 2{,}5$ это $60\,\%$. Для кори $R_0$ обычно оценивают в $12$–$18$, и порог коллективного иммунитета выходит около $92$–$94\,\%$. Но если эпидемия уже идёт, она не останавливается на пороге: сколько человек переболеет в итоге, отвечает следующая теорема.

Пусть в начале $R(0) = 0$, $S(0) = s_0$ и $I(0) = 1 - s_0 > 0$. Тогда $S(t)$ стремится к числу $s_\infty > 0$, для которого $s_\infty = s_0\,e^{-R_0(1 - s_\infty)}$, а $I(t) \to 0$.

Идея: исключить время и следить за эпидемией как за точкой на плоскости $(S, I)$.

Поделим второе уравнение на первое: $\frac{dI}{dS} = \frac{\beta SI - \gamma I}{-\beta SI} = \p1{-1 + \frac{1}{R_0 S}}$. Пока $I > 0$, доля $S$ строго убывает, и траекторию можно считать графиком функции $I(S)$ с таким наклоном.

Первообразная: $I = -S + \frac{1}{R_0}\ln S + C$. В начале $S = s_0$, $I = 1 - s_0$, откуда $C = 1 - \frac{1}{R_0}\ln s_0$. Вся эпидемия идёт по одной кривой $\p2{I = 1 - S + \frac{1}{R_0}\ln\frac{S}{s_0}}$, справа налево.

Наклон равен нулю при $S = \frac{1}{R_0}$: пик эпидемии наступает ровно тогда, когда доля восприимчивых падает до порога коллективного иммунитета. При $s_0 \approx 1$ и $R_0 = 2{,}5$ на пике болеют $\p3{1 - \frac{1}{R_0} - \frac{\ln R_0}{R_0}} \approx 23\,\%$.

Кривая пересекает ось $I = 0$ левее пика, в точке $s_\infty > 0$. Точка $(S, I)$ туда не проскочит: $I$ остаётся положительным, ведь уравнение $I' = I(\beta S - \gamma)$ линейно по $I$, и $I = I(0)\,e^{\int(\beta S - \gamma)dt} > 0$. Убывающая и ограниченная снизу доля $S$ сходится; предел не может быть правее $s_\infty$, иначе $I$ оставалось бы больше некоторого $\varepsilon > 0$, а $S' = -\beta SI$ — меньше отрицательного числа, и $S$ ушла бы в минус. Значит, $S \to s_\infty$ и $I \to 0$.

В точке $s_\infty$ кривая даёт $0 = 1 - s_\infty + \frac1{R_0}\ln\frac{s_\infty}{s_0}$, то есть $s_\infty = s_0 e^{-R_0(1 - s_\infty)}$. При $R_0 = 2{,}5$ и $s_0 \approx 1$ корень $s_\infty \approx 0{,}107$: переболеют $89\,\%$, намного больше порога в $60\,\%$. Эпидемия по инерции «перелетает» через порог.

Меняйте $R_0$, длительность болезни и долю привитых. Включите ограничение контактов с 30-го дня: пик станет ниже и позже. Что будет, если ограничения снять слишком рано?

Для болезни с $R_0 = 4$ какую долю населения нужно привить заранее, чтобы вспышка не началась? Ответ дайте дробью.

Нужно $R_0(1 - p) < 1$, то есть $1 - p < \frac14$ и $p > \frac34$. Порог — $\frac34$, или $75\,\%$ населения.

Правила лаборатории

До сих пор у каждой задачи Коши было ровно одно решение. Это молчаливое правило лаборатории: запустите опыт дважды с одинаковых начальных условий — и получите одно и то же. Всегда ли так? Нет.

Возьмём уравнение $y' = 3\sqrt[3]{y^2}$ с условием $y(0) = 0$. Ему удовлетворяет $y = 0$ — и $y = t^3$ тоже: $(t^3)' = 3t^2 = 3\sqrt[3]{t^6}$. Более того, годится функция, которая стоит в нуле сколько угодно долго, а потом при $t = c$ трогается по кривой $(t - c)^3$. Решение не единственно, и лаборатория с таким законом непредсказуема в самом прямом смысле. Другая беда: у $y' = y^2$ с $y(0) = 1$ решение $y = \frac{1}{1 - t}$ уходит в бесконечность при $t = 1$, хотя правая часть — безобидная парабола. Решение существует лишь на конечном промежутке. Всё дело в том, насколько резко $f(t, y)$ может меняться при изменении $y$.

Функция $f(t, y)$ удовлетворяет условию Липшица по $y$ с константой $K$, если $|f(t, y_1) - f(t, y_2)| \le K|y_1 - y_2|$ для всех рассматриваемых $t$, $y_1$, $y_2$. Так будет, например, если $f$ имеет по $y$ производную, ограниченную по модулю числом $K$: это теорема Лагранжа. У $3\sqrt[3]{y^2}$ около нуля наклон бесконечен, и условия Липшица нет.

Пусть функция $f(t, y)$ непрерывна в прямоугольнике $|t - t_0| \le a$, $|y - y_0| \le b$ и удовлетворяет в нём условию Липшица по $y$. Тогда на некотором промежутке $|t - t_0| \le h$ задача Коши $y' = f(t, y)$, $y(t_0) = y_0$ имеет решение, и притом только одно.

Идея: заменить уравнение равносильным интегральным и получать решение последовательными приближениями, каждое из которых подставляется в правую часть. Непрерывная функция $y$ решает задачу Коши тогда и только тогда, когда $y(t) = y_0 + \int_{t_0}^{t} f(s, y(s))\,ds$: в одну сторону — формула Ньютона — Лейбница, в другую — дифференцирование интеграла по верхнему пределу (глава 28). Пусть $|f| \le M$ в прямоугольнике; возьмём $h = \min\left(a, \frac bM\right)$ и рассмотрим $t \ge t_0$ (влево всё так же). Любое решение при $t_0 \le t \le t_0 + h$ не выходит из прямоугольника: пока оно внутри, $|y'| \le M$, и за время $h$ оно не может отойти от $y_0$ дальше чем на $Mh \le b$.

Единственность. Пусть $y$ и $z$ — два решения, $u(t) = |y(t) - z(t)|$ и $U$ — наибольшее значение $u$ на $[t_0;\,t_0 + h]$. Вычтем интегральные равенства и применим условие Липшица: $u(t) \le K\int_{t_0}^t u(s)\,ds$. Отсюда $u(t) \le KU(t - t_0)$; подставим это в интеграл снова: $u(t) \le K^2U\frac{(t - t_0)^2}{2}$, и по индукции $u(t) \le U\frac{(K(t - t_0))^n}{n!}$ при любом $n$. Правая часть — общий член сходящегося ряда для $Ue^{K(t - t_0)}$, она стремится к нулю (глава 30). Значит, $u = 0$: решения совпадают.

Существование. Положим $y_0(t) = y_0$ и $y_{n+1}(t) = y_0 + \int_{t_0}^t f(s, y_n(s))\,ds$. Из $|f| \le M$ и выбора $h$ следует, что все приближения остаются в прямоугольнике: $|y_{n+1}(t) - y_0| \le M(t - t_0) \le b$. Сначала $|y_1(t) - y_0(t)| \le M(t - t_0)$, а дальше каждый шаг добавляет множитель $K$ по условию Липшица и интеграл: $|y_{n+1}(t) - y_n(t)| \le K\int_{t_0}^t |y_n(s) - y_{n-1}(s)|\,ds$. По индукции $|y_{n+1}(t) - y_n(t)| \le \frac{M K^n (t - t_0)^{n+1}}{(n + 1)!}$. Поэтому ряд $y_0 + (y_1 - y_0) + (y_2 - y_1) + \dots$ сходится при каждом $t$ по признаку сравнения, и приближения стремятся к некоторой функции $y(t)$, причём отличаются от неё меньше чем на хвост числового ряда $\sum\frac{MK^nh^{n+1}}{(n+1)!}$, одинаковый для всех $t$. Остаётся перейти к пределу в равенстве $y_{n+1} = y_0 + \int f(s, y_n)$: интеграл отличается от $\int f(s, y)$ не больше чем на $Kh$, умноженное на этот хвост. Для этого нужно ещё знать, что предельная функция непрерывна. Так и есть, потому что сходимость одинаково быстрая при всех $t$; её называют равномерной, и соответствующую теорему докажет глава 53. Вместе с ней доказательство полное.

Итерации Пикара можно увидеть в деле. Для $y' = y$, $y(0) = 1$ получается $y_1 = 1 + t$, $y_2 = 1 + t + \frac{t^2}{2}$, $y_3 = 1 + t + \frac{t^2}{2} + \frac{t^3}{6}$ — многочлены Тейлора экспоненты из прошлой главы, которые сходятся к $e^t$.

Могут ли две разные интегральные кривые уравнения $y' = t - y$ пересечься?

Функция $t - y$ удовлетворяет условию Липшица с $K = 1$. Если бы две кривые имели общую точку $(t_0, y_0)$, обе были бы решениями одной задачи Коши, а по теореме Пикара решение одно. Так что интегральные кривые не пересекаются и не касаются — ровно это видно на первом стенде.

Теорема Пикара — математическая форма детерминизма. Пьер-Симон Лаплас в 1814 году писал о воображаемом разуме, который знал бы положения и скорости всех частиц мира и потому видел бы будущее, как прошлое. По теореме Пикара такой разум действительно мог бы вычислить будущее для любой системы с гладкими законами. Но только если бы знал начальные условия абсолютно точно. Последний стенд показывает, чего стоит эта оговорка.

Стенд седьмой: бабочка

В 1961 году метеоролог Эдвард Лоренц из Массачусетского технологического института гонял на компьютере упрощённую модель атмосферы. Чтобы повторить расчёт с середины, он ввёл промежуточное значение из распечатки, где оно было округлено: $0{,}506$ вместо $0{,}506127$. Разница — около одной десятитысячной, меньше обычной ошибки измерений. Сначала повторный расчёт шёл след в след за первым, а потом модельная погода пошла совсем по-другому. В 1963 году Лоренц опубликовал ещё более простую систему из трёх уравнений, где видно то же самое:

$$x' = \sigma(y - x), \qquad y' = x(\rho - z) - y, \qquad z' = xy - \beta z,$$

с числами $\sigma = 10$, $\rho = 28$, $\beta = \frac83$. Правые части — многочлены, условие Липшица выполнено, и теорема Пикара гарантирует: каждое начальное состояние определяет будущее однозначно. Посмотрите, что это значит на деле.

Два запуска системы Лоренца, начальные значения $x$ отличаются на $\delta$. Уменьшайте $\delta$ в десять, сто, миллион раз и следите, насколько позже расходятся кривые.

Каждое уменьшение начальной ошибки в десять раз отодвигает расхождение лишь примерно на $2{,}5$ единицы времени. Ошибка растёт экспоненциально, как решение уравнения $y' = ky$, только $k$ здесь задаёт не внешний закон, а сама система. Каждая новая верная цифра начальных данных прибавляет к сроку прогноза одно и то же время: точность измерений приходится наращивать экспоненциально, а прогноз удлиняется лишь линейно.

Первая часть ответа. Уравнения, по которым движется атмосфера, детерминированы: по теореме Пикара настоящее однозначно задаёт будущее. Но однозначно — не значит предсказуемо. В системах вроде системы Лоренца малая неточность начальных данных растёт экспоненциально, и каждая новая верная цифра в измерениях продлевает прогноз лишь на фиксированный срок. Почему так происходит, как выглядит странный аттрактор Лоренца и откуда берётся предел предсказуемости погоды около двух недель — в главе 59.

Куда дальше

На стенде волков и зайцев решения шли по линиям уровня функции $V(x, y)$ — как туристы, которые обходят гору по горизонтали, не поднимаясь и не спускаясь. Функции нескольких переменных встречались нам в этой главе на каждом шагу: $f(t, y)$ в поле направлений, $V(x, y)$ у Вольтерры, SIR с тремя долями. Но смотрели мы на них лишь краем глаза. Представьте горный пейзаж: высота — функция двух координат. Как измерить крутизну склона, если она зависит от направления? Куда идти, чтобы подняться быстрее всего? Где у такой функции вершины, седловины и впадины? Ответ даёт глава о функциях многих переменных.