Неявный метод эйлера для решения дифференциальных уравнений

Неявные методы решения дифференциальных уравнений

1. Неявные методы

Методы численного решения ОДУ, с которыми мы познакомились в первом разделе этого курса (метод Эйлера, метод средней точки и т. п.), называются «явными» методами. Однако иногда система ОДУ может стать «жесткой», а решать такие системы явными методами неудобно. В этом случае желательно изменить формулировку задачи так, что не пришлось иметь дело с жесткой системой. Но это не всегда возможно, поэтому вы должны уметь решать жесткие ОДУ. Для этого, как правило, используют методы решения, которые называются «неявными».

2. Пример жесткого ОДУ

Во-первых, в чем смысл и причина появления жестких уравнений? Давайте рассмотрим пример, который часто возникает в динамике. Предположим, что у нас есть частица, с координатами \((x(t),y(t))\) , и предположим, что мы хотим, чтобы ее \(y\) -координата всегда оставалась равной нулю. Один из способов добиться этого — добавить слагаемое \(-ky(t)\) , к производной \(\dot(t)\) , где \(k\) — большая положительная постоянная. Если \(k\) достаточно велико, то частица никогда не уйдет далеко от \(y(t)=0\) , так как слагаемое \(-ky(t)\) всегда приведет \(y(t)\) обратно к нулю. Предположим далее, что мы хотим, чтобы пользователь мог перемещать частицу как угодно вдоль оси \(x\) . Дифференциальное уравнение, описывающее движение данной системы, будет иметь вид

(Кроме того, мы предполагаем, что частица не запускается из \(y_0=0\) ). В результате частица будет сильно притягиваться к прямой \(y = 0\) , и менее сильно — к \(x = 0\) . Если решать ОДУ на достаточно продолжительном интервале времени, то частица рано или поздно попадет в точку \((0, 0)\) и останется в ней.

Теперь предположим, что для решения уравнения мы используем метод Эйлера. Если сделать шаг размера \(h\) , то получим

Если мы посмотрим на \(у\) -компоненту этого уравнения, то увидим, что при \(|1-hk|>1\) , вычисленное нами \(y_\) будет по модулю больше, чем \(|у_0|\) . Другими словами, при \(|1-hk|>1\) метод Эйлера будет неустойчив: каждый шаг приводит к увеличению \(y_\) по сравнению с предыдущим значением и приближенное решение будут все дальше отклоняться от нуля. Таким образом, для обеспечения устойчивости метода нужно, чтобы \(1-hk>-1\) или \(hk . Самый большой шаг, который мы можем сделать не нарушив устойчивости, должен быть меньше \(2/k\) .

Теперь, если \(k\) – большое число, мы вынуждены делать очень маленькие шаги. Это означает, что частица будет двигаться к \((0,0)\) мучительно медленно. Даже если взять \(y_0\) очень близким к нулю, то придется делать настолько маленькие шаги, что изменение \(x\) -координаты будет практически незаметно. Вот так выглядит жесткое ОДУ. В данном случае жесткость возникает из-за слишком большого \(k\) , призванного удержать частицу возле прямой \(у = 0\) . Позже, когда мы будем рассматривать частицы, соединенные пружинами, мы увидим то же самое явление: от жесткости пружины и происходит термин «жесткое» ОДУ. Даже если мы используем более совершенный численный метод, такой как метод Рунге-Кутта 4-го порядка, это лишь слегка улучшит ситуацию с выбором величины шага, но мы все равно будем иметь серьезные вычислительные проблемы.

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

3. Решение жесткого ОДУ

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

Формула явного метода Эйлера

продвигает систему вперед на шаг \(h\) во времени. Для жестких систем, однако, удобнее заменить эту формулу на следующую

То есть, нам нужно вычислить \(f\) в точке, в которую мы стремимся попасть ( \(\textbf_\) ), а не в исходной ( \(\textbf_0\) ). (Если представить, что время может двигаться в обратном направление, то смысл этого уравнения очевиден. Оно говорит: «если вы находитесь в \(\textbf_\) , и сделаете шаг \(-hf (\textbf_)\) , то попадете в \(\textbf_0\) ». Так что если ваше ДУ представляет собой систему, движущуюся вспять во времени, то этот шаг имеет смысл. Это просто поиск точки \(\textbf_\) такой, что если запустить время вспять, вы в конечном итоге окажетесь в \(\textbf_0\) .) Таким образом, мы ищем точку \(\textbf_\) такую, что \(f\) , вычисленная в ней и умноженная на \(h\) , приводит к исходной точке \(\textbf_0\) . К сожалению, найти \(\textbf_\) из уравнения \eqref в общем случае невозможно, если только \(f\) не является линейной функцией.

Чтобы справиться с этим, заменим \(f (\textbf_)\) линейной аппроксимацией, основанной на разложении \(f\) в ряд Тейлора. Введем обозначение \(\Delta\textbf=\textbf_-\textbf_0\) . Подставив его в уравнение \eqref, получим

Теперь заменим \(f(\textbf_0 + \Delta\textbf)\) следующим приближением

(Заметим, что поскольку \(f(\textbf_0)\) является вектором, то производная \(f^\prime(\textbf_0)\) является матрицей.) Используя это приближение, мы можем записать \(\Delta\textbf\) как

Разделим обе части последнего соотношения на \(h\) и перепишем результат в виде

где \(I\) — единичная матрица.

Разрешая это соотношение относительно \(\Delta\textbf\) , получим

Вычисление \(\textbf_=\textbf_0+\Delta\textbf\) для неявного метода очевидно требует больших вычислительных затрат, чем при использовании явного метода, так как мы должны решать систему линейных уравнений на каждом шаге. Хотя это может показаться серьезным недостатком (в вычислительном плане), не отчаивайтесь (пока). Для многих классов задач, матрица \(f^\prime\) будет разреженной — например, для «решетки» частиц, соединенных пружинами, \(f^\prime\) будет иметь структуру, которая соответствует связям между частицами. (Обсуждение разреженности и методов решений, применимых в этом случае, см. Baraff и Witkin [1]. Основной материал в Press et al. [2] также будет полезен.) В результате, как правило, можно решить уравнение \eqref в линейном времени (т. е. за время, пропорциональное размерности \(\textbf\) ). Выигрыш в таких случаях будет весьма существенным: мы, как правило, можем делать большие шаги по времени без потери устойчивости (т. е. без расхождения, как это происходит для явного метода, если длина шага слишком велика). Время, необходимое для решения каждой линейной системы, таким образом, более чем компенсируется тем, шагов можно зачастую сделать на порядки больше, чем при использовании явных методов. (Конечно, код, необходимый для реализации всего этого, гораздо сложнее, чем в явном случае; как мы уже говорили: переделывайте свои задачи в нежесткие, а если не получается, то платите положенную цену.)

Применим теперь неявный метод для решения уравнения \eqref. В нашем случае \(f(\textbf(t))\) равно

Дифференцирование по \(\textbf\) дает

Тогда матрица \((1/h)\textbf — f^<\prime>(\textbf_0)\) будет равна

Обращая эту матрицу и умножая на \(f(\textbf_0)\) , получим

Какова же предельная длина шага в этом случае? Ответ: нет предела! Если позволить \(h\) расти до бесконечности, мы получаем следующее

Это означает, что мы достигнем \(\textbf_=\textbf_0 + (-\textbf_0)=0\) за один шаг! В общем случае для жесткого ОДУ мы не сможем сделать шаг произвольного размера, но мы сможем сделать его гораздо большим, чем если бы использовали явный метод. Дополнительные расходы на решение системы линейных уравнений компенсируются экономией, возникающий благодаря возможности сделать меньше шагов.

4. Решение уравнений второго порядка

Большинство задач динамики записывается в виде ДУ 2-го порядка

Это уравнение легко преобразуется систему ДУ 1-го порядка, добавлением новых переменных. Если мы определим \(\textbf=\dot<\textbf>\) , то сможем переписать уравнение \eqref в виде

что представляет собой систему ДУ 1-го порядка. Однако, применяя обратный (неявный) метод Эйлера к уравнению \eqref, получим линейную систему размера \(2n \times 2n\) где \(n\) — размерность \(\textbf\) . Простое преобразование позволяет уменьшить размер линейной системы до \(n \times n\) . Важно отметить, что обе системы — исходная и преобразованная — будут иметь одинаковую степень разреженности. Таким образом, решение системы меньшего размера будет выполняться быстрее.

Система \(n \times n\) , которая должна быть решена, получается следующим образом. Введем для краткости следующие обозначения \(\textbf_0=\textbf(t_0)\) и \(\textbf_0=\textbf(t_0)\) . Определяем \(\Delta\textbf\) и \(\Delta\textbf\) как \(\Delta\textbf = \textbf(t_0+h)-\textbf(t_0)\) и \(\Delta\textbf = \textbf(t_0+h)-\textbf(t_0)\) . Очередное приближение по неявному методу Эйлера, примененному к уравнению \eqref, дает в результате

Применяя к \(f\) разложение в ряд Тейлора, которое в данном контексте является функцией и \(\textbf\) и \(\textbf\) , получим приближение 1-го порядка

В этом уравнении, производная \(\partial f / \partial\textbf\) оценивается в точке \((\textbf_0, \textbf_0)\) и аналогично вычисляется \(\partial f / \partial\textbf\) . Подставляя это приближение в уравнение \eqref, получим линейную систему

Подставив во вторую строку последнего равенства соотношение \(\Delta\textbf = h(\textbf_0 + \Delta\textbf)\) (т. е. первую строку того же равенства), придем к

Вводя единичную матрицу \(I\) и перегруппировав члены, получим соотношение

из которого находим \(\Delta\textbf\) . А зная \(\Delta\textbf\) легко вычислить \(\Delta\textbf = h(\textbf_0 + \Delta\textbf)\) .

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

Ссылки

1. D. Baraff and A. Witkin. Large steps in cloth simulation. Computer Graphics (Proc. SIGGRAPH), 1998.
2. W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling. Numerical Recipes. Cambridge University Press, 1986.

Читайте также

Комментарии

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

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

Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений $$ \begin \tag <1>\frac &= 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 $$

Используя векторные обозначения, задачу (1), (2) можно записать как задачу Коши $$ \begin \tag <3>\frac> &= \pmb(t, \pmb), \quad t > 0, \\ \tag <4>\pmb(0) &= \pmb_0 \end $$ В задаче Коши необходимо по известному решению в точке \( t = 0 \) необходимо найти из уравнения (3) решение при других \( t \).

Численные методы решения задачи Коши

Существует большое количество методов численного решения задачи (3), (4). Вначале рассмотрим простейший явный метод Эйлера и его программную реализацию. Затем будут представлены методы Рунге—Кутта и многошаговые методы.

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

Идея численных методов решения задачи (3), (4) состоит из четырех частей:

1. Вводится расчетная сетка по переменной \( t \) (время) из \( N_t + 1 \) точки \( t_0 \), \( t_1 \), \( \ldots \), \( t_ \). Нужно найти значения неизвестной функции \( \pmb \) в узлах сетки \( t_n \). Обозначим через \( \pmb^n \) приближенное значение \( \pmb(t_n) \).

2. Предполагаем, что дифференциальное уравнение выполнено в узлах сетки.

3. Аппроксимируем производные конечными разностями.

4. Формулируем алгоритм, который вычисляет новые значения \( \pmb^ \) на основе предыдущих вычисленных значений \( \pmb^k \), \( k 0 \) при \( \tau\to 0 \).

Явный метод Эйлера

Проиллюстрируем указанные шаги. Для начала введем расчетную сетку. Очень часто сетка является равномерной, т.е. имеет одинаковое расстояние между узлами \( t_n \) и \( t_ \): $$ \omega_\tau = \< t_n = n \tau, n = 0, 1, \ldots, N_t \>. $$

Затем, предполагаем, что уравнение выполнено в узлах сетки, т.е.: $$ \pmb^\prime (t_n) = \pmb(t_n, u(t_n)), \quad t_n \in \omega_\tau. $$

Заменяем производные конечными разностями. С этой целью, нам нужно знать конкретные формулы, как производные могут быть аппроксимированы конечными разностями. Простейший подход заключается в использовании определения производной: $$ \pmb^\prime(t) = \lim_ <\tau \to 0>\frac<\pmb(t+\tau) — \pmb(t)><\tau>. $$

В произвольном узле сетки \( t_n \) это определение можно переписать в виде: $$ \begin \pmb^\prime(t_n) = \lim_ <\tau \to 0>\frac<\pmb(t_n+\tau) — \pmb(t_n)><\tau>. \end $$ Вместо того, чтобы устремлять шаг сетки к нулю, мы можем использовать малый шаг \( \tau \), который даст численное приближение \( u^\prime(t_n) \): $$ \begin \pmb^\prime(t_n) \approx \frac<\pmb^ — \pmb^><\tau>. \end $$ Такая аппроксимация известна как разностная производная вперед и имеет первый порядок по \( \tau \), т.е. \( O(\tau) \). Теперь можно использовать аппроксимацию производной. Таким образом получим явный метод Эйлера: $$ \begin \tag <5>\frac<\pmb^ — \pmb^n> <\tau>= \pmb(t_n, \pmb^). \end $$

Четвертый шаг заключается в получении численного алгоритма. Из (5) следует, что мы должны знать значение \( y^n \) для того, чтобы решить уравнение (5) относительно \( y^ \) и получить формулу для нахождения приближенного значения искомой функции на следующем временном слое \( t_ \): $$ \begin \tag <6>\pmb^ = \pmb^n + \tau \pmb(t_n, \pmb^) \end $$

При условии, что у нас известно начальное значение \( \pmb^0 = \pmb_0 \), мы можем использовать (6) для нахождения решений на последующих временных слоях.

Программная реализация явного метода Эйлера

Выражение (6) может быть как скалярным так и векторным уравнением. И в скалярном и в векторном случае на языке Python его можно реализовать следующим образом

При решении системы (векторный случай), 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:

Строка F_ = lambda . требует пояснений. Для пользователя, решающего систему ОДУ, удобно задавать функцию правой части в виде списка компонент. Можно, конечно, требовать чтобы пользователь возвращал из функции массив numpy , но очень легко осуществлять преобразование в самой функции решателе. Чтобы быть уверенным, что результат F будет нужным массивом, который можно использовать в векторных вычислениях, мы вводим новую функцию F_ , которая вызывает пользовательскую функцию F «прогоняет» результат через функцию assaray модуля numpy .

Неявный метод Эйлера

При построении неявного метода Эйлера значение функции \( F \) берется на новом временном слое, т.е. для решении задачи (5) используется следующий метод: $$ \begin \tag <7>\frac<\pmb^ — \pmb^n> <\tau>= \pmb(t_, \pmb^). \end $$

Таким образом для нахождения приближенного значения искомой функции на новом временном слое \( t_ \) нужно решить нелинейное уравнение относительно \( \pmb^ \): $$ \begin \tag <8>\pmb^ — \tau \pmb(t_, \pmb^) — y^n = 0. \end $$

Для решения уравнения (8) можно использовать, например, метод Ньютона.

Программная реализация неявного метода Эйлера

Функция backward_euler решения системы уравнений реализована в файле euler.py:

Отметим, что для нахождения значения u[n+1] используется функция fsolve модуля optimize библиотеки scipy . В качестве начального приближения для решения нелинейного уравнения используется значение искомой функции с предыдущего слоя u[n] .

Методы Рунге—Кутта

Одношаговый метод Рунге—Кутта в общем виде записывается следующим образом: $$ \begin \tag <9>\frac<\pmb^ — \pmb^n> <\tau>= \sum_^s b_i \pmb_i, \end $$ где $$ \begin \tag <10>\pmb_i = \pmb\left( t_n + c_i\tau, \pmb^n + \tau \sum_^s a_\pmb_j \right), \quad i = 1, 2, \ldots, s. \end $$ Формула (9) основана на \( s \) вычислениях функции \( \pmb \) и называется \( s \)-стадийной. Если \( a_ = 0 \) при \( j \geq i \) имеем явный метод Рунге—Кутта. Если \( a_ = 0 \) при \( j > i \) и \( a_ \ne 0 \), то \( \pmb_i \) определяется неявно из уравнения $$ \begin \tag <11>\pmb_i = \pmb\left( t_n + c_i\tau, \pmb^n + \tau \sum_^ a_\pmb_j + \tau a_ \pmb_i \right), \quad i = 1, 2, \ldots, s. \end $$ О таком методе Рунге—Кутта говорят как о диагонально-неявном.

Одним из наиболее распространенных является явный метод Рунге-Кутта четвертого порядка: $$ \begin \tag <12>\pmb_1 & = \pmb(t_n, \pmb^n), &\quad \pmb_2 &= \pmb\left( t_n + \frac<\tau><2>, \pmb^n + \tau \frac<\pmb_1> <2>\right),\\ \pmb_3 &= \pmb\left( t_n + \frac<\tau><2>, \pmb^n + \tau \frac<\pmb_2> <2>\right), &\quad \pmb_4 &= \pmb\left( t_n + \tau, \pmb^n + \tau \pmb_3 \right),\\ \frac<\pmb^ -\pmb^n> <\tau>&= \frac<1> <6>(\pmb_1 + 2\pmb_2 + 2\pmb_3 + \pmb_4) & & \end $$

Многошаговые методы

В методах Рунге—Кутта в вычислениях участвуют значения приближенного решения только в двух соседних узлах \( \pmb^n \) и \( \pmb^ \) — один шаг по переменной \( t \). Линейный \( m \)-шаговый разностный метод записывается в виде $$ \begin \tag <13>\frac<1> <\tau>\sum_^m a_i \pmb^ = \sum_^ b_i \pmb(t_, \pmb^), \quad n = m-1, m, \ldots \end $$ Вариант численного метода определяется заданием коэффициентов \( a_i \), \( b_i \), \( i = 0, 1, \ldots, m \), причем \( a_0 \ne 0 \). Для начала расчетов по рекуррентной формуле (13) необходимо задать \( m \) начальных значений \( \pmb^0 \), \( \pmb^1 \), \( \dots \), \( \pmb^ \) (например, можно использовать для их вычисления метод Эйлера).

Различные варианты многошаговых методов (методы Адамса) решения задачи с начальными условиями для систем обыкновенных дифференциальных уравнений могут быть получены на основе использования квадратурных формул для правой части равенства $$ \begin \tag <14>\pmb(t_) — \pmb(t_n) = \int_^> \pmb(t, \pmb) dt \end $$

Для получения неявного многошагового метода используем для подынтегральной функции интерполяционную формулу по значениям функции \( \pmb^ = \pmb(t_, \pmb^) \), \( \pmb^n \), \( \dots \), \( \pmb^ \), т.е. $$ \begin \tag <15>\frac<\pmb^ — \pmb^n> <\tau>= \sum_^ b_i \pmb(t_, \pmb^) \end $$

Для интерполяционного метода Адамса (15) наивысший порядок аппроксимации равен \( m+1 \).

Для построения явных многошаговых методов можно использовать процедуру экстраполяции подынтегральной функции в правой части (14). В этом случае приближение осуществляется по значениям \( \pmb^n \), \( \pmb^ \), \( \dots \), \( \pmb^ \) и поэтому $$ \begin \tag <16>\frac<\pmb^ — \pmb^n> <\tau>= \sum_^ b_i \pmb(t_, \pmb^) \end $$

Для экстраполяционного метода Адамса (16) погрешность аппроксимации имеет \( m \)-ый порядок.

На основе методов Адамса строятся и схемы предиктор–корректор. На этапе предиктор используется явный метод Адамса, на этапе корректора — аналог неявного метода Адамса. Например, при использовании методов третьего порядка аппроксимации в соответствии с (18) для предсказания решения положим $$ \frac<\pmb^ — \pmb^n> <\tau>= \frac<1> <12>(23 \pmb^ -16\pmb^ + 5\pmb^). $$ Для уточнеия решения (см. (17)) используется схема $$ \frac<\pmb^ — \pmb^n> <\tau>= \frac<1> <24>(9\pmb^ + 19\pmb^ — 5\pmb^ + \pmb^). $$ Аналогично строятся и другие классы многошаговых методов.

Жесткие системы ОДУ

При численном решении задачи Коши для систем обыкновенных дифференциальных уравнений (3), (4) могут возникнуть дополнительные трудности, порожденные жесткостью системы. Локальные особенности поведения решения в точке \( u = w \) передаются линейной системой $$ \begin \frac