Решение задачи коши методом рунге кутта

Формулировка задачи Коши для ДУ первого порядка:

Дано ОДУ первого порядка, разрешенное относительно производной

Необходимо найти решение (1), удовлетворяющее начальному условию

То есть в задаче Коши необходимо найти кривую , проходящую через заданную точку .

Решение задачи Коши является частным решением уравнения (1) при условии (2).

Достаточные условия существования и единственности задачи Коши содержатся в следующей теореме.

Теорема. Пусть функция — правая часть ДУ (1) непрерывна вместе со своей частной производной в некоторой области D на плоскости . Тогда при любых начальных значениях задача Коши (1) (2) имеет единственное решение .

При выполнении условий теоремы через точку проходит единственная кривая.

В дальнейшем будем считать, что условия теоремы выполнены.

2. Метод Рунге

Это простейший численный метод. Рассмотрим задачу Коши

Зададим равномерную сетку

Введём обозначения . Тогда имеем вычислительную формулу для метода Рунге Кутта первого порядка

Данная формула позволяет, начиная от начального условия , найти последовательно величины с шагом h и, таким образом, решить задачу Коши.

Пример. Построить вычислительную формулу для решение ДУ

методом Эйлера с шагом на отрезке .

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

Здесь — точное решение; — численное решение;

Метод Эйлера есть метод первого порядка точности.

3. Метод Рунге

На равномерной сетке имеем формулу Рунге Кутта второго порядка точности:

4. Метод Рунге

Это наиболее распространенный метод решения задачи Коши. Вычислим интеграл в (3) по формуле Симпсона. Получим вычислительную формулу:

Схемы Рунге Кутта имеют ряд достоинств:

1) все они (кроме метода Эйлера) имеют хорошую точность;

2) они являются явными, т. е. значения вычисляются по ранее найденным значениям ;

3) схемы допускают введение переменного шага .

5. Правило Рунге оценки погрешности

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

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

Изменение шага для методов Рунге — Кутты сложности не представляет. Оценить погрешность достаточно сложно, так как простые способы оценки погрешности отсутствуют.

Поскольку предельное значение погрешности для методов Рунге Кутта пропорционально величине ( k порядок точности метода) или

то здесь, как и в случае численного интегрирования методами Ньютона Котеса, возможно использование правила Рунге.

Пусть , — решение методом Рунге Кутта в точке ; с шагом 2 h ; , — решение в точке с шагом h , — точное решение в точке с. Тогда погрешность

Вывод формулы аналогичен выводу формулы для правила Рунге в методах численного интегрирования.

Замечание. Данная оценка достаточно грубая и тем точнее, чем выше порядок точности метода.

6. Решение систем ОДУ первого порядка методом Рунге — Кутта

Формулировка задачи Коши на сетке , , :

Решить систему ОДУ первого порядка:

где , — неизвестные функции, удовлетворяющие начальным условиям

Задача Коши в векторной форме:

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

Условия существования и единственности решения данной задачи Коши такие же , как и для ОДУ первого порядка:

Функции должны быть непрерывными вместе со своими частными производными , , в некоторой области D в пространстве .

При выполнении этих условий через точку проходит единственная совокупность кривых:

Величина погрешности при решении систем ОДУ оценивается значением

где — погрешность решения на сетке с шагом h в точке :

Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений $$ egin ag <1>frac &= f_i (t, u_1, u_2, ldots, u_n), quad t > 0\ ag <2>u_i(0) &= u_i^0, quad i = 1, 2, ldots, m. end $$

Используя векторные обозначения, задачу (1), (2) можно записать как задачу Коши $$ egin ag <3>frac> &= pmb(t, pmb), quad t > 0, \ ag <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 ) при ( au o 0 ).

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

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

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

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

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

Четвертый шаг заключается в получении численного алгоритма. Из (5) следует, что мы должны знать значение ( y^n ) для того, чтобы решить уравнение (5) относительно ( y^ ) и получить формулу для нахождения приближенного значения искомой функции на следующем временном слое ( t_ ): $$ egin ag <6>pmb^ = pmb^n + au 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) используется следующий метод: $$ egin ag <7>frac<pmb^ — pmb^n> < au>= pmb(t_, pmb^). end $$

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

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

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

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

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

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

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

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

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

В методах Рунге—Кутта в вычислениях участвуют значения приближенного решения только в двух соседних узлах ( pmb^n ) и ( pmb^ ) — один шаг по переменной ( t ). Линейный ( m )-шаговый разностный метод записывается в виде $$ egin ag <13>frac<1> < au>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
e 0 ). Для начала расчетов по рекуррентной формуле (13) необходимо задать ( m ) начальных значений ( pmb^0 ), ( pmb^1 ), ( dots ), ( pmb^ ) (например, можно использовать для их вычисления метод Эйлера).

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

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

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

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

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

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

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

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

= sum_^ frac<partial f_i> <partial u_j>(t, w) v + ar(t), quad t > 0. end $$

Пусть ( lambda_i(t) ), ( i = 1, 2, ldots, m ) — собственные числа матрицы $$ egin A(t) = < a_(t) >, quad a_(t) = frac<partial f_i><partial u_j>(t, w). end $$ Система уравнений (3) является жесткой, если число $$ egin 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) область его устойчивости содержит угол $$ egin |arg(-mu)| —>

В работе [2] рассмотрена постановка и решение задачи Коши методами Эйлера и Рунге–Кутта. В данной работе мы ограничимся программами расчета в рамках пакета Mathcad для решения задачи Коши дифференциального уравнения первого порядка, разрешенного относительно производной.

y ′ = f ( x , y ), x [ a , b ], y ( a ) = y 0.

В качестве тестовой задачи выберем несложную задачу Коши, имеющую аналитическое решение.

y ′ = 0,5 e x y 2 , x [0,1], y (0) = 1.

Это легко интегрирующее уравнение с разделяющимися переменными. Общее решение получается следующим образом: