【Runge-Kutta】龙格库塔不同步长积分到终点不一样

参考链接：
1、https://blog.csdn.net/xiaokun19870825/article/details/78763739

解析式和微分式

疑问步长不一样最终结果相同吗？
对于解析式y=sqrt(1+2x)，可以写成下面常微分形式：

使用RK4（ode45）

下面是自己编写的matlab代码和ode45：
test_fun.m

function dy=test_fun(x,y)
dy = zeros(1,1);%init data
dy(1) = y - 2*x/y;
% dy(1) = 1./y;

runge_kutta1.m

function [x,y]=runge_kutta1(ufunc,y0,h,a,b)
% The order of the parameter list is the function name, 
% initial value vector, step size, time starting point and time 
% ending point of the differential equation system 
% (parameter form refers to ode45 function)

n=floor((b-a)/h);%Step number
x(1)=a;%Starting point of time
y(:,1)=y0;%Initial value can be vector, but dimension should be paid attention to.
for ii=1:n
x(ii+1)=x(ii)+h;

k1=ufunc(x(ii),y(:,ii));
k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);
k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);
k4=ufunc(x(ii)+h,y(:,ii)+h*k3);

y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;
%Numerical solution based on Runge Kutta method
end

运行脚本run_main.m

[T,F] = ode45(@test_fun1,[0 1],[1]);
subplot(131)
plot(T,F)%Matlab's own ode45 function results
title('Ode45 function effect')

[T1,F1]=runge_kutta1(@test_fun1,[1],0.25,0,1);%When testing, change the dimension of test_fun function. Don't forget to change the dimension of the initial value
subplot(132)
plot(T1,F1)%Self compiled function of Runge Kutta function
title('Self compiled Runge Kutta function,step: 0.25')

[T2,F2]=runge_kutta1(@test_fun1,[1],0.1,0,1);%When testing, change the dimension of test_fun function. Don't forget to change the dimension of the initial value
subplot(133)
plot(T2,F2)%Self compiled function of Runge Kutta function
title('Self compiled Runge Kutta function,step: 0.01')

从0到1步长不一样，得到的结果不一样

ode45：内部自动插值步长0.025，计算出来y(1) = 1.732050816766531
RK4: 步长：0.250，计算出来y(1) = 1.732274190526737
RK4: 步长：0.100， 计算出来y(1) = 1.732056365165567
RK4: 步长：0.025， 计算出来y(1) = 1.732050828604835
RK4: 步长：0.010， 计算出来y(1) = 1.732050808103274 (精度最高)
真实y(1) = 1.732050807568877
结论：针对不含误差的方程而言：从0积分到1步长越小积分越准确。有误差的微分方程还未做实验。

计算图像大致一样