Метод интегрирования дифференциальных уравнений эйлера коши
Дифференциальное уравнение Эйлера и методы его решения
Более общее уравнение Эйлера имеет вид: . Это уравнение подстановкой t = ax+b приводится к более простому виду, которое мы и будем рассматривать.
Приведение дифференциального уравнения Эйлера к уравнению с постоянными коэффициентами.
Рассмотрим уравнение Эйлера: (1) . Оно сводится к линейному уравнению с постоянными коэффициентами подстановкой: x = e t . Действительно, тогда ; ; ;
; ; .
Таким образом, множители, содержащие x m , сокращаются. Остаются члены с постоянными коэффициентами. Однако на практике, для решения уравнений Эйлера, можно применять методы решения линейных ДУ с постоянными коэффициентами без использования указанной выше подстановки.
Решение однородного уравнения Эйлера
Рассмотрим однородное уравнение Эйлера: (2) . Ищем решение уравнения (2) в виде . ; ; . . Подставляем в (2) и сокращаем на x k . Получаем характеристическое уравнение: . Решаем его и получаем n корней, которые могут быть комплексными.
Рассмотрим действительные корни. Пусть ki – кратный корень кратности m . Этим m корням соответствуют m линейно независимых решений: .
Рассмотрим комплексные корни. Они появляются парами вместе с комплексно сопряженными. Пусть ki – кратный корень кратности m . Выразим комплексный корень ki через действительную и мнимую части: . Этим m корням и m комплексно сопряженным корням соответствуют 2 m линейно независимых решений: ; ; . .
После того как получены n линейно независимых решений, получаем общее решение уравнения (2): (3) .
Примеры
Решение неоднородного уравнения Эйлера
Рассмотрим неоднородное уравнение Эйлера: . Метод вариации постоянных (метод Лагранжа) также применим и к уравнениям Эйлера.
Сначала мы решаем однородное уравнение (2) и получаем его общее решение (3). Затем считаем постоянные функциями от переменной x . Дифференцируем (3) n – 1 раз. Получаем выражения для n – 1 производных y по x . При каждом дифференцировании члены, содержащие производные приравниваем к нулю. Так получаем n – 1 уравнений, связывающих производные . Далее находим n -ю производную y . Подставляем полученные производные в (1) и получаем n -е уравнение, связывающее производные . Из этих уравнений определяем . После чего интегрируя, получаем общее решение уравнения (1).
Пример
Неоднородное уравнение Эйлера со специальной неоднородной частью
Рассмотрим уравнение Эйлера со специальной неоднородной частью: (4) , где – многочлены от степеней и , соответственно.
Наиболее простой способ решения такого уравнения заключается в том, чтобы сделать подстановку , и решать линейное уравнение с постоянными коэффициентами со специальной неоднородной частью.
Автор: Олег Одинцов . Опубликовано: 14-08-2013 Изменено: 24-10-2020
Задачи с начальными условиями для систем обыкновенных дифференциальных уравнений
Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений $$ \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. Предполагаем, что дифференциальное уравнение выполнено в узлах сетки.
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 $$ О таком методе Рунге—Кутта говорят как о диагонально-неявном.
В методах Рунге—Кутта в вычислениях участвуют значения приближенного решения только в двух соседних узлах \( \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
= \sum_^ \frac<\partial f_i> <\partial u_j>(t, w) v + \bar(t), \quad t > 0. \end $$
Пусть \( \lambda_i(t) \), \( i = 1, 2, \ldots, m \) — собственные числа матрицы $$ \begin A(t) = \< a_(t) \>, \quad a_(t) = \frac<\partial f_i><\partial u_j>(t, w). \end $$ Система уравнений (3) является жесткой, если число $$ \begin S(t) = \frac <\max_<1 \leq i \leq m>|Re \lambda_i(t)|> <\min_<1 \leq i \leq m>|Re \lambda_i(t)|> \end $$ велико. Это означает, что в решении присутствуют составляющие с сильно различающимися масштабами изменения по переменной \( t \).
Для численное решения жестких задач используются вычислительные алгоритмы, которые имеют повышенный запас устойчивости. Необходимо ориентироваться на использование \( A \)-устойчивых или \( A(\alpha) \)-устойчивых методов.
Метод называется \( A \)-устойчивым, если при решении задачи Коши для системы (3) область его устойчивости содержит угол $$ \begin |\arg(-\mu)| —>
Метод интегрирования дифференциальных уравнений эйлера коши
последовательностью точек x0,x1. в соответствующие моменты времени t0,t1. Значения точек должны удоволетворять приближенному равенству
Если специально не оговорено иное, то предполагается, что моменты времени выбираются через равные интервалы с величиной шага h>0, то есть
Аппроксимируем производную в момент времени tk соотношением
При такой аппроксимации уравнение (1) примет вид:
Формула (2) известна как прямой метод Эйлера.
На рис.1(a) показана графическая интерпретация прямого метода Эйлера. На (k+1)-ом шаге векторное поле предполагается (локально) постоянным со значением f(xk,tk).
Рис.1 Иллюстрация алгоритмов (а) прямого метода Эйлера, (b) обратного метода Эйлера
Меньшее значение величины шага h в итоге дает точки аппроксимации чаще и, как демонстрирует рис.2, приводит к большей точности интегрирования, что приобретает математический смысл, поскольку (2) стремится к (1) при h->0.
Рис.2 Влияние величины шага. Уравнение dx/dt=-6x+5t -t интегрируется от x=1 прямым методом Эйлера при h=0.3 (а) и при h=0.1 (b). Точное решение показано штриховой линией.
Обратный метод Эйлера подобен прямому, но есть одно отличие в аппроксимации для производной
Такая аппроксимация дает формулу обратного метода Эйлера:
На рис.1(b) показана геометрическая интерпретация обратного метода Эйлера. На (k+1)-ом шаге векторное поле предполагается (локально) постоянным со значением f(xk+1,tk+1).
Обратный метод Эйлера — это пример неявного алгоритма интегрирования , где xk+1 является функцией от самой себя. И напротив, прямой метод Эйлера представляет собой явный алгоритм. В неявных алгоритмах для определения xk+1 требуются дополнительные вычисления, но они по сравнению с аналогичными прямыми алгоритмами более устойчивы и дают более высокую точность вычислений (см. рис.3). Возможно это обусловлено наличием члена xk+1 в правой части формулы, что может рассматриваться как вид обратной связи.
Рис.3 Та же система, что и на рис.2 проинтегрирована от x0=1.0 с h=0.3 (a) прямым методом Эйлера, (b) обратным методом Эйлера. Точное решение показано штриховой линией.
Трапецеидальный алгоритм имеет вид:
Имеется целое семейство уравнений Рунге-Кутта второго порядка. Мы рассмотрим модифицированный алгоритм Эйлера-Коши, заданный соотношением:
Из этой формулы следует, что модифицированный алгоритм Эйлера-Коши включает два этапа. На первом этапе с помощью прямого метода Эйлера происходит перемещение на пол шага вперед к моменту времени (tk+h/2):
На втором этапе это промежуточное значение используется для аппроксимации векторного поля с помощью итераций Эйлера прямого типа:
Как и в случае алгоритма второго порядка метод Рунге-Кутта четвертого порядка относится к явным алгоритмам. Он использует промежуточные моменты времени для для вычисления состояния в момент времени tk+1. Следующие формулы определяют алгоритм Рунге-Кутта четвертого порядка:
Для определенных выше алгоритмов на каждом шаге требуется только одна начальная точка xk. Такие алгоритмы называются одношаговыми. Одношаговые алгоритмы высокого порядка имеют высокую точность, но они не эффективны, если велики затраты на вычисление f. Например, алгоритм Рунге-Кутта четвертого порядка требует на шаге четыре значения f. Кроме того, на текущем шаге не выполняются оценки функций с целью их использования на последующих шагах. В отличие от одношаговых алгоритмов, многошаговые алгоритмы повторно используют предыдущую информацию о траектории. В m-шаговом алгоритме для определения xk+1 используют m предыдущих точек xk, xk-1. xk-m+1 и значения f в этих точках. Общая формула m-шагового алгоритма имеет вид
Локальная ошибка определяется как ошибка на шаге алгоритма:
Для m-шаговых алгоритмов предполагается, что предыдущие m точек xk-i при i=0. m-1 являются точно заданными, т.е.
Ошибка усечения — это локальная ошибка, которая получилась бы в результате выполнения алгоритма на компьютере с бесконечной точностью. Другими словами, подразумевается, что эта локальная ошибка возникает помимо ошибки округления. Также важно помнить, что для m-шаговых алгоритмов предыдущие m точек xk. xk-m+1 предполагаются точно заданными. Ошибка усечения берет свое название от алгоритмов, основанных на рядах Тейлора (например, Рунге-Кутта). Эти алгоритмы были бы точными, если бы использовались полные (бесконечные) ряды. Ошибка возникает при усечении ряда до конечного числа членов. Ошибка усечения зависит только от алгоритма. Она не зависит от используемого компьютера и, следовательно, может быть проанализирована. Для алгоритмов Рунге-Кутта K-го порядка, при подходящих условиях, локальная ошибка усечения представляет собой
где «альфа» зависит от K, f и xk, но не зависит от h. Для многошагового алгоритма K-го порядка локальная ошибка усечения имеет вид
Глобальнае ошибка округления — это простое накопление локальных ошибок округления. Если локальная ошибка составляет «эпсилон», то ошибка округления на единичном интервале будет
Подобно ошибке округления, локальная ошибка усечения также накапливается с каждым шагом. Для одношаговых алгоритмов K-го порядка локальная ошибка усечения составляет
Если пренебречь зависимостью ak от xk, то на единичном интервале времени ошибка усечения будет:
