Численное решение обыкновенных дифференциальных уравнений

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

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

Постановка задачи:

Требуется найти функцию у = у(х), удовлетворяющую уравнению

(3)

и принимающую при х = х₀ заданное значение у₀.

(4)

При этом решение необходимо получить в интервале х₀ ≤ х ≤ х_к. Из теории дифференциальных уравнений известно, что решение у(х) задачи Коши (3), (4) существует, единственно и является гладкой функцией, если правая часть F(х, у) удовлетворяет некоторым условиям гладкости. Численное решение задачи Коши методом Рунге-Кутта 4-го порядка заключается в следующем. На заданном интервале [ х₀, х_к ] выбираются узловые точки. Значение решения в ну­левой точке известно у(х ₀) = у₀. В следующей точке у(х₁) определяется по фор­муле

(5)

здесь

k₀ = hF(х₀,у₀); k₁=hF(х₀+h/2, у₀+k₀/ 2);

k₂ = hF(х₀ + h /2, у₀+k₁/ 2); к₃ = hF(х₀ +h, у₀ + k ₂); h - шаг сетки, (6)

т. е. данный вариант метода Рунге-Кутта требует на каждом шаге четырехкратного вычисления правой части уравнения (3). Этот алгоритм реализован в программе ode45. Кроме этой программы МАТLАВ располагает обширным набором аналогичных программ, позволяющих успешно решать обыкновенные дифференциальные уравнения.

Пример 5. Решить задачу Коши

на [0,1] при y(0)=1 (7)

Точное решение имеет вид

Выполним решение данной задачи с помощью программы ode45. Вначале в М-файл записываем правую часть уравнения (7), сам М-файл оформляется как файл - функция, даем ему имя F:

Function dydx = F(x,y)

dydx = zeros (1,1);

dydx( 1) = 2*(х^{^}2+у(1));

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

Протокол программы.

>>[X Y]=ode45 (@ F, [0 1], [1]);

% Дескриптор @ обеспечивает связь с файлом - функцией правой части

% [0 1 ] - интервал на котором необходимо получить решение

% [1] - начальное значение решения

>>р1оt (Х,У);

>> % Построение графика численного решения задачи Коши (7)

>> hold on; gtext (' у(х) ')

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

>> [X Y]

>>% Последняя команда выводит таблицу численного решения задачи.

Результаты решения. График решения задачи Коши (7) показан на рис. 3. Численное решение представлено в таблице 4, где приведены только отдельные узловые точки. В программе ode45по умолчанию интервал разбивается на 40 точек с шагом h = 1/40 = 0.025.

Рис. 3 – График решения задачи Коши

Таблица 4

Как следует из таблицы 4 численное решение программой ode45явля­ется точным.

Варианты заданий. Построить график и вывести в виде таблицы решение задачи Коши на интервале [0; 1] методом Рунге-Кутта 4-го порядка. Данные взять из таблицы 5.

Таблица 5

№ п/п	f (х,у)	y0
1.		0.0
2.		0.1
3.	е+х+3у	2.0
4.		0.3
5.		0.4
6.	1/(1 + у2)+ x2	0.0
7.	1/(1 + у2)+ху	0.1
8.	cos у + xу	0.2
9.	х2 соs у + 0.1	0.3
10.	x 3 cos у + 0.1	0.4
11.	cos(ху)- 0.5	0.5
12.	e-у +ех-2	0.0
13.		0.5
14.	e-xy+1	0.4
15.		0.3
16.		0.2
17.		0.1
18.		0.0
19.		0.1
20.		0.2
21.		0.3
22.		0.4
23.		0.5
24.		0.6
25.		0.7
26.		0.0
27.		0.1
28.		0.2
29.		0.3
30.		0.4