Метод Рунге-Кутты 4-го порядка не дает ожидаемой ошибки

Я написал программу для сравнения метода Эйлера и RK4 в простой задаче о мяче, брошенном под известным углом и скоростью с единственной силой тяжести. Я вычислил абсолютную ошибку и построил график с шагом (h) по оси x. Абсолютная ошибка Эйлера оказывается верной (пропорциональна h ^ 2), но абсолютная ошибка RK4 дает значения от 0 до 10 ^ -15. Однако я думаю, что ожидаемый результат пропорционален h ^ 4. Вот код:

hr=0.1:0.05:1;h=0.1/4:0.05/4:1/4; %Chosen step for each method
d1=zeros(1,19);d2=zeros(1,19); %difference between calculated and actual value for the two methods
ne=0;nrk=0;
for j=1:19%j counting the number of the step-sizes (19 different step-sizes)
n=2500*hr(j);
x=zeros(1,n);xrunge=zeros(1,n);
y=zeros(1,n);yrunge=zeros(1,n);
uy=zeros(1,n);uyrunge=zeros(1,n);
t=zeros(1,n);tr=zeros(1,n);
%Initializing
uy(1)=20*sin(pi/3);ux=20*cos(pi/3);
uyrunge(1)=20*sin(pi/3);
a=-10;
%Theoritical equations of motion for position and velocity to the y axis
function [y]=f(t)
    y=20*sin(pi/3).*t-5*t.^2;
end
%Euler
for i=2:n  
    if y(i-1)+h(j)*uy(i-1)<0
        break;
    end
    if h(j)==0.1
        ne=ne+1;%number of calculations for uy with step-size 0.1
    end
    t(i)=t(i-1)+h(j);
    y(i)=y(i-1)+h(j)*uy(i-1);
    x(i)=x(i-1)+h(j)*ux;
    uy(i)=uy(i-1)+h(j)*a;  
    if i==2
        d1(j)=abs(y(i)-f(t(i))); %number of operations done to calculate y with step-size 0.1
    end
end
%RK4
for i=2:n
    tr(i)=tr(i-1)+hr(j);
    xrunge(i)=xrunge(i-1)+hr(j)*ux;
    uyrunge(i)=uyrunge(i-1)+hr(j)*a;
    k1y=uyrunge(i-1);
    k2y=uyrunge(i-1)+hr(j)*a/2;
    k3y=uyrunge(i-1)+hr(j)*a/2;
    k4y=uyrunge(i-1)+hr(j)*a;
    if  yrunge(i-1)+hr(j)*(k1y+2*k2y+2*k3y+k4y)/6<0
        break;
    end
    yrunge(i)=yrunge(i-1)+hr(j)*(k1y+2*k2y+2*k3y+k4y)/6;
    if i==2
        d2(j)=abs(yrunge(i)-f(tr(i))); %difference at each step-size for the first iteration
    end
end
end
%function to fit |y-f(t)| according to least squares
function [d]=f2(h)
d=5*h.^2-5.621*10^(-7)*h+3.125*10^(-8);
end
figure(3)
plot(h,d1,'r.',hr,d2,'k.',h,f2(h))
xlabel('h')
legend('Euler', 'RK4','Least squares fitting')
ylabel( 'y-f(t)')
title('Accuracy with step')
disp(d2)

И график: абсолютная ошибка с шагом времени


runge-kutta
person JimmyNik    schedule 13.01.2019    source источник

Ответы (1)


arrow_upward
0
arrow_downward

В этом нет ничего удивительного, вы интегрируете постоянную функцию (дважды) с методом порядка 4, полученные линейные и квадратичные функции точны в пределах арифметики с плавающей запятой.

Чтобы получить нетривиальную ошибку метода, вам нужно будет построить пример, который не является простым интегрированием, как текущий, самый простой - x'=x, следующий стандартный пример - гармонический осциллятор x''+x=0.

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

person Lutz Lehmann    schedule 13.01.2019