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

Данная статья объясняет дифференциальные уравнения без необходимости в знаниях в высшей математике

Вступление

Какой язык у природы? Физик наверняка сказал бы: дифференциальные уравнения (ДУ). Ими описано всё: от движения атомов до ураганов и орбит планет. Даже второй закон Ньютона – ДУ (ускорение – это вторая производная от положения). Существует больша́я проблема: большинство ДУ нельзя решить точно. Секреты природы скрыты за их сложностью. Значит ли это, что уравнения не работают? Нет! Просто точных решений нет. Но можно найти приближенные решения с помощью итераций и других численных методов.

Определение ОДУ

Для начала, разберём определение дифференциального уравнения.

Примеры: $y’(x) + y(x) + x = 0$, $f’’(x, y) + f(x, y) + y = 0$.
Введём определение обыкновенных дифференциальных уравнений, с которыми мы и будем работать.

Под функцией $F$ имеется в виду любое выражение, содержащее функции, их производные, и независимую переменную, от которой они зависят. Максимальная (по порядку) производная в уравнении определяет порядок ОДУ.
Решением дифференциального уравнения являются функции. Дифференциальное уравнение обычно имеет бесконечно много решений. Иногда, решение не является непрерывным на всём промежутке.
Важно уточнить зачем вообще нужны численные методы. Не каждое дифференциальное уравнение можно решить точно. Большинство из них решаются только численными методами, которые лишь приблизительно считают значение в некоторых точках.

Численные методы

Все численные методы приблизительно решают задачу Коши.

Вот задача Коши, которую мы будем решать в качестве примера:

$$ \begin{cases} 8e^{-x} \cdot (1 + x) -2y - y' = 0 \\ y(0) = 1 \end{cases} $$

Чтобы решить уравнение численным методом, нужно его представить в виде:

$$ y^{(n)} = \dots $$

То есть узнать значение максимальной производной. Переделаем задачу в следующую систему:

$$ \begin{cases} f(x, y) = 8e^{-x} \cdot (1 + x) -2y \\ y' = f(x, y) \\ y(0) = 1 \end{cases} $$

Таким образом, если мы знаем точку $(x, y)$, то мы можем посчитать в ней производную. Мы можем узнать значение производной в начальной точке: $y’(0) = f(0, 1) = 6$ (Можете убедится в этом сами, подставив числа в $f$). Это пригодится нам при применении численных методов.
Данная задача Коши имеет точное решение, получение которого не будет рассматриваться в данной статье, однако доказать данное решение легко (достаточно подставить функцию и её производную):

$$ y = e^{-2x} \cdot (8xe^x + 1) $$

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

Метод Эйлера(явный метод Рунге-Кутты первого порядка)

Для определения метода Эйлера, вспомним физическое определение производной:

А теперь посмотрим на формулы в физике, связанные с прямолинейным равномерным движением:

Тем не менее, в нашей задаче $y’ \neq \text{const}$, а поэтому мы не можем использовать формулу напрямую.
Пусть $h$ — некоторое число, назовём его шагом. Разделим нашу функцию на отдельные точки, разница по $x$-координате между каждой точкой будет $h$, начальная точка будет $(x_0, y_0)$.
Тогда при $h \to 0$:

$$ \begin{cases} x_i = x_{i-1} + h \\ y_i = y_{i-1} + h \cdot f(x_{i-1}, y_{i-1}) \end{cases} $$

Можно представить $h$, как шаг, мы ходим по точкам:

$$(x_0, y_0), \, (x_0 + h, y_1), \, (x_0+2h, y_2), \, \dots, \, (x_0 + n \cdot h, y_n)$$

и считаем значение в них. В этом и заключается метод Эйлера.
Давайте попробуем посчитать значения функции при помощи данного метода.
Пусть шаг $h$ равен $0.1$. Начальная точка у нас считается “нулевой” (по индексу), так что теперь считаем первую:

$$ \begin{cases} x_1= 0 + h = 0.1 \\ y_1 = 1 + h \cdot f(x_{0}, y_{0}) = 1 + 0.1 \cdot 6 = 1.6 \end{cases} $$

А теперь давайте посчитаем точное значение при $x = x_1$:

$$ y_1 = e^{-2\cdot 0.1} \cdot (8 \cdot {0.1} \cdot e^{0.1} + 1) \approx 1.5426 $$

Как видите, мы уже потеряли точность в $0.06$! Происходит это из-за относительно большого $h$. Тем не менее продолжаем:

$$ \begin{cases} x_2= x_1 + h = 0.2 \\ y_2 = y_1 + h \cdot f(x_{1}, y_{1}) \approx 1.6 + 0.1 \cdot 4.76 \approx 2.076 \end{cases} $$

Здесь к $f(x_1, y_1)$ применено округление. И снова посчитаем точное значение:

$$ y_1 = e^{-2\cdot 0.2} \cdot (8 \cdot {0.2} \cdot e^{0.2} + 1) \approx 1.98 $$

Теперь погрешность равна $0.09$.
Давайте посмотрим на график решения, и решения методом Эйлера (вычислено компьютером): Здесь хорошо видно, что метод Эйлера куда хуже справляется, если производная начинает быстро изменятся, впоследствии это будет заметно при применении других численных методов. При шаге в $0.01$ ошибка становится довольно незначительной.

Локальная и глобальная ошибка

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

Глобальная ошибка возникает из-за суммирования локальной ошибки. Чтобы сравнивать ошибки методов, нам нужно ввести обозначение “О-большое”.

Итак, локальная ошибка метода Эйлера равна $O(h^2)$, глобальная — $O(h)$. Доказательство данного факта требует знания рядов Тейлора.

Время работы

Хотелось бы также понимать, с какой скоростью компьютер вычисляет данный метод. Для этого введём понятия временной сложности алгоритма.

Правильное определение элементарной операции находится в машине Тьюринга, но мы для простоты будем под элементарной операцией понимать вычисление дифференциального уравнения (то есть $f$).
Заметим, что в каждой $i$-той точке мы вычисляем $f$ один раз. Определим количество вычисляемых точек как $n$, не включая начальную точку. Тогда временная сложность метода Эйлера равна $n = O(n)$.
Чтобы рассчитать время исполнения, нужно понять, сколько элементарных операций компьютер может сделать за секунду. В данной методичке мы возьмём скорость $2 \cdot 10^5$ операций в секунду. Значит, при помощи метода Эйлера мы можем обрабатывать $2\cdot 10^5$ точек в секунду.

Обратный метод Эйлера (неявный метод Рунге-Кутты первого порядка)

Модифицируем метод Эйлера следующим образом:

$$ \begin{cases} x_i = x_{i-1} + h \\ y_i = y_{i-1} + h \cdot f(x_{i}, y_{i}) \end{cases} $$

Возникает логичный вопрос, а как это решать? Как мы можем вычислить $y_i$, если нам для этого нужно знать $y_i$? Методы, где $y_i$ зависит от самого себя, называются неявными. Чтобы посчитать неявным методом, нужно превратить его в явный. Для этого приблизительно посчитаем $y_i$ методом Эйлера:

$$ \begin{cases} x_i = x_{i-1} + h \\ \tilde y_i = y_{i-1} + h \cdot f(x_{i-1}, y_{i-1}) \\ y_i = y_{i-1} + h \cdot f(x_{i}, \tilde y_{i}) \end{cases} $$

Знак тильды над $y_i$ обозначает предсказание. Мы предсказываем значение методом Эйлера и корректируем его методом Обратного Эйлера. Это является результатом обратного метода Эйлера через модель предиктора-корректора. Отметим, что вычисления можно повторять несколько раз для повышения точности.
Локальная ошибка обратного метода Эйлера равна $O(h^2)$, глобальная — $O(h)$. Зачем тогда им пользоваться? Я не указал ещё один способ сравнения методов: стабильность (иногда называется устойчивостью). Во многом потому, что данная тема довольна сложна. Обратный метод Эйлера куда более стабилен чем метод Эйлера. Разницу в стабильности можно увидеть на быстро сменяющихся функциях.

Метод Гейна (Модифицированный метод Эйлера, метод Рунге-Кутты второго порядка, RK2)

Можно ли соединить метод Эйлера и обратный метод Эйлера? Давайте просто возьмём из них среднее арифметическое:

$$ \begin{cases} x_i = x_{i-1} + h \\ \tilde y_i = y_{i-1} + h \cdot f(x_{i-1}, y_{i-1}) \\ y_i = y_{i-1} + h \cdot \frac {f(x_{i-1}, y_{i-1}) + f(x_{i}, \tilde y_{i})} {2} \end{cases} $$

Так мы получаем Метод Гейна.
Локальная ошибка метода Гейна равна $O(h^3)$, глобальная — $O(h^2)$.
Временная сложность алгоритма равна $2n = O(n)$.

Классический метод Рунге-Кутты четвёртого порядка (RK4)

Данный метод является наиболее используемым. Вот его формула:

$$ \begin{cases} x_i = x_{i-1} + h \\ y_i = y_{i-1} + \frac h 6 \cdot (k_1 + 2k_2 + 2k_3 + k_4) \end{cases} $$

Где:

$$ \begin{aligned} k_{1} &= f(x_{i-1},y_{i-1}),\\ k_{2} &= f\!\left(x_{i-1}+{\frac {h}{2}},y_{i-1}+h{\frac {k_{1}}{2}}\right),\\ k_{3} &= f\!\left(x_{i-1}+{\frac {h}{2}},y_{i-1}+h{\frac {k_{2}}{2}}\right),\\ k_{4} &= f\!\left(x_{i-1}+h,y_{i-1}+hk_{3}\right) \end{aligned} $$

Для того чтобы посчитать следующую точку, мы подсчитываем скорость в разных точках: $k_1$ — скорость в точке $x_{i-1}$
$k_2$ — скорость в середине
$k_3$ — скорректированная скорость в середине
$k_4$ — скорость в точке $x_i$
Далее мы берём среднее арифметическое с приоритетом (большим весом) на $k_2$ и $k_3$. Локальная ошибка RK4 равна $O(h^5)$, глобальная — $O(h^4)$.
Временная сложность: $4n = O(n)$

Система ОДУ

Далее все функции будут зависимы от $t$. Как и с обычными уравнениями, дифференциальные уравнение также можно составить в систему. Представим простейшую систему дифференциальных уравнений:

$$ \begin{cases} x'(t) &= f\left(t,\ x,\ y,\ y'\right) \\ y'(t) &= g\left( t,\ x,\ y,\ x' \right) \end{cases} $$

В таком случае, чтобы вычислить производную нам нужно знать ещё и решение в соседнем уравнении. Решение достаточно простое: давайте брать решения из прошлого! Вот решение методом Эйлера:

$$ \Large \begin{cases} t_i = t_{i-1} + h \\ x_i' = f(t_{i-1}, x_{i-1}, y_{i-1}, y_{i-1}') \\ x_i = x_{i-1} + h \cdot x_{i-1}' \\ y_i' = f(t_{i-1}, x_{i-1}, y_{i-1}, x_{i-1}' ) \\ y_i = y_{i-1} + h \cdot y_{i-1}' \end{cases} $$

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

ОДУ второго порядка

Все рассмотренные ранее уравнения были первого порядка, что же со вторым порядком?
Рассмотрим простейшее дифференциальное уравнение второго порядка:

$$ t - y(t) - y'(t) + y''(t) = 0 $$

Представим его в виде $y’’ = \dots$:

$$ y'' = y + y' - t $$

Сделаем замену $z = y’$:

$$ \begin{cases} y'(t) = z \\ z'(t) = y + z - t = f(t, y, z) \end{cases} $$

Таким образом получаем систему ОДУ первого порядка, которую мы уже можем решить.
Если избавится от замены, то вот как будет выглядеть решение методом Эйлера:

$$ \Large \begin{cases} t_i = t_{i-1} + h \\ y_i'(t) = y_{i-1}' + h \cdot f\left(t_{i-1}, y_{i-1}, y_{i-1}'\right) \\ y_i = y_{i-1} + h \cdot y_{i-1}' \end{cases} $$

Аналогично и с бо́льшими порядками.

Примеры решений

Решение рассматриваемой задачи Коши можно посмотреть на Google Colab.
Решение задачи про эластичный маятник также опубликовано на Colab.