參考連結：
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步長越小積分越準確。有誤差的微分方程還未做實驗。

計算影象大致一樣