MATLAB — преобразование метода Эйлера в метод Хойна

Я видел этот метод Эйлера, который вычисляет ошибки:

function [t,le,ge] = euler_errors(h)

f=@(u) u*(2-u); % this is the function for the IVP
t0=0;
tn=1;
t=t0:h:tn;%we want to find the errors along this solutions
%here is the exact solution of the IVP
u_exact=(0.2*exp(2*t))./(2+0.1*(exp(2*t)+1)); %the with initial value u(0)=0.1

n=length(t);
u_e=zeros(1,n);
u_g=zeros(1,n);

i=1;
u_e(i)=0.1;
u_g(i)=0.1;

%u_e and u_g are both values given by Euler method
%u_e is for the local error
%u_g is for the global error
while (i<n)
    u_e(i+1)=u_e(i)+h*f(u_e(i));
    u_g(i+1)=u_g(i)+h*f(u_exact(i));
    i=i+1;
end;
%le1 is the local error
%ge1 is the global error
le=abs(u_e-u_exact);
ge=abs(u_g-u_exact);
end

Я попытался преобразовать метод в метод Хойна — вот моя попытка:

function [t,le,ge] = heun_errors(h)

f=@(u) u*(2-u); % this is the function for the IVP
t0=0;
tn=1;
t=t0:h:tn;%we want to find the errors along this solutions
%here is the exact solution of the IVP
u_exact=(0.2*exp(2*t))./(2+0.1*(exp(2*t)+1)); %the with initial value u(0)=0.1

n=length(t);
u_e=zeros(1,n);
u_g=zeros(1,n);

i=1;
u_e(i)=0.1;
u_g(i)=0.1;

%u_e and u_g are both values given by Euler method
%u_e is for the local error
%u_g is for the global error
while (i<n)
    u_e1(i+1)=u_e(i)+h*f(u_e(i));
u_e(i+1)=u_e(i)+(h/2)*(f(u_e(i))+f(u_e1(i+1)));
u_g1(i+1)=u_g(i)+h*f(u_exact(i));
u_g(i+1)=u_g(i)+(h/2)*(f(u_exact(i))+f(u_g1(i+1)));

    i=i+1;
end;
%le1 is the local error
%ge1 is the global error
le=abs(u_e-u_exact);
ge=abs(u_g-u_exact);
end

Но теперь ошибка действительно увеличилась. Может кто-нибудь сказать мне, что я сделал неправильно и как я могу это исправить?


matlab numerical-methods differential-equations runge-kutta
person peter    schedule 10.12.2015    source источник
comment
Ваш код правильный, насколько я могу судить. Возможно, по совпадению именно этот ОДУ лучше работает с методом Эйлера. Как вы сравниваете ошибку? Лучший способ проверить это — определить сходимость каждого метода по мере уменьшения h. Метод Эйлера должен сходиться линейно, но метод Хойна должен улучшаться быстрее (квадратично).   -  person David    schedule 10.12.2015
comment
Имеет ли метод Хеуна большую ошибку усечения, потому что в расчетах больше частей? два шага вместо одного в методе Эйлера?   -  person peter    schedule 10.12.2015
comment
Нет, у него должна быть меньшая ошибка усечения, но большая ошибка округления, я думаю.   -  person David    schedule 11.12.2015

Ответы (1)


arrow_upward
0
arrow_downward

Исправление точного решения

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

(e^(2t)/u)'=e^(2t)*(-u'/u^2 + 2/u) = e^(2t)

e^(2t)/u(t)-1/u0 = (e^(2t)-1)/2  

u(t) = 2*u0*exp(2*t) ./ ( 2 + u0*(exp(2*t)-1) )

О вычислении локальных и глобальных ошибок

В организации кода Эйлера для вычисления локальной и глобальной ошибки уже есть что-то концептуально перепутанное. Полное нормальное интегрирование дает глобальный результат и, следовательно, также глобальную ошибку в отличие от точного решения. Таким образом, то, что называется u_e, должно быть u_g.

Локальная ошибка сравнивает точное решение u_i(t), начатое в (t_i,u_i), с шагом метода из той же точки. Поскольку в игре есть только одно точное решение, он должен выполнить шаг для (t_i,u_i) = (t(i), u_exact(i)), чтобы затем сравнить шаг метода, результат которого сохраняется в u_e(i+1), с точным значением u_exact(i+1).

Обе выданные адресуются с изменением цикла на

while (i<n)
  u_g(i+1)=u_g(i)+h*f(u_g(i));
  u_e(i+1)=u_exact(i)+h*f(u_exact(i));
  i=i+1;
end;

Точно так же эти проблемы также должны быть решены в коде Heun, где цикл, таким образом, изменяется на (с использованием промежуточных переменных для этапов, как в общем методе Рунге-Кутты)

while (i<n)
  k1=f(u_g(i));
  k2=f(u_g(i)+h*k1);
  u_g(i+1)=u_g(i)+(h/2)*(k1+k2);
  k1=f(u_exact(i));
  k2=f(u_exact(i)+h*k1);
  u_e(i+1)=u_exact(i)+(h/2)*(k1+k2);
  i=i+1;
end;

Графики профилей ошибок

Затем построение профилей ошибок le/h^(p+1) локальной и lg/h^p глобальной ошибки во времени с порядками ошибок p=1 для Эйлера и p=2 для Гойна и наложение точечных диаграмм с различными размерами шага дает

графики профиля ошибки

которые визуально инвариантны относительно размера шага.

В качестве проверки на непротиворечивость локальная ошибка Эйлера составляет около 15_, а константа Липшица ОДУ около 2 в соответствующей области, поэтому глобальная ошибка соответствует лемме Грёнвалла примерно закону (exp(2*t)-1)*0.15*h для малых t. Точно так же метод Heun имеет локальную ошибку около 0.15*h^3, что дает глобальную ошибку (exp(2*t)-1)*0.075*h^2 для небольшого t, что совместимо с глобальным графиком ошибок.

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

person Lutz Lehmann    schedule 01.01.2019