Практическое введение в решение дифференциальных уравнений в Python

Н. М. Ершов

Практическое 
введение в решение 
дифференциальных 
уравнений в Рython

Москва, 2022

УДК 517.9
ББК 22.161.61
Е80

Ершов Н. М.
Е80  Практическое введение в решение дифференциальных уравнений 
в Python. – М.: ДМК Пресс, 2022. – 176 с.: ил. 

ISBN 978-5-93700-147-4

Книга посвящена вопросам практического применения символьных вычислений 
для решения различных прикладных задач, приводящих к дифференциальным 
уравнениям и их системам, с использованием библиотеки символьных 
вычислений SymPy языка Python.
Издание ориентировано на школьников старших классов, студентов, преподавателей 
и всех, интересующихся тематикой математического моделирования. 
Может служить дополнением к классическим учебникам по теории обыкновенных 
дифференциальных уравнений.

УДК 517.9
ББК 22.161.61

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

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

 
© Н. М. Ершов, 2022
ISBN 978-5-93700-147-4 
©  Оформление, издание, ДМК Пресс, 2022

Оглавление

Введение. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .5

Глава

1

Вращение жидкости . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
Задача (7) Модель (9) Упражнения и замечания (18)

Глава

2

Водяные часы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Задача (21) Модель (22) Упражнения и замечания (29)

Глава

3

Элементарные химические реакции . . . . . . . . . . . . . . . . . . . . . . . . . . . .33
Закон действующих масс (33) Задача (35) Модель (35)
Упражнения и замечания (38)

Глава

4

Задача о четырех жуках . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
Кривые погони (41) Задача (41) Модель (43) Упражнения
и замечания (47)

Глава

5

Барометрическая формула . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
Задача (51) Модель (53) Упражнения и замечания (58)

Глава

6

Модели роста. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .61
Модель естественного роста (61) Модель Ферхюльста (62)
Модель (62) Упражнения и замечания (68)

Глава

7

Табулирование функций . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
Метод Эйлера (73) Задача (74) Модель (75) Упражнения
и замечания (81)

Глава

8

Ортогональные траектории . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
Ортогональные семейства кривых (85) Задача (85)
Модель (87) Упражнения и замечания (90)

Глава

9

Математическая вышивка . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
Задача (95) Уравнение Клеро (96) Модель (97)
Упражнения и замечания (101)

Глава

10

Криволинейные зеркала . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .105
Задача (105) Модель (106) Упражнения и замечания (112)

Оглавление

Глава

11

Из пушки на луну . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
Задача (115) Модель (116) Упражнения и замечания (120)

Глава

12

Метроном . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .123
Задача (123) Модель (124) Упражнения и замечания (130)

Глава

13

Пружинный маятник . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
Задача (133) Модель (134) Упражнения и замечания (140)

Глава

14

Модель Лотки–Вольтерры. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .143
Задача (143) Модель (144) Упражнения и замечания (149)

Глава

15

Системы реакций первого порядка . . . . . . . . . . . . . . . . . . . . . . . . . . . 153
Задача (153) Преобразование Лапласа (154) Модель (156)
Упражнения и замечания (160)

Глава

16

Геодезические линии . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
Задача (163) Уравнение Эйлера—Лагранжа (165)
Модель (166) Упражнения и замечания (171)

Введение

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

РИС.
1 SymPy — библиотека

Python для символьных вычислений


Программные системы, позволяющие пользователю 
работать с математическими формулами, выполняя 
над ними те или иные символьные преобразования, 
называются системами символьных вычислений,
или системами компьютерной алгебры. Первые такие
системы появились еще в 60-х годах прошлого века.
В настоящее время пользователям доступны десятки 
подобного рода систем — от коммерческих (Maple,
Mathematica) до систем с открытым исходным кодом
(Reduce, Maxima), некоторые из которых способны работать 
уже и на мобильных устройствах.
Система SymPy, которой и посвящена данная книга,
является по сути обычной библиотекой языка Python.
Такой подход к построению систем символьных вычислений 
имеет ряд преимуществ. Во-первых, система
оказывается открытой и доступной всем без исключения 
пользователям. Во-вторых, работа с библиотекой 
SymPy в среде Jupyter Notebook позволяет проводить 
символьные вычисления прямо в браузере либо
локально, либо удаленно с помощью какого-нибудь об-

Введение

лачного сервиса1, например Google Colab. В-третьих,
1 При таком сценарии пользователю 
вообще не нужно устанавливать 
локально никакое программное 
обеспечение.
символьные вычисления в рамках данной библиотеки 
оказываются интегрированными в обычную программу 
на языке Python, при этом пользователю системы 
оказывается доступным весь функционал языка 
Python, а также вся его инфраструктура в виде бесчисленного 
набора разнообразнейших пакетов и библиотек.

2
Заметим, что автором не ставилась 
цель дать полное описание 
библиотеки SymPy. Такая задача, 
с одной стороны, является 
просто неподъемной, библиотека 
SymPy состоит из большого числа 
модулей, только часть из которых 
посвящена решению дифференциальных 
уравнений. Кроме 
того, данная библиотека все
еще активно развивается, ее функционал 
постоянно расширяется и
модифицируется. Более полную и
актуальную информацию по работе

с
библиотекой
SymPy
читателю 
рекомендуется находить либо

на
официальном
сайте
библиотеки 
https://www.sympy.org, либо 
на профильных ресурсах в сети 
Интернет, например на сайте
https://stackoverflow.com.

Предлагаемая читателю книга устроена достаточно 
просто. Каждая ее глава посвящена рассмотрению
одной прикладной модели из физики, химии, биологии 
и т. д. После теоретического рассмотрения модели 
и вывода соответствующего ей дифференциального 
уравнения максимально детально описывается процесс 
формализации модели и решения возникающих
в ней дифференциальных уравнений с помощью библиотеки 
SymPy2. Каждая глава сопровождается набором 
упражнений как технического характера — посчитать 
интеграл, решить уравнение, построить график,
так и исследовательского — построить и исследовать
по описанной схеме аналогичную модель.
Автор с благодарностью примет все конструктивные 
отзывы, комментарии и замечания относительно 
структуры и содержания настоящего издания, которые 
читатель может оставить или в телеграм-чате
https://t.me/odesinsympy, или прислать автору на электронную 
почту по адресу ershovnm@gmail.com.

ГЛАВА 1

Вращение жидкости

Задача • В сосуд, имеющий форму прямого кругового 
цилиндра, налита жидкость, например вода.
Сосуд вращается с постоянной угловой скоростью ω
относительно оси цилиндра (рис. 1). Требуется опре-

y

x

ω

РИС.
1 Вращение
сосуда
с
жидкостью

1 «Достаточно долго» понимается 
здесь в том смысле, что вся
жидкость в сосуде должна прийти
в стационарное состояние относительно 
самого сосуда, т. е. каждый
элементарный ее объем будет совершать 
только общее вращательное 
движение с заданной угловой
скоростью ω.

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

α

α

x

y

y0

P

mg

N

Nx

Ny

O

РИС.
2 Силы,
действующие
на малый объем жидкости ν
массы m в точке P(x, y)

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

Введем в модель систему координат, как это показано 
на рис. 1. Ось Oy направим по оси цилиндра,
ось Ox — перпендикулярно оси Oy вдоль основания
цилиндра, начало координат, таким образом, оказывается 
в центре основания цилиндра.
Рассмотрим малый объем ν жидкости на ее поверхности (
точка P на рис. 2). Пусть его масса равна m. На
этот объем действуют две силы: сила тяжести F = mg
и сила реакции опоры N со стороны всей остальной
жидкости. Сила тяжести, как обычно, действует вертикально 
вниз, а реакция опоры — по нормали к нашей
искомой кривой, т. е. перпендикулярно касательной к
этой кривой в точке P.
Обозначим угол наклона данной касательной к оси
Ox через α и разложим силу N на горизонтальную Nx

Глава 1

и вертикальную Ny составляющие (рис. 2):

Nx = N sin α, Ny = N cos α.

Запишем теперь 2-й закон Ньютона отдельно для каждой 
координаты. Вдоль вертикальной оси наш объем
ν находится в состоянии равновесия, поэтому суммарная 
сила, действующая на него в вертикальном направлении, 
должна быть равна нулю:

N cos α − mg = 0.
(1)

Единственная сила, действующая на объем ν вдоль горизонтального 
направления, — это проекция реакции
опоры Nx. Так как данный объем совершает круговое
движение с радиусом вращения x, то он испытывает
центростремительное ускорение a, которое направлено 
горизонтально к оси вращения (рис. 3). Величина

a

ω

P

O

x

РИС. 3 Центростремительное
ускорение для точки P (вид
на сосуд сверху)

этого ускорения вычисляется по формуле (квадрат угловой 
скорости на радиус вращения):

a = ω2x.
(2)

Таким образом, вдоль оси Ox второй закон Ньютона
записывается в следующем виде:

−N sin α = −ma = −mω2x,
(3)

обе части взяты с отрицательным знаком, потому что
и сила Nx, и ускорение a направлены против положительного 
направления оси Ox.
Запишем формулы (1) и (3) в виде системы двух
уравнений:

Ox : N sin α = mω2x,
Oy : N cos α = mg.
(4)

Поделим первое уравнение системы на второе:

tg α = ω2x

g .

Теперь учтем тот факт, что согласно геометрическому 
смыслу производной функции y(x) величина производной 
в точке x0 равна тангенсу угла наклона касательной 
к соответствующей кривой y(x) в заданной
точке x0. То есть мы можем заменить tg α на y′, что и
даст нам искомое дифференциальное уравнение:

y′ = ω2x

g .
(5)

Вращение жидкости
9

Если дополнительно обозначить высоту жидкости
в точке x = 0 (на оси вращения) через y0, то можно
поставить для полученного дифференциального уравнения 
задачу Коши:

y′ = ω2x

g ,

y(0) = y0.
(6)

Модель • Рассмотрим теперь процесс решения поставленной 
начальной задачи с помощью библиотеки
SymPy. Сначала мы построим общее решение дифференциального 
уравнения (5), затем, используя начальное 
условие y(0) = y0, решим начальную задачу (6).
Следующим шагом определим параметр y0 из условия
постоянства объема вращаемой жидкости, рассмотрев
при этом два случая — подкритический, когда поверхность 
вращаемой жидкости не касается дна цилиндра,
и надкритический, когда такое касание происходит.

1 Построение модели начинаем с подключения библиотеки 
SymPy и функции display, определенной в модуле 
IPython.display.

1 from sympy import *

2 from IPython.display import display

2 С помощью команды symbols создаем символы x,
g, ω, y0 и C1, которые нам потребуются для задания и
решения дифференциального уравнения. Аргументом
2
Информация о положительности 
того или иного символа может 
быть использована библиотекой 
SymPy для упрощения формул,
в которые входит этот символ, например, 
при извлечении квадратных 
корней.

3 Для краткости вывод команды
display в данном случае мы поместили 
в одну строку, на самом
деле эта команда выводит каждое
из переденных ей выражений в отдельной 
строке.

этой команды является строка, содержащая вводимые
символы и определяющая, как эти символы будут выглядеть 
при печати. Возвращает эта команда кортеж
с объектами SymPy, которые нужно запомнить в соответствующих 
переменных для последующего использования. 
Так как параметр g, представляющий собой
ускорение свободного падения, является всегда строго 
положительным, то при создании соответствующего 
символа укажем этот факт, используя опциональный 
аргумент positive2. Используем команду display,
чтобы посмотреть, как выглядят определенные нами
символы. Например, видно, что SymPy корректно отображает 
греческие буквы, а цифры после символов преобразует 
в нижние индексы этих символов3.

1 x, omega, y0, C1 = symbols("x omega y0 C1")

2 g = symbols("g", positive = True)

3 display(x, g, omega, y0, C1)

Глава 1

x g ω y0 C1

3 Наше дифференциальное уравнение

y′ = ω2x

g

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

Создадим переменную ode_rhs, в которую и запишем
правую часть решаемого уравнения.

1 ode_rhs = omega ** 2 * x / g

2 display(ode_rhs)

ω2x

g

4 Общее решение простейшего дифференциального
уравнения y′ = f(x) — это неопределенный интеграл
от правой части f(x) этого уравнения. В библиотеке
SymPy неопределенные интегралы вычисляются с помощью 
команды integrate(f, x), где f — подынтегральное 
выражение, x — переменная интегрирования. Константа 
интегрирования C1 при этом автоматически командой 
integrate к найденному интегралу не прибавляется, 
это надо делать вручную5, поэтому мы прос-
5 Если это необходимо.
то к результату работы функции integrate прибавим
символ C1. Найденное общее решение дифференциального 
уравнения сохраним в переменной dsol.

1 dsol = integrate(ode_rhs, x) + C1

2 display(dsol)

C1 + ω2x2

2g

5 Стандартный способ определения константы интегрирования 
C1 в общем решении дифференциального 
уравнения из начального условия заключается
в подстановке в это решение величин из начального
условия, что приводит к алгебраическому уравнению
для переменной C1. Решив это уравнение и подставив

Вращение жидкости
11

найденное решение вместо C1 обратно в формулу общего 
решения, мы найдем решение соответствующей
начальной задачи для заданного дифференциального
уравнения. В нашем случае начальное условие имеет
вид y(0) = y0, поэтому мы сначала подставляем в общее 
решение dsol вместо x значение 0, это делается с
помощью метода subs(f, g), где f — выражение, которое 
мы хотим заменить на выражение g. Полученное
в результате подстановки выражение приравниваем к
величине y0, используя команду Eq(lhs, rhs), которая
создает равенство вида lhs = rhs6. Созданное уравне-
6
Python не позволяет переопределить 
команду присваивания =,
поэтому для создания уравнений
приходится
использовать
специальную 
функцию Eq.

ние запоминаем в переменной eq1.

1 eq1 = Eq(dsol.subs(x, 0), y0)

2 display(eq1)

C1 = y0

6 Следующим шагом решаем построенное уравнение
с помощью команды solveset(eq, x), где первый аргумент 
eq — это само уравнение, второй аргумент x
— переменная, относительно которой данное уравнение 
должно решаться7. Найденное решение сохраним
7 В рассматриваемом случае уравнение 
тривиально, но ровно такие
же действия нужно будет выполнять 
и в более сложных случаях.
в переменной sol1.

1 sol1 = solveset(eq1, C1)

2 display(sol1)

{y0}

7 Как мы видим, команда solveset выдает решения
уравнения в виде множества, чтобы выбрать какое-
то одно решение, нужно преобразовать это множество
в кортеж (tuple) или список (list) и выбрать соответствующий 
элемент, указав его индекс. В нашем случае
множество содержит всего одно решение, поэтому для
его извлечения используем нулевой индекс. Запоминаем 
найденное значение в переменной C2.

1 C2 = tuple(sol1)[0]

2 display(C2)

y0

8 Последним действием подставляем найденное значение 
С2 вместо символа C1 в общее решение dsol, что
и дает нам искомое решение начальной задачи. Сохраним 
это решение в переменной dsol_y0.

1 dsol_y0 = dsol.subs(C1, C2)

2 display(dsol_y0)

Глава 1

y0 + ω2x2

2g

9 Как видно из последней формулы, найденная зависимость 
высоты y поверхности жидкости от расстояния 
до оси вращения x оказывается квадратичной,
т. е. интегральные кривые (графики решений дифференциального 
уравнения) должны быть параболами8.
Убедимся в этом, построив графики решения. Библиотека 
SymPy имеет несколько простых функций для построения 
графиков9, простейшей из которых является 
функция plot. Построим с помощью этой функции
график найденного нами решения начальной задачи
при следующих значениях входящих в это решение
констант:
8 Таким образом, искомая поверхность 
будет парабалоидом вращения.

9
Являющихся,
по
сути,
надстройкой 
над существенно более
мощной графической библиотекой
Matplotlib.

y0 = 0, ω = π/2, g = 9.8.

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

1 ds = dsol_y0.subs({y0: 0, omega: pi/2, g: 9.8})

2 p1 = plot(ds)

РИС. 4 График решения начальной 
задачи

Результат выполнения10 этого фрагмента кода пока-

10 Переменная pi является предопределенной 
константой в SymPy.

зан на рис. 4.

10 Построим теперь графики семейства решений начальной 
задачи для различных значений угловой скорости 
ω ∈ [0, 2π] с шагом π/10. Для формирования
списка значений переменной omega будем использовать
функцию arange из библиотеки NumPy. Для построения
нескольких кривых на одном графике нужно сначала
создать пустой график и запомнить его в некоторой
переменной, в нашем случае это будет переменная p2.
Для каждой кривой следует создавать свой график
(переменная p) и затем добавлять его к общему графику 
p2 с помощью метода extend. По умолчанию любой
создаваемый график в среде Jupyter сразу же визуа-
лизируется, чтобы избежать этого, нужно в команде
plot использовать опцию show = False. Построенный
график можно визуализировать с помощью вызова метода 
show.

1 from numpy import arange

2 p2 = plot(show = False)

3 for om in arange(0, 2 * pi, pi / 10):