Алгоритм решения задачи коши для дифференциального уравнения

Данко П. Е., Попов А. Г., Кожевникова Т. Я. Высшая математика в упражнениях и задачах…
Часть II. Глава IV. Обыкновенные дифференциальные уравнения

§ 1. Дифференциальные уравнения первого порядка

1. Основные понятия. Дифференциальным уравнением называется уравнение, связывающее независимые переменные, их функцию и производные (или дифференциалы) этой функции. Если независимая переменная одна, то уравнение называется обыкновенным; если же независимых переменных две или больше, то уравнение называется дифференциальным уравнением в частных производных.

Наивысший порядок производной, входящей в уравнение, называется порядком дифференциального уравнения. Например:

1) х²у’ + 5xy = у² – обыкновенное дифференциальное уравнение первого порядка;

2) – обыкновенное дифференциальное уравнение второго порядка;

3) y’³ + y»y»’ = х – обыкновенное дифференциальное уравнение третьего порядка;

4) F (х, у, у’, у») = 0 – общий вид обыкновенного дифференциального уравнения второго порядка;

5) – уравнение в частных производных первого порядка.

В этом параграфе рассматриваются обыкновенные дифференциальные уравнения первого порядка, т. е. уравнения вида F (х, у, у’) = 0 или (в разрешенном относительно у’ виде) y’ = f(х, у).

Решением дифференциального уравнения называется такая дифференцируемая функция у = φ (x), которая при подстановке в уравнение вместо неизвестной функции обращает его в тождество. Процесс нахождения решения дифференциального уравнения называется интегрированием дифференциального уравнения.

Общим решением дифференциального уравнения первого порядка у’ = f(x, у) в области D называется функция у = φ(x, C), обладающая следующими свойствами: 1) она является решением данного уравнения при любых значениях произвольной постоянной С, принадлежащих некоторому множеству; 2) для любого начального условия у(х0) = у0 такого, что (x0; y0) ∈ 0, существует единственное значение С = С0, при котором решение у = φ(x, C0) удовлетворяет заданному начальному условию.

Всякое решение у = φ(x, C0), получающееся из общего решения у = φ (x, C) при конкретном значении С = С0, называется частным решением.

Задача, в которой требуется найти частное решение уравнения y’ = f(х, у) удовлетворяющее начальному условию у(х0) = y0, называется задачей Коши.

Построенный на плоскости хОу график всякого решения у = φ(х) дифференциального уравнения называется интегральной кривой этого уравнения. Таким образом, общему решению у = φ(х, С) на плоскости хОу соответствует семейство интегральных кривых, зависящее от одного параметра – произвольной постоянной С, а частному решению, удовлетворяющему начальному условию y(x0) = y0, – кривая этого семейства, проходящая через заданную точку М0(x0; у0).

Если функция f(х, у) непрерывна и имеет непрерывную производную в области D, то решение дифференциального уравнения у’= f (х, у) при начальном условии у(х0) = у0 существует и единственно, т. е. через точку (x0; y0) проходит единственная интегральная кривая данного уравнения (теорема Коши).

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

Особые решения не получаются из общего решения дифференциального управления ни при каких значениях произвольной постоянной С (в том числе и при С = ± ∞).

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

Например, общее решение уравнения записывается в виде у = sin (х + С). Это семейство интегральных кривых имеет две огибающие: у = 1 и у = -1, которые и будут особыми решениями.

2. Дифференциальные уравнения с разделяющимися переменными. Дифференциальное уравнение вида

относится к типу уравнений с разделяющимися переменными. Если ни одна из функций f1(x), f2(y), φ1(x), φ2(y) не равна тождественно нулю, то в результате деления исходного уравнения на f2 (x) φ1 (y) оно приводится к виду

Почленное интегрирование последнего уравнения приводит к соотношению

которое и определяет (в неявной форме) решение исходного уравнения. (Решение дифференциального уравнения, выраженное в неявной форме, называют интегралом этого уравнения.)

507. Решить уравнение х(у²-4)dx + y dy = 0.

△ Разделив обе части уравнения на у² – 4 ≠ 0, имеем

x² + ln|у² – 4| = ln|C|, или у² – 4 = Сe -λ²

Это общее решение данного дифференциального уравнения.

Пусть теперь у² – 4 = 0, т. е. у = ± 2. Непосредственной подстановкой убеждаемся, что у = ±2 – решение исходного уравнения. Но оно не будет особым решением, так как его можно получить из общего решения при С = 0. ▲

508. Найти частный интеграл уравнения у’ cos х = у / ln у, удовлетворяющий начальному условию y(0) = l.

△ Полагая , перепишем данное уравнение в виде

Проинтегрируем обе части уравнения:

, или

Используя начальное условие у = 1 при х = 0, находим С = 0. Окончательно получаем

▲

509. Найти общий интеграл уравнения у’ = tg x tg y.

△ Полагая и разделяя переменные, приходим к уравнению ctg у dy = tg х dx. Интегрируя, имеем

, или ln|sin у| = -ln|cos x| + ln С.

Отсюда находим sin y = C/cos x, или sin y / cos x = С (общий интеграл). ▲

510. Найти частное решение дифференциального уравнения (l + x²)dy + y dx = 0 при начальном условии у(1) = 1.

△ Преобразуем данное уравнение к виду . Интегрируя, получим

, или ln |y| = – arctg x + С

Это и есть общий интеграл данного уравнения.

Теперь, используя начальное условие, найдем произвольную постоянную С; имеем ln 1 = — arctg 1 + С, т. е. С = π/4. Следовательно,

ln у = – arctg х + π/4,

откуда получаем искомое частное решение y = e π/4 – arctg x . ▲

Данко П. Е., Попов А. Г., Кожевникова Т. Я. Высшая математика в упражнениях и задачах… Ч. II. Стр. 117-119.

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

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

= \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)| —>
Порядок дифференциального уравнения и его решения, задача Коши
Обыкновенным дифференциальным уравнением называется уравнение, связывающее независимую переменную, неизвестную функцию этой переменной и её производные (или дифференциалы) различных порядков.
Порядком дифференциального уравнения называется порядок старшей производной, содержащейся в нём.
Кроме обыкновенных изучаются также дифференциальные уравнения с частными производными. Это уравнения, связывающие независимые переменные , неизвестную функцию этих переменных и её частные производные по тем же переменным. Но мы будем рассматривать только обыкновенные дифференциальные уравнения и поэтому будем для краткости опускать слово «обыкновенные».
Примеры дифференциальных уравнений:
(1) ;
(2) ;
(3) ;
(4) ;
(5) .
Уравнение (1) — четвёртого порядка, уравнение (2) — третьего порядка, уравнения (3) и (4) — второго порядка, уравнение (5) — первого порядка.
Дифференциальное уравнение n-го порядка не обязательно должно содержать явно функцию, все её производные от первого до n-го порядка и независимую переменную. В нём могут не содержаться явно производные некоторых порядков, функция, независимая переменная.
Например, в уравнении (1) явно нет производных третьего и второго порядков, а также функции; в уравнении (2) — производной второго порядка и функции; в уравнении (4) — независимой переменной; в уравнении (5) — функции. Только в уравнении (3) содержатся явно все производные, функция и независимая переменная.
Решением дифференциального уравнения называется всякая функция y = f(x), при подстановке которой в уравнение оно обращается в тождество.
Процесс нахождения решения дифференциального уравнения называется его интегрированием.
Пример 1. Найти решение дифференциального уравнения .
Решение. Запишем данное уравнение в виде . Решение состоит в нахождении функции по её производной. Изначальная функция, как известно из интегрального исчисления, есть первообразная для , т. е.
.
Это и есть решение данного дифференциального уравнения. Меняя в нём C, будем получать различные решения. Мы выяснили, что существует бесконечное множество решений дифференциального уравнения первого порядка.
Общим решением дифференциального уравнения n-го порядка называется его решение, выраженное явно относительно неизвестной функции и содержащее n независимых произвольных постоянных, т. е.
Решение дифференциального уравнения в примере 1 является общим.
Частным решением дифференциального уравнения называется такое его решение, в котором произвольным постоянным придаются конкретные числовые значения.
Пример 2. Найти общее решение дифференциального уравнения и частное решение при .
Решение. Проинтегрируем обе части уравнения такое число раз, которому равен порядок дифференциального уравнения.
,
,
.
В результате мы получили общее решение —
данного дифференциального уравнения третьего порядка.
Теперь найдём частное решение при указанных условиях. Для этого подставим вместо произвольных коэффициентов их значения и получим
.
Если кроме дифференциального уравнения задано начальное условие в виде , то такая задача называется задачей Коши. В общее решение уравнения подставляют значения и и находят значение произвольной постоянной C, а затем частное решение уравнения при найденном значении C. Это и есть решение задачи Коши.
Пример 3. Решить задачу Коши для дифференциального уравнения из примера 1 при условии .
Решение. Подставим в общее решение значения из начального условия y = 3, x = 1. Получаем
.
Записываем решение задачи Коши для данного дифференциального уравнения первого порядка:
.
При решении дифференциальных уравнений, даже самых простых, требуются хорошие навыки интегрирования и взятия производных, в том числе сложных функций. Это видно на следующем примере.
Пример 4. Найти общее решение дифференциального уравнения .
Решение. Уравнение записано в такой форме, что можно сразу же интегрировать обе его части.
.
Применяем метод интегрирования заменой переменной (подстановкой). Пусть , тогда .
Требуется взять dx и теперь — внимание — делаем это по правилам дифференцирования сложной функции, так как x и есть сложная функция («яблоко» — извлечение квадратного корня или, что то же самое — возведение в степень «одна вторая», а «фарш» — самое выражение под корнем):
Возвращаясь к переменной x, получаем:
.
Это и есть общее решение данного дифференциального уравнения первой степени.
Не только навыки из предыдущих разделов высшей математики потребуются в решении дифференциальных уравнений, но и навыки из элементарной, то есть школьной математики. Как уже говорилось, в дифференциальном уравнении любого порядка может и не быть независимой переменной, то есть, переменной x. Помогут решить эту проблему не забытые (впрочем, у кого как) со школьной скамьи знания о пропорции. Таков следующий пример.
Пример 5. Найти общее решение дифференциального уравнения .
Решение. Как видим, переменная x в уравнении отсутствует. Вспоминаем из курса дифференциального исчисления, что производная может быть записана также в виде . В результате уравнение приобретает вид
,
то есть, в нём в некотором виде появился x.
Теперь вспомнаем одно из свойств пропорции: из пропорции выткают следующие пропорции:
,
то есть в пропорции можно менять местами крайние и средние члены или те и другие одновременно.
Применяя это свойство, преобразуем уравнение к виду
,
после чего интегрируем обе части уравнения:
.
Оба интеграла — табличные, находим их:
и получаем решение данного дифференциалного уравнения первого порядка:
.
Эта статья представила необходимый минимум сведений о дифференциальных уравнениях и их решениях и должна помочь вам уверенно и увлечённо перейти к изучению различных видов дифференциальных уравнений.

источники:
http://slemeshevsky.github.io/num-mmf/ode/html/._ode-FlatUI001.html
http://function-x.ru/differential_equations1.html