Astronomical Software Club

从这里下载该文件

预览：

 
function [X,T]=myRKF7(x,t,rfun,tol,stepmin,stepmax,tmax,ani_flag)
 
% ****help info for myRKF7adv *****
 
% 
 
% [X,T]=myRKF7(x,t,rfun,tol,stepmin,stepmax,tmax,ani_flag)
 
% solving ODEs with RKF78 method (adaptive Runge-Kutta-Fehlberg 7/8th order),精确终点
 
%  ODEs:       dX/dt=rfun(X,T)
 
% x,t: initial vectors of X and T
 
% tol: error tolerence per unit time step for each dimension
 
% stepmin: minimum integral step for t
 
% stepmax: maximum integral step permitted
 
% tmax: maximum integral time 
 
% ani_flag: animation on/off,specified with the value 0/1. 
%If ani_flag=1, it will open a figure to plot the integration process
 
%
 
% demo parameters:
 
%   x=[1 0 0.2,0.2 0.7 0.2],t=0,h=0.05,
 
%   rfun=@sin
 
%   step=0.01,stepmax=1
 
%   tol=1e-8
 
%   tmax=10
 
%   ani_flag=1
%
% Author: 
% Han Jiaxin, <a href="mailto:hanjiaxin@gmail.com">hanjiaxin@gmail.com</a> 
% written 2005 in NJU,Nanjing
% http://asc.2dark.org

关于本程序中变步长算法的说明：

变步长算法：
1)计算ki及误差率delta, 令del=tolerance/delta
2)误差率小于容限时计算x'=x+dx
3)当delta<0.5tolerence,即del>2时将步长调大
4)如果误差率大于容限，步长调小

因方法误差~O（h^8）,则单位步长误差
del~O(h^7)=C*h^7，
为将误差率控制在容限处，取步长h'使得
del'=C*h'^7=tol
则
h'/h=(del'/del)^(1/7)=(tol/del)^(1/7)
可得
h'=h*(tol/del)^(1/7)
本程序中调整h使得del'=0.8*tol

龙格－库塔法的基础可参考南大周济林老师的课件：http://astronomy.nju.edu.cn/~zjl/cm/appdix-3.pdf