Runge-Kutta 四阶方法没有给出预期的错误

Question

我编写了一个程序来比较 Euler 和 RK4 方法在以已知角度和速度抛出球的简单问题中的唯一重力。我计算了绝对误差，并用 x 轴的时间步长 (h) 绘制了它。 Euler 的绝对误差结果是正确的（与 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)

图表： absolute error with time-step

Answer 1

没有什么好奇怪的，您将常数函数（两次）与 4 阶方法相结合，得到的线性函数和二次函数正好在浮点运算的范围内。

为了得到一个非平凡的方法错误，你需要构建一个不像当前例子那样简单积分的例子，最简单的是 x'=x，下一个标准例子是谐振子x''+x=0。

您也可以通过空气摩擦进行弹道射击，但是没有可以比较的精确解决方案。您将需要更精确的数值解，例如以一半的步长再次积分。

Runge-Kutta 四阶方法没有给出预期的错误

Runge-Kutta 4th order method doesn't give the expected error

runge-kutta