Оценка погрешностей методов. Как уже отмечалось выше, при использовании метода Эйлера на каждом шаге положение

Как уже отмечалось выше, при использовании метода Эйлера на каждом шаге положение следующей точки решения, вычисляется с использованием информации о правой части f (x, y) уравнения только в крайней левой точке отрезка [ x_k, x_{k +1}]. В результате его абсолютная погрешность на каждом шаге пропорциональна величине h ²

,

где

.

В формуле для C_{k +1} в качестве функции используется некая промежуточная функция, график которой располагается между графиками приближённого и неизвестного точного решений.

Аналогичные соотношения для оценки абсолютных погрешностей на каждом шаге можно привести для усовершенствованного метода Эйлера

,

где

и метода Рунге–Кутта 4-го порядка

,

где формула для C_{k +1} не приводится из-за своей громоздкости.

Так как на отрезке интегрирования дифференциального уравнения требуется выполнить несколько шагов, то общая абсолютная погрешность вычисления приближённого решения в точке x_n будет складываться из погрешностей каждого отдельного шага и погрешностей, которые они порождают на последующих шагах. Она имеет оценку

,

где m = 2 для метода Эйлера, m = 3 для усовершенствованного метода Эйлера и m = 5 для метода Рунге–Кутта 4-го порядка, а параметры C _max и q определяются следующим образом

, .

Такое вычисление абсолютных погрешностей решения задачи Коши затруднено, так как требует решения дополнительной, зачастую даже более сложной, задачи поиска максимума высших производных правой части уравнения. Поэтому чаще для вычисления погрешности методов используют апостериорные оценки, базирующиеся на правиле Рунге (правило двойного счёта), идея которого описана в разделе «Вычисление определённых интегралов». Его основное соотношение, справедливое для всех способов приближённого решения задачи Коши, имеет вид

,

где y (x_k, h) и y (x_k, 2 h) – приближённые значения решения, вычисленные в точке x_k при шагах интегрирования, отличающихся друг от друга в два раза.

В качестве относительной погрешности решения задачи Коши на отрезке [ x ₀, x_n ] его построения, часто используют её интервальную оценку, равную отношению абсолютной погрешности к максимальному значению модуля решения на данном отрезке

.

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

с начальным условием y (0) = 0 на отрезке [0, 0.4] с шагом h = 0.1.

1. Решение задачи методом Эйлера. В соответствии с алгоритмом метода рас­чётная схема может быть представлена в виде следующих соотношений

.

Поэтому, пошаговый процесс построения решения с заданным шагом интегрирования h = 0.1 будет выглядеть следующим образом

,

,

,

,

.

Для оценки погрешности полученного решения необходимо повторить про­деланные расчёты с удвоенным шагом интегрирования h = 0.2

,

,

.

Полученные результаты позволяют оценить абсолютную и относительную погрешности решения с шагом интегрирования h = 0.1

,

,

,

.

2. Решение задачи усовершенствованным методом Эйлера. В данном случае рас­чётная схема метода представляется в виде соотношений

В соответствии с этим процесс построения решения с шагом интегрирова­ния h = 0.1 будет иметь вид

,

,

Аналогично предыдущему примеру, где реализован метод Эйлера, для вы­числения погрешности решения с шагом h = 0.1 повторяются расчёты с шагом h = 0.2

,

,

Отсюда абсолютная и относительная погрешности решения с шагом h = 0.1будут соответственно равны

,

,

,

,

.

Программное обеспечение

Система Matlab позволяет решать задачу Коши для обыкновенных дифференциальных уравнений и их систем посредством обращения к функциям ode23 и ode45. Первая из них реализует методы Рунге–Кутта второго и третьего порядков, а вторая – четвертого и пятого порядков. Выбор шага интегрирования в них осуществляется автоматически на основе задаваемой абсолютной погрешности решения. Решение задачи эти функции возвращают в виде массива значений аргумента и искомой функции.

Обращение и функциям ode23 и ode45 осуществляется следующим образом

[t,x]=ode23('fun', t0, tf, x0)

[t,x]=ode45('fun', t0, tf, x0)

или

[t,x]=ode23('fun', t0, tf, x0, tol, trace)

[t,x]=ode45('fun', t0, tf, x0, tol, trace)

где fun – имя m-файла, с помощью которого вычисляются правые части уравнений или системы дифференциальных уравнений, t0 – начальное значение аргумента, левая граница отрезка построения решения, tf – конечное значение аргумента, правая граница построения решения, х0 – вектор начальных условий, tol – задаваемая абсолютная погрешность построения решения, по умолчанию для ode23 она принята равной 0.001, а для ode45 – 0.00000l, trace – признак выдачи промежуточных результатов.

Пример решения на ПЭВМ в среде Matlab

Пусть дана задача Коши для дифференциального уравнений 1-го порядка

,

для которой надо построить решение на отрезке [0, 5] с относительной погрешностью не более 2%.

Для применения одной из описанных выше функций Matlab’а сначала надо сформировать m-файл, например с именем S.m, для вычисления правой части уравнения

function v=S(x,y)

v=sin(x*y);

После этого можно обратиться к функции ode45, задавая имя созданного m-файл, начальное условие, границы интервала интегрирования и требующуюся абсолютную погрешность интегрирования, например 10^–12. Относительная погрешность решения при этом не находится. Построим также график полученного решения. Описанные действия реализуются следующим набором команд

t0=0;

x0=[1];

tf=5;

[t,P]=ode45('S',t0,tf,x0,1e-12)

plot(t,P);

После запуска этой программы на выполнение в командном окне появятся результаты решения в виде таблицы

t = P =

0 1.0000

0.1250 1.0078

0.2500 1.0316

0.3750 1.0719

0.5000 1.1297

0.6250 1.2058

0.7500 1.3006

0.8750 1.4119

................

4.9600 0.6636

5.0000 0.6576

а в графическом окне (см. рис.6) будет показан график полученного решения.

Рис.6.

Для прямой реализации метода Эйлера и оценки его относительной погрешности по правилу Рунге следует написать две программы. В первой программе надо реализовать сам метод Эйлера и оформить её как функцию, в которой можно задавать шаг интегрирования, начальное условие и границу интервала. Текст этой программы может выглядеть следующим образом

function [y,x]=eiler(a,b,n,y0)

h=(b-a)/n

y(1)=y0; x(1)=a; ←задание начального условия

for i=2:n+1

x(i)=a+(i-1)*h;

y(i)=y(i-1)+h*sin(x(i-1)*y(i-1));

end

При такой организации процесса вычислений вторая программа должна реализовывать правило Рунге. Она должна дважды обратиться к первой программе. Первый раз с n шагами по отрезку интегрирования рассматриваемой задачи Коши, а второй раз – с 2 n шагами. Результаты такого обращения будут храниться в векторах y и y1 соответственно. Для расчета относительной погрешности решения, хранящегося в векторе y1, следует вычислить максимальное значение этого решения. После этого можно найти его относительную погрешность, сохранив её в виде вектора dy. Текст такой программы при 10-и шагах по отрезку построения решения может иметь вид

n=10

[y,x]=eiler(0,5,n,1);

[y1,x1]=eiler(0,5,2*n,1);

maxy(1:n)=y1(2:2:2*n+1);

plot(x,y,x1,y1), grid ←вывод графиков решений

for i=2:n+1

maxy(i)=max(y1(3:2:2*i-1));←поиск максимума решения

end

dy=abs(y(1:n+1)-y1(1:2:2*n+1))./maxy;

hold on ←использование того же графического окна

Y=[x',y1(1:2:2*n+1).',dy.'] ←построение матрицы решений и погрешности

plot(x,dy),grid

Эту программу следует сохранить в виде m-файла с именем odeqv.m в той же папке, что и функцию, реализующую метод Эйлера.

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

Y =

0 1.0000 0

0.5000 1.0619 0.0582

1.0000 1.3829 0.1035

1.5000 1.8520 0.0753

2.0000 1.8778 0.0563

..........................

5.0000 0.6555 0.0031

При этом в графическом окне отобразится рисунок с графиками решения с 10-ю шагами по отрезку интегрирования рассматриваемой задачи Коши, с 20-ю шагами и график относительной погрешности последнего решения. Он показан на рис.7. На этом рисунке можно видеть, что для 20-и шагов интегрирования заданная точность решения задачи Коши не достигнута. Наибольшая относительная погрешность превышает 10%. Поэтому расчёт надо повторить с числом шагов 20 и 40, 40 и 80 и так далее. Как показывают расчёты при 80 и 160 шагах получаемое решение имеет относительную погрешность при 160 шагах h = 0.03125, которая не превышает 1,6%, то есть удовлетворяет заданному условию. График полученного решения совпадает с графиком, показанным на рис.6.

При необходимости решения поставленной задачи Коши усовершенствованным методом Эйлера можно воспользоваться предыдущеей программой. Для этого надо переделать m-файл, в котором хранится программа решения задачи Коши методом Эйлера. Его содержимое, сохраняемое под именем u_eiler, показано ниже

function [y,x]=u_eiler(a,b,n,y0)

h=(b-a)/n

y(1)=y0; x(1)=a; ←задание начального условия

for i=2:n+1

xs(i)=x(i-1)+h/2; ←средняя точка

ys(i)=y(i-1)+h/2*sin(x(i-1)*y(i-1));

x(i)=a+(i-1)*h;

y(i)=y(i-1)+h*sin(xs(i)*ys(i));

end

При этом в программе odeqv.m следует только заменить имя функции eiler на u_eiler. Её запуск на выполнение для 10-шагов по отрезку интегрирования даёт в командном окне следующий результат, в графическом окне отобразится график, показнный на рис.8.

Y =

0 1.0000 0

0.5000 1.1282 0.0024

1.0000 1.5334 0.0027

1.5000 1.9068 0.0099

2.0000 1.8100 0.0263

...........................

5.0000 0.6595 0.0208

Рис.8.

Точность полученного решения – неудовлетворительна: относительная погрешность составляет 3%. Поэтому надо пересчитать решение с более мелким разбиением отрезка интегрирования. В данном случае оказывается, что достаточно разбить отрезок интегрирования на 40 участков, чтобы получить относительную погрешность решения не превышающую 1.1%. При этом график этого решения совпадёт с графиком, показанным на рис.6.

Дата добавления: 2015-07-08