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

Метод Рунге-Кутта применим и к задаче Коши для системы m дифференциальных уравнений с m неизвестными функциями

(7.18)

(7.19)

Приведем для задачи (7.18), (7.19) расчетные формулы метода Рунге-Кутта четвёртого порядка. Пусть требуется найти систему m функций удовлетворяющих в интервале (x₀, X) дифференциальным уравнениям (7.18), а в точке x_{0 —} начальным условиям (7.19). Предположим, что отрезок [x₀, X] разбит на N частей

Тогда каждую l-ую функцию y_l(x) можно приближенно вычислять в точках x_i+1 по формулам Рунге-Кутта

(7.20)

Здесь через y_l,i мы обозначили приближенное значение функции y_l(x) в точке x_i.

Обращаем внимание на порядок вычислений по формулам (7.20). На каждом шаге сначала вычисляются коэффициенты K_l,i в следующем порядке

K_1,1, K_2,1, …, K_m,1,

K_1,2, K_2,2, …, K_m,2,

K_1,3, K_2,3, …, K_m,3,

K_1,4, K_2,4, …, K_m,4,

и лишь затем приближенные значения функций y_1,i+1, y_2,i+1,…, y_m,i+1.

Пример 7.4. Решить на отрезке [1, 2] задачу Коши

Решение в программе Excel.Точным решением задачи является система функций .

Разделим отрезок [1, 2] на N = 10 частей. Тогда шаг равен h = 0,1.

Число неизвестных функций m = 3. Число коэффициентов K_l,i равно 12. Создадим функции для вычисления правых частей f_l(x, y₁, y₂, y₃) системы дифференциальных уравнений и коэффициентов K_l,i:

Function f_1(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3)

f_1 = 2 * x * y_1 / y_2

End Function

Function f_2(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3)

f_2 = 8 * y_3 / (3 * y_2)

End Function

Function f_3(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3)

f_3 = 3 * y_3 / y_1

End Function

Function K_1_1(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3)

K_1_1 = f_1(x, y_1, y_2, y_3)

End Function

Function K_2_1(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3)

K_2_1 = f_2(x, y_1, y_2, y_3)

End Function

Function K_3_1(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3)

K_3_1 = f_3(x, y_1, y_2, y_3)

End Function

Function K_1_2(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3, ByVal h)

K11 = K_1_1(x, y_1, y_2, y_3)

K21 = K_2_1(x, y_1, y_2, y_3)

K31 = K_3_1(x, y_1, y_2, y_3)

K_1_2 = f_1(x + h / 2, y_1 + h * K11 / 2, y_2 + h * K21 / 2, y_3 + h * K31 / 2)

End Function

Function K_2_2(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3, ByVal h)

K11 = K_1_1(x, y_1, y_2, y_3)

K21 = K_2_1(x, y_1, y_2, y_3)

K31 = K_3_1(x, y_1, y_2, y_3)

K_2_2 = f_2(x + h / 2, y_1 + h * K11 / 2, y_2 + h * K21 / 2, y_3 + h * K31 / 2)

End Function

Function K_3_2(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3, ByVal h)

K11 = K_1_1(x, y_1, y_2, y_3)

K21 = K_2_1(x, y_1, y_2, y_3)

K31 = K_3_1(x, y_1, y_2, y_3)

K_3_2 = f_3(x + h / 2, y_1 + h * K11 / 2, y_2 + h * K21 / 2, y_3 + h * K31 / 2)

End Function

Function K_1_3(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3, ByVal h)

K12 = K_1_2(x, y_1, y_2, y_3, h)

K22 = K_2_2(x, y_1, y_2, y_3, h)

K32 = K_3_2(x, y_1, y_2, y_3, h)

K_1_3 = f_1(x + h / 2, y_1 + h * K12 / 2, y_2 + h * K22 / 2, y_3 + h * K32 / 2)

End Function

Function K_2_3(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3, ByVal h)

K12 = K_1_2(x, y_1, y_2, y_3, h)

K22 = K_2_2(x, y_1, y_2, y_3, h)

K32 = K_3_2(x, y_1, y_2, y_3, h)

K_2_3 = f_2(x + h / 2, y_1 + h * K12 / 2, y_2 + h * K22 / 2, y_3 + h * K32 / 2)

End Function

Function K_3_3(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3, ByVal h)

K12 = K_1_2(x, y_1, y_2, y_3, h)

K22 = K_2_2(x, y_1, y_2, y_3, h)

K32 = K_3_2(x, y_1, y_2, y_3, h)

K_3_3 = f_3(x + h / 2, y_1 + h * K12 / 2, y_2 + h * K22 / 2, y_3 + h * K32 / 2)

End Function

Function K_1_4(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3, ByVal h)

K13 = K_1_3(x, y_1, y_2, y_3, h)

K23 = K_2_3(x, y_1, y_2, y_3, h)

K33 = K_3_3(x, y_1, y_2, y_3, h)

K_1_4 = f_1(x + h, y_1 + h * K13, y_2 + h * K23, y_3 + h * K33)

End Function

Function K_2_4(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3, ByVal h)

K13 = K_1_3(x, y_1, y_2, y_3, h)

K23 = K_2_3(x, y_1, y_2, y_3, h)

K33 = K_3_3(x, y_1, y_2, y_3, h)

K_2_4 = f_2(x + h, y_1 + h * K13, y_2 + h * K23, y_3 + h * K33)

End Function

Function K_3_4(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3, ByVal h)

K13 = K_1_3(x, y_1, y_2, y_3, h)

K23 = K_2_3(x, y_1, y_2, y_3, h)

K33 = K_3_3(x, y_1, y_2, y_3, h)

K_3_4 = f_3(x + h, y_1 + h * K13, y_2 + h * K23, y_3 + h * K33)

End Function

Введем эти функции в программу Excel. Для этого с помощью меню «Сервис — Макросы — Редактор Visual Basic» откроем окно редактора и выполним команду «Insert — Module» и далее вводим описания функций.

После этого переходим в окно программы Excel и введем данные и формулы, как показано в таблице 7.9. Строки 1 и 2 и столбец A заполняем как в таблице. В ячейки B3, C3 и D3 вводим начальные значения функций.

Ячейке B1 присвоим имя h. В следующие ячейки введем указанные формулы и маркером заполнения копируем их вниз до строки 13 (табл. 7.8):

Таблица 7.8

В результате вычислений в столбцах B, C и D получаем приближенные значения функций y₁, y₂, y₃.

Таблица 7.9

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

7.1.6. Задачи Коши для дифференциальных уравнений
высших порядков

Задача Коши для дифференциального уравнения n-го порядка

, (7.21)

(7.22)

легко сводится к задаче Коши для системы дифференциальных уравнений с помощью замены переменных

(7.23)

Учитывая (7.23) из уравнения (7.21) получим систему дифференциальных уравнений

(7.24)

Начальные условия (7.22) для функций z_l переписываются в виде

. (7.25)

Запишем для полученной системы метод Рунге-Кутта

. (7.26)

Для вычисления коэффициентов K_l,1, K_l,2, K_l,3 и K_l,4 имеем следующие формулы:

Пример 7.5. Решить задачу Коши для дифференциального уравнения третьего порядка

,

Решение в Excel. Точное решение данной задачи известно . Преобразуем эту задачу с помощью замены к задаче для системы уравнений:

Точным решением системы теперь будут функции

.

Выберем шаг h = 0,1. В программе Excel мы используем функции и алгоритм решения, созданные для решения примера 7.4.

Откроем файл примера 7.4 и внесем необходимые изменения. Результаты расчетов представлены в таблице 7.10.

1) Заменить функции f_1(), f_2() и f_3() на следующие:

Function f_1(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3)

f_1 = 2 * x * y_1 / y_2

End Function

Function f_2(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3)

f_2 = 8 * y_3 / (3 * y_2)

End Function

Function f_3(ByVal x, ByVal y_1, ByVal y_2, ByVal y_3)

f_3 = 3 * y_3 / y_1

End Function

Функции для вычисления коэффициентов K_l,I не изменяем.

2) Изменим начальные условия в ячейках B3, C3, D3 на 0, 1, 0.

В столбце значений переменной x_i вводим числа 0; 0,1; 0,2; …; 1.

В ячейках E3, F3, G3 вводим соответственно формулы =SIN(A3), =COS(A3) и =-SIN(A3).

Таблица 7.10

Сравнивая столбцы B и E, видим, что метод Рунге-Кутта дает хорошие приближения к значениям точного решения y1=sinx.