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

В этой лабораторной работе мы рассмотрим некоторые приближенные методы решения задачи Коши, состоящей в отыскании решения y(x) дифференциального уравнения y' = f(x,y), удовлетворяющего заданному начальному условию y(x₀) = y₀.

Задачу приближенного решения задачи Коши будем понимать как задачу построения на заданном отрезке [x₀, b] функции φ(x), которая «близка» к решению y(x) задачи Коши с заданной точностью ε в том смысле, что |y(x) – φ(x)| ≤ ε (x₀ ≤ x≤ b).

Метод Эйлера.

Рассмотрим задачу Коши для уравнения 1-го порядка:

Будем искать численное решение уравнения на отрезке [x₀, b]. Зададим на этом отрезке сетку {x_i, i = 0,1,…,N}, таким образом, чтобы x₀ < x₁ < ... < x_N = b. Введем обозначение для шага сетки h_i = x_i+1 – x_i, i = 0,…, N–1. Заменив производную в уравнении правой разностью, получим

Известно, что y₀ = y(x₀) = η. Откуда можно найти все остальные значения y_i по рекуррентной формуле:

Данный метод нахождения численного решения называется методом Эйлера (или методом ломаных). Схемы, в которых значение функции явно выражается через уже найденные значения, называются явными, иначе – неявными. Таким образом, схема Эйлера является явной. Оценка погрешности для данного метода дает O(max(h_i)), что предполагает малый шаг сетки для получения удовлетворительного решения.

На каждом из отрезков [x_i, x_i+1] полученное решение будет представлять собой отрезок прямой, проведенной через точку (x_i, y_i) с угловым коэффициентом f(x_i, y_i). Такая геометрическая интерпретация решения объясняет второе название метода (метод ломаных).

Пример П. 1. Найти численное решение следующей задачи Коши на отрезке x ϵ [0,1] методом Эйлера (на равномерной сетке с шагом h = 0.1) и сравнить его с аналитическим:

Реализуем метод Эйлера в виде файл-функции:

function [yy,xx]=euler(f,x0,y0,xe,h)

xx = x0:h:xe; % значения координат х для расчета

% выделяем память для значений функции:

yy = zeros(length(y0),length(xx));

yy(:,1) = y0; % начальное значение y

% последовательное вычисление значений в цикле

for i=1:length(xx)-1

yy(:,i+1) = yy(:,i) + h*f(xx(i),yy(:,i));

end

end

Нам понадобится также вспомогательный файл для функции в правой части уравнения:

function f = f(x,y)

f=x.^2;

end

С помощью созданных функций найдем решение задачи Коши:

>> [y,x]=euler(@f,0,1,1,0.1)

y =

1.0000 1.0000 1.0010 1.0050 1.0140 1.0300 1.0550 1.0910 1.1400 1.2040 1.2850

x =

0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000

Отобразим на графике полученное приближенное решение и точное решение, найденное аналитически (y = x³/3 + 1). На рис.П. 1 сплошная кривая – аналитическое решение, пунктирная ломаная со звездочками – приближенное решение.

Рис.П. 1. Точное и приближенное решения задачи Коши

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

1) Найти решение задачи на отрезке xϵ[0,1] методом Эйлера, используя равномерную сетку с шагом h = 0.1.

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

3) Построить графики точного и численного решений (на одном рисунке) и график ошибки вычислений (на отдельном рисунке).

4) Повторить решение задачи для шага сетки h = 0.05. К графикам, построенным в пункте 3, добавить график нового численного решения и график соответствующей ему ошибки.

2. Решение уравнений p-го порядка и систем обыкновенных дифференциальных уравнений.

Рассмотрим обыкновенное дифференциальное уравнение (ОДУ) p-го порядка:

Путем введения замены {y^(k) = y_k, k = 1,…, p – 1} данное уравнение можно свести к системе дифференциальных уравнений 1-го порядка:

Для получения единственного решения системы нужно наложить p дополнительных условий на функции y_k(x). Для задачи Коши данные условия задаются в одной точке:

Для нахождения решения полученной задачи Коши может быть применен метод Эйлера, рассмотренный ранее. Представим полученную задачу в векторном виде

введя обозначения:

В данных обозначениях формулы метода Эйлера могут быть представлены в виде

Обратите внимание, что функция для метода Эйлера, описанная в примере П. 1, изначально была рассчитана на работу с переменными y, y₀, заданными в виде вектора-столбца. Вносить дополнительные изменения в эту функцию для решения систем уравнений не требуется. Достаточно лишь правильно задать параметры и реализовать вычисление правой части.

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

1) Найти решение задачи на отрезке x ϵ [0,1] методом Эйлера, используя равномерную сетку с шагом h = 0.1 и предварительно преобразовав задачу в задачу Коши для системы уравнений 1-го порядка.

2) Найти ошибку вычислений как разность между точным решением задачи (можно найти вручную или в символьном виде с помощью функции dsolve) и полученным численным решением в каждой точке сетки. Определить погрешность решения как максимум модуля ошибки вычислений.

3) Построить графики точного и численного решений (на одном рисунке) и график ошибки вычислений (на отдельном рисунке).

4) Повторить решение задачи для шага сетки h = 0.05. К графикам, построенным в пункте 3, добавить график нового численного решения и график соответствующей ему ошибки.

Методы Рунге – Кутты.

Для погрешности метода Эйлера справедлива оценка ε = O(h), т.е. метод Эйлера имеет первый порядок точности. Это означает, что для уменьшения погрешности вычислений в 100 раз шаг сетки также необходимо уменьшить в 100 раз. Для большинства практических задач, к сожалению, такой низкий порядок точности не достаточен.

Для решения поставленной задачи с более высоким порядком точности было придумано целое семейство методов, получивших название методы Рунге – Кутты. Наиболее популярный из этих методов – метод Рунге – Кутты 4-го порядка точности, который мы рассмотрим далее. На практике также часто применяют метод 2-го порядка и уже знакомый нам метод Эйлера (он же метод Рунге – Кутты 1-го порядка). В литературе можно найти методы Рунге – Кутты до 8-го порядка точности включительно, но практического распространения они не получили.

Для поиска решения задачи Коши с помощью методов Рунге – Кутты в области решения вводится равномерная сетка, и значения функции вычисляются последовательно в узлах сетки, начиная с известного значения в точке x₀. В общем виде формула для вычисления нового значения функции задается как:

Каждый конкретный метод семейства Рунге – Кутты определяется числом промежуточных этапов (стадий) s и фиксированными значениями коэффициентов a_j,l, b_j, c_j. Значения коэффициентов подбираются таким образом, чтобы при заданном порядке точности число требуемых стадий было минимальным.

В частности, метод Эйлера имеет одну стадию и коэффициенты b₁ = 1, c₁ = 0. Для достижения 2-го порядка точности необходимо две стадии (используемые для этого значения коэффициентов можно легко найти в литературе). Для достижения 4-го порядка точности требуется использовать метод Рунге – Кутты с четырьмя стадиями:

Метод с таким набором коэффициентов получил название метода Рунге – Кутты 4-го порядка. Для его погрешности справедлива оценка ε = O(h⁴), т.е. при уменьшении шага сетки в 10 раз, погрешность уменьшается в 10000 раз.

Данный метод можно применять как для одного уравнения, так и для системы уравнений, используя переход в векторную форму по правилам, описанным ниже в примере П. 2.

Пример П.2. Написать файл-функцию, реализующую метод Рунге – Кутты 4-го порядка. Предусмотреть возможность использования метода для системы уравнений.

function [yy,xx]=runge(f,x0,y0,xe,h)

xx = x0:h:xe;

% yy создадим в виде матрицы, i-я строка которой соответсвует i-й функции

yy = zeros(length(y0),length(xx));

yy(:,1) = y0;

for i=1:length(xx)-1

k1 = h * f(xx(i), yy(:,i));

k2 = h * f(xx(i)+h/2, yy(:,i)+k1/2);

k3 = h * f(xx(i)+h/2, yy(:,i)+k2/2);

k4 = h * f(xx(i)+h, yy(:,i)+k3);

yy(:,i+1)=yy(:,i)+(k1+2*k2+2*k3+k4)/6;

end

end

Упражнение П.3.Выполнить задания упражнений 1 и 2, применяя вместо метода Эйлера метод Рунге – Кутты 4-го порядка. Сравнить погрешности, полученные при использовании разных методов.