Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений \begin{align} \tag{1} \frac{d u_i}{d t} &= f_i (t, u_1, u_2, \ldots, u_n), \quad t > 0\\ \tag{2} u_i(0) &= u_i^0, \quad i = 1, 2, \ldots, m. \end{align}
Используя векторные обозначения, задачу (1), (2) можно записать как задачу Коши \begin{align} \tag{3} \frac{d \pmb{u}}{d t} &= \pmb{F}(t, \pmb{u}), \quad t > 0, \\ \tag{4} \pmb{u}(0) &= \pmb{u}_0 \end{align} В задаче Коши необходимо по известному решению в точке t = 0 необходимо найти из уравнения (3) решение при других t .
Существует большое количество методов численного решения задачи (3), (4). Вначале рассмотрим простейший явный метод Эйлера и его программную реализацию. Затем будут представлены методы Рунге—Кутта и многошаговые методы.
При построении численных алгоритмов будем считать, что решение этой дифференциальной задачи существует, оно единственно и обладает необходимыми свойствами гладкости.
Идея численных методов решения задачи (3), (4) состоит из четырех частей:
1. Вводится расчетная сетка по переменной t (время) из N_t + 1 точки t_0 , t_1 , \ldots , t_{N_t} . Нужно найти значения неизвестной функции \pmb{u} в узлах сетки t_n . Обозначим через \pmb{y}^n приближенное значение \pmb{u}(t_n) .
2. Предполагаем, что дифференциальное уравнение выполнено в узлах сетки.
3. Аппроксимируем производные конечными разностями.
4. Формулируем алгоритм, который вычисляет новые значения \pmb{y}^{n+1} на основе предыдущих вычисленных значений \pmb{y}^k , k < n .
Метод сходится в точке t_n , если |\pmb{y}^n - \pmb{u}(t_n)| \to 0 при \tau \to 0 . Метод имеет p -ый порядок точности, если |\pmb{y}^n - \pmb{u}(t_n)| = O(\tau^p) , p > 0 при \tau\to 0 .
Проиллюстрируем указанные шаги. Для начала введем расчетную сетку. Очень часто сетка является равномерной, т.е. имеет одинаковое расстояние между узлами t_n и t_{n+1} : \omega_\tau = \{ t_n = n \tau, n = 0, 1, \ldots, N_t \}.
Затем, предполагаем, что уравнение выполнено в узлах сетки, т.е.: \pmb{u}^\prime (t_n) = \pmb{F}(t_n, u(t_n)), \quad t_n \in \omega_\tau.
Заменяем производные конечными разностями. С этой целью, нам нужно знать конкретные формулы, как производные могут быть аппроксимированы конечными разностями. Простейший подход заключается в использовании определения производной: \pmb{u}^\prime(t) = \lim_{\tau \to 0} \frac{\pmb{u}(t+\tau) - \pmb{u}(t)}{\tau}.
В произвольном узле сетки t_n это определение можно переписать в виде: \begin{equation*} \pmb{u}^\prime(t_n) = \lim_{\tau \to 0} \frac{\pmb{u}(t_n+\tau) - \pmb{u}(t_n)}{\tau}. \end{equation*} Вместо того, чтобы устремлять шаг сетки к нулю, мы можем использовать малый шаг \tau , который даст численное приближение u^\prime(t_n) : \begin{equation*} \pmb{u}^\prime(t_n) \approx \frac{\pmb{u}^{n+1} - \pmb{u}^{n}}{\tau}. \end{equation*} Такая аппроксимация известна как разностная производная вперед и имеет первый порядок по \tau , т.е. O(\tau) . Теперь можно использовать аппроксимацию производной. Таким образом получим явный метод Эйлера: \begin{equation} \tag{5} \frac{\pmb{y}^{n+1} - \pmb{y}^n}{\tau} = \pmb{F}(t_n, \pmb{y}^{n}). \end{equation}
Четвертый шаг заключается в получении численного алгоритма. Из (5) следует, что мы должны знать значение y^n для того, чтобы решить уравнение (5) относительно y^{n+1} и получить формулу для нахождения приближенного значения искомой функции на следующем временном слое t_{n+1} : \begin{equation} \tag{6} \pmb{y}^{n+1} = \pmb{y}^n + \tau \pmb{F}(t_n, \pmb{y}^{n}) \end{equation}
При условии, что у нас известно начальное значение \pmb{y}^0 = \pmb{u}_0 , мы можем использовать (6) для нахождения решений на последующих временных слоях.
Выражение (6) может быть как скалярным так и векторным уравнением. И в скалярном и в векторном случае на языке Python его можно реализовать следующим образом
u[n+1] = u[n] + tau*F_(t[n], u[n])
При решении системы (векторный случай), u[n] — одномерный массив
numpy длины m+1 ( m — размерность задачи), а функция F должна
возвращать numpy-массив размерности m+1 , t[n] — значение в
момент времени t_n .
Таким образом численное решение на отрезке [0, T] должно быть
представлено двумерным массивом, инициализируемым нулями
u = np.zeros((N_t+1, m+1)). Первый индекс соответствует временному
слою, а второй компоненте вектора решения на соответствующем временном
слое. Использование только одного индекса, u[n] или, что то же
самое, u[n, :], соответствует всем компонентам вектора решения.
Функция euler решения системы уравнений реализована в файле
euler.py:
import numpy as np
import matplotlib.pyplot as plt
def euler(F, u0, tau, T):
N_t = int(round(T/tau))
F_ = lambda t, u: np.asarray(F(t, u))
t = np.linspace(0, N_t*tau, N_t+1)
u = np.zeros((N_t+1, len(u0)))
u[0] = np.array(u0)
for n in range(N_t):
u[n+1] = u[n] + tau*F_(t[n], u[n])
return u, t
Строка F_ = lambda ... требует пояснений. Для пользователя,
решающего систему ОДУ, удобно задавать функцию правой части в виде
списка компонент. Можно, конечно, требовать чтобы пользователь
возвращал из функции массив numpy, но очень легко осуществлять
преобразование в самой функции решателе. Чтобы быть уверенным, что
результат F будет нужным массивом, который можно использовать в
векторных вычислениях, мы вводим новую функцию F_, которая вызывает
пользовательскую функцию F «прогоняет» результат через функцию
assaray модуля numpy.
При построении неявного метода Эйлера значение функции F берется на новом временном слое, т.е. для решении задачи (5) используется следующий метод: \begin{equation*} \tag{7} \frac{\pmb{y}^{n+1} - \pmb{y}^n}{\tau} = \pmb{F}(t_{n+1}, \pmb{y}^{n+1}). \end{equation*}
Таким образом для нахождения приближенного значения искомой функции на новом временном слое t_{n+1} нужно решить нелинейное уравнение относительно \pmb{y}^{n+1} : \begin{equation} \tag{8} \pmb{y}^{n+1} - \tau \pmb{F}(t_{n+1}, \pmb{y}^{n+1}) - y^n = 0. \end{equation}
Для решения уравнения (8) можно использовать, например, метод Ньютона.
Функция backward_euler решения системы уравнений реализована в файле
euler.py:
def backward_euler(F, u0, tau, T):
from scipy import optimize
N_t = int(round(T/tau))
F_ = lambda t, u: np.asarray(F(t, u))
t = np.linspace(0, N_t*tau, N_t+1)
u = np.zeros((N_t+1, len(u0)))
u[0] = np.array(u0)
def Phi(z, t, v):
return z - tau*F_(t, z) - v
for n in range(N_t):
u[n+1] = optimize.fsolve(Phi, u[n], args=(t[n], u[n]))
return u, t
Отметим, что для нахождения значения u[n+1] используется функция
fsolve модуля optimize библиотеки scipy. В качестве начального
приближения для решения нелинейного уравнения используется значение
искомой функции с предыдущего слоя u[n].
Одношаговый метод Рунге—Кутта в общем виде записывается следующим образом: \begin{equation} \tag{9} \frac{\pmb{y}^{n+1} - \pmb{y}^n}{\tau} = \sum_{i = 1}^s b_i \pmb{k}_i, \end{equation} где \begin{equation} \tag{10} \pmb{k}_i = \pmb{F}\left( t_n + c_i\tau, \pmb{y}^n + \tau \sum_{j=1}^s a_{ij}\pmb{k}_j \right), \quad i = 1, 2, \ldots, s. \end{equation} Формула (9) основана на s вычислениях функции \pmb{F} и называется s -стадийной. Если a_{ij} = 0 при j \geq i имеем явный метод Рунге—Кутта. Если a_{ij} = 0 при j > i и a_{ii} \ne 0 , то \pmb{k}_i определяется неявно из уравнения \begin{equation} \tag{11} \pmb{k}_i = \pmb{F}\left( t_n + c_i\tau, \pmb{y}^n + \tau \sum_{j=1}^{i-1} a_{ij}\pmb{k}_j + \tau a_{ii} \pmb{k}_i \right), \quad i = 1, 2, \ldots, s. \end{equation} О таком методе Рунге—Кутта говорят как о диагонально-неявном.
Одним из наиболее распространенных является явный метод Рунге-Кутта четвертого порядка: \begin{align*} \tag{12} \pmb{k}_1 & = \pmb{F}(t_n, \pmb{y}^n), &\quad \pmb{k}_2 &= \pmb{F}\left( t_n + \frac{\tau}{2}, \pmb{y}^n + \tau \frac{\pmb{k}_1}{2} \right),\\ \pmb{k}_3 &= \pmb{F}\left( t_n + \frac{\tau}{2}, \pmb{y}^n + \tau \frac{\pmb{k}_2}{2} \right), &\quad \pmb{k}_4 &= \pmb{F}\left( t_n + \tau, \pmb{y}^n + \tau \pmb{k}_3 \right),\\ \frac{\pmb{y}^{n+1} -\pmb{y}^n}{\tau} &= \frac{1}{6} (\pmb{k}_1 + 2\pmb{k}_2 + 2\pmb{k}_3 + \pmb{k}_4) & & \end{align*}
В методах Рунге—Кутта в вычислениях участвуют значения приближенного решения только в двух соседних узлах \pmb{y}^n и \pmb{y}^{n+1} — один шаг по переменной t . Линейный m -шаговый разностный метод записывается в виде \begin{equation} \tag{13} \frac{1}{\tau} \sum_{i=0}^m a_i \pmb{y}^{n+1-i} = \sum_{i=0}^{m} b_i \pmb{F}(t_{n+1-i}, \pmb{y}^{n+1-i}), \quad n = m-1, m, \ldots \end{equation} Вариант численного метода определяется заданием коэффициентов a_i , b_i , i = 0, 1, \ldots, m , причем a_0 \ne 0 . Для начала расчетов по рекуррентной формуле (13) необходимо задать m начальных значений \pmb{y}^0 , \pmb{y}^1 , \dots , \pmb{y}^{m-1} (например, можно использовать для их вычисления метод Эйлера).
Различные варианты многошаговых методов (методы Адамса) решения задачи с начальными условиями для систем обыкновенных дифференциальных уравнений могут быть получены на основе использования квадратурных формул для правой части равенства \begin{equation} \tag{14} \pmb{u}(t_{n+1}) - \pmb{u}(t_n) = \int_{t_n}^{t_{n+1}} \pmb{F}(t, \pmb{u}) dt \end{equation}
Для получения неявного многошагового метода используем для подынтегральной функции интерполяционную формулу по значениям функции \pmb{F}^{n+1} = \pmb{F}(t_{n+1}, \pmb{y}^{n+1}) , \pmb{F}^n , \dots , \pmb{F}^{n+1-m} , т.е. \begin{equation} \tag{15} \frac{\pmb{y}^{n+1} - \pmb{y}^n}{\tau} = \sum_{i=0}^{m} b_i \pmb{F}(t_{n+1-i}, \pmb{y}^{n+1-i}) \end{equation}
Для интерполяционного метода Адамса (15) наивысший порядок аппроксимации равен m+1 .
Для построения явных многошаговых методов можно использовать процедуру экстраполяции подынтегральной функции в правой части (14). В этом случае приближение осуществляется по значениям \pmb{F}^n , \pmb{F}^{n-1} , \dots , \pmb{F}^{n+1-m} и поэтому \begin{equation} \tag{16} \frac{\pmb{y}^{n+1} - \pmb{y}^n}{\tau} = \sum_{i=1}^{m} b_i \pmb{F}(t_{n+1-i}, \pmb{y}^{n+1-i}) \end{equation}
Для экстраполяционного метода Адамса (16) погрешность аппроксимации имеет m -ый порядок.
Примерами методов Адамса (15), (16) при m=3 являются \begin{align} \tag{17} \frac{\pmb{y}^{n+1} - \pmb{y}^n}{\tau} &= \frac{1}{24} (9\pmb{F}^{n+1} + 19\pmb{F}^{n} - 5\pmb{F}^{n-1} + \pmb{F}^{n-2}),\\ \tag{18} \frac{\pmb{y}^{n+1} - \pmb{y}^n}{\tau} &= \frac{1}{12} (23 \pmb{F}^{n} -16\pmb{F}^{n-1} + 5\pmb{F}^{n-2}), \end{align} соответственно.
На основе методов Адамса строятся и схемы предиктор–корректор. На этапе предиктор используется явный метод Адамса, на этапе корректора — аналог неявного метода Адамса. Например, при использовании методов третьего порядка аппроксимации в соответствии с (18) для предсказания решения положим \frac{\pmb{y}^{n+1} - \pmb{y}^n}{\tau} = \frac{1}{12} (23 \pmb{F}^{n} -16\pmb{F}^{n-1} + 5\pmb{F}^{n-2}). Для уточнеия решения (см. (17)) используется схема \frac{\pmb{y}^{n+1} - \pmb{y}^n}{\tau} = \frac{1}{24} (9\pmb{F}^{n+1} + 19\pmb{F}^{n} - 5\pmb{F}^{n-1} + \pmb{F}^{n-2}). Аналогично строятся и другие классы многошаговых методов.
При численном решении задачи Коши для систем обыкновенных дифференциальных уравнений (3), (4) могут возникнуть дополнительные трудности, порожденные жесткостью системы. Локальные особенности поведения решения в точке u = w передаются линейной системой \begin{equation*} \frac{d v_i}{dt} = \sum_{j=1}^{m} \frac{\partial f_i}{\partial u_j} (t, w) v + \bar{f_i}(t), \quad t > 0. \end{equation*}
Пусть \lambda_i(t) , i = 1, 2, \ldots, m — собственные числа матрицы \begin{equation*} A(t) = \{ a_{ij}(t) \}, \quad a_{ij}(t) = \frac{\partial f_i}{\partial u_j}(t, w). \end{equation*} Система уравнений (3) является жесткой, если число \begin{equation*} S(t) = \frac{\max_{1 \leq i \leq m} |Re \lambda_i(t)|}{\min_{1 \leq i \leq m} |Re \lambda_i(t)|} \end{equation*} велико. Это означает, что в решении присутствуют составляющие с сильно различающимися масштабами изменения по переменной t .
Для численное решения жестких задач используются вычислительные алгоритмы, которые имеют повышенный запас устойчивости. Необходимо ориентироваться на использование A -устойчивых или A(\alpha) -устойчивых методов.
Метод называется A -устойчивым, если при решении задачи Коши для системы (3) область его устойчивости содержит угол \begin{equation*} |\arg(-\mu)| < \alpha, \quad \mu = \lambda \tau. \end{equation*}
Среди A -устойчивых методов можно выделить чисто неявные многошаговые методы (методы Гира), когда \begin{equation} \tag{19} \frac{1}{\tau} \sum_{i=0}^{k} a_i \pmb{y}^{n+1-i} = \pmb{F}(t_n, \pmb{y}^{n+1}) \end{equation}
Написать программу для численного решения задачи Коши для системы обыкновенных дифференциальных уравнений явным методом Рунге—Кутта четвертого порядка. Продемонстрировать работоспособность этой программы при решении задачи Коши (построить график зависимости решения от t ) \frac{d^2 u}{d t^2} = -\sin u, \quad 0 < t < 4\pi, u(0) = 1, \quad \frac{d u}{d t} = 0.
Написать программу, которая реализует метод предиктор–корректор: \begin{align*} \frac{\bar{\pmb{y}}^{n+1} - \pmb{y}^n}{\tau} &= \frac{1}{24} (55 \pmb{F}(t_n, \pmb{y}^{n}) -59\pmb{F}(t_{n-1}, \pmb{y}^{n-1}) + 37 \pmb{F}(t_{n-2}, \pmb{y}^{n-2}) - 9 \pmb{F}(t_{n-3}, \pmb{y}^{n-3})), \\ \frac{\pmb{y}^{n+1} - \pmb{y}^n}{\tau} &= \frac{1}{24} (9\pmb{F}(t_{n+1}, \bar{\pmb{y}}^{n+1}) + 19\pmb{F}(t_n, \pmb{y}^{n}) - 5\pmb{F}(t_{n-1}, \pmb{y}^{n-1}) + \pmb{F}(t_{n-2}, \pmb{y}^{n-2}). \end{align*} С помощью этой программы решить задачу \begin{align*} \frac{d u_1}{dt} &= u_2 - u_3, \\ \frac{d u_2}{dt} &= u_1 + a u_2,\\ \frac{d u_3}{dt} &= b + u_3(u_1 - c), \quad 0 < t \leq 100,\\ u_1(0) &= 1, \quad u_2(0) = 1, \quad u_3(0) = 1. \end{align*} при a, b, c = 0.2, 0.2, 2.5 и a, b, c = 0.2, 0.2, 5 .
Написать программу для приближенного решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка методом Гира четвертого порядка точности: \begin{equation*} \frac{25 \pmb{y}^{n+1} - 48 \pmb{y}^n + 36 \pmb{y}^{n-1} - 16 \pmb{y}^{n-2} + 3\pmb{y}^{n-3}}{12\tau} = \pmb{F}(t_{n+1}, \pmb{y}^{n+1}) \end{equation*} С ее помощью найти решение задачи \begin{align*} \frac{du_1}{dt} &= u_2, \\ \frac{d u_2}{dt} &= \mu(1-u_1^2)u_2 - u_1, \quad 0 < t \leq 100,\\ u_1(0) &= 2, \quad u_2(0) = 0. \end{align*} при \mu = 50 .