Метод Эйлера

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

Для получения единственного решения системы нужно наложить 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, добавить график нового численного решения и график соответствующей ему ошибки.

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-го порядка. Сравнить погрешности, полученные при использовании разных методов.