Эйлера модифицированного (ЭМ) и Эйлера усовершенствованного (ЭУ).

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

Решаем уравнение методом Эйлера . Индексы i и i+1 означают, соответственно. моменты времени t и t+Δ t. Разложим x(t+Δ t) в ряд Тейлора в точке t:

(4.6)

Сравнивая алгоритм Эйлера и алгоритм ряда Тейлора, видим совпадение 2^х слагаемых, следовательно, порядок точности метода Эйлера 0(Δ t). В данном случае первое отброшенное слагаемое равно .

Запишем обобщенную формулу алгоритмов метода Рунге-Кутта в виде , где (4.7)

записана как полусумма производных в точках t и t+Δ t, т.е. в соответствии с методом ЭМ. При этом значение x (t+Δ t) получено методом Эйлера:

. Шаг сетки по времени равномерный: .

Используя обобщенный алгоритм запишем функцию в соответствии с методом ЭУ: . (4.8)

Производная взята в средней точке интервала t+Δ t, а значение x в этой точке получено методом Эйлера:

.

Представим слагаемые ряда (4.6) в следующем виде:

, (4.9)

где введены следующие обозначения: (x_i, t_i)=f, .

Чтобы оценить значения f для методов ЭМ и ЭУ разложим f(x, t) в ряд Тейлора:

. (4.10)

Используя вместо x и t значения x_i+1 и t_i+1 = t_i+Δ t, определим значение . И, подставив ее в (4.7), получим функцию для метода ЭУ: .

Подставим ее в (4.) и получим . Сравним полученное выражение с разложением в ряд Тейлора (4.9). Совпадают три слагаемых ряда Тейлора.

Аналогичным образом получим функцию для метода ЭМ. , Подставив ее в обобщенный алгоритм, увидим, что имеется совпадение с тремя слагаемыми ряда Тейлора, поэтому точность методов ЭМ и ЭУ имеет порядок 0(Δ t²).

Приведем программы методов ЭМ и ЭУ на MATLAB:

function EylerY

x=0; t=0;

display(x)

display(t)

dt=1;

for m=1: 6

dx=FunkEMY(x);

x1=x+dx*dt;

dx1=FunkEMY(x1);

x=x+0.5*(dx+dx1)*dt

t=t+dt

end

function dx=FunkEMY(x)

k=1; u=1; Tau=1;

dx=0;

dx=(k*u-x)/Tau;

function EylerM

x=0; t=0;

display(x)

display(t)

dt=1;

for m=1: 6

dx=FunkEMY(x);

x1=x+0.5*dx*dt;

dx1=FunkEMY(x1);

x=x+dx1*dt

t=t+dt

end