因研究卫星轨道的预报而需要研究常微分方程的数值解法，关于龙格库塔积分算法的证明公式可从代码的参考文献中入手。个人感觉龙格库塔积分的最大创新意义在于它实现了泰勒展开公式中高阶导数无法求解的困难，而是采用不同点的一阶导数替代，然后与泰勒展开公式对比系数实现！

1.四阶龙格库塔积分算法MATLAB代码

%--------------------------------------------------------------------------

%功能：            四阶龙格-库塔积分算法

%参考文献：      《GPS理论、算法与应用》—许国昌（GFZ）

%          作者：何成文       2016年5月30日

%I/O口说明：

%       t——时间节点！

%       X0——初始值!

%       F——求导方程右边的含t和X的函数!——外部函数调用！

%       h——每次积分的间隔，即h=X2-X1!

%       X——输出为龙格库塔积分的单步结果！

%--------------------------------------------------------------------------

function X=Runge_Kutta4(t,X0,h)

K1=h*F(t,X0);

K2=h*F(t+h/3,X0+K1/3);

K3=h*F(t+2*h/3,X0+K2-K1/3);

K4=h*F(t+h,X0+K1/2+K3/2);

X=X0+(K1+K2+K3+K4)/4;

%--------------------------------------------------------------------------

2.八阶龙格库塔算法MATLAB代码

%--------------------------------------------------------------------------

%功能：            8阶龙格-库塔积分算法

%参考文献：      《GPS理论、算法与应用》—许国昌（GFZ）

%          作者：何成文       2016年5月30日

%I/O口说明：

%       t——时间节点！

%       X0——初始值!

%       F——求导方程右边的含t和X的函数!——外部函数调用！

%       h——每次积分的间隔，即h=X2-X1!

%       X——输出为龙格库塔积分的单步结果！

%--------------------------------------------------------------------------

function X=Runge_Kutta8(t,X0,h)

K1=h*F(t,X0);

K2=h*F(t+h*4/27,X0+4*K1/27);

K3=h*F(t+2*h/9,X0+K1/18+K2/6);

K4=h*F(t+h/3,X0+K1/12+K3/4);

K5=h*F(t+h/2,X0+K1/8+3*K4/8);

K6=h*F(t+2*h/3,X0+(13*K1-27*K3+42*K4+8*K5)/54);

K7=h*F(t+h/6,X0+(389*K1-54*K3+966*K4-824*K5+243*K6)/4320);

K8=h*F(t+h/6,X0+(-231*K1+81*K3-1164*K4+656*K5-122*K6+800*K7)/20);

K9=h*F(t+5*h/6,X0+(-127*K1+18*K3-678*K4+456*K5-9*K6+576*K7+4*K8)/288);

K10=h*F(t+h,X0+(1481*K1-81*K3+7104*K4-3376*K5+72*K6-5040*K7-60*K8+720*K9)/820);

X=X0+(41*K1+27*K4+272*K5+27*K6+216*K7+216*K9+41*K10)/840;

%--------------------------------------------------------------------------

X=X0+(K1+K2+K3+K4)/4;

%--------------------------------------------------------------------------

3.算例1仿真结果曲线（仿真求导函数y'=cos(t),初值y(0)=0）

子函数：

function f=F(t,x)

f=cos(t);

仿真结果如下：

误差曲线：

结论：四阶龙格库塔比八阶龙格库塔算法效果好，误差低！

4.算例2（y'=1.0/sqrt(x*x+y*y),y(0)=0.1）

主函数：

%--------------------------------------------------------------------------

%功能：            验证龙格库塔积分法的正确性（y'=cos(t),y0=0)

%               作者：何成文     2016年5月30日

%--------------------------------------------------------------------------

clear;clc;

t=[0:0.01:40];

y_4=zeros(1,length(t));y_4(1,1)=0.1;

y_8=zeros(1,length(t));y_8(1,1)=0.1;

h=0.01;

for i=2:length(t)

   wuyon4=y_4(i-1);

   t0=t(i-1);

   y_4(i)=Runge_Kutta4(t0,wuyon4,h);

   wuyon8=y_8(i-1);

   y_8(i)=Runge_Kutta8(t0,wuyon8,h);

end

%--------------------------------------------------------------------------

figure;

subplot(1,3,1);plot(t,y_4,'ob');hold on;plot(t,y_4,'b');grid on;title('4阶龙格库塔积分结果');

subplot(1,3,2);plot(t,y_8,'og');hold on;plot(t,y_8,'b');grid on;title('8阶龙格库塔积分结果');

subplot(1,3,3);plot(t,y_8-y_4,'or');grid on;title('两种龙格库塔积分方法的误差曲线');xlabel('t');ylabel('误差');

子函数：

function f=F(t,x)

f=1.0/sqrt(t*t+x*x);

仿真图形：