7-8阶变步长Runge-Kutta方法（Matlab版)
Sun, 11/25/2007 - 22:33 — kambrian
从这里下载该文件
预览：
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
- Add new comment
- 1316 reads
- quote
在这里有从rk2到rk7的所有例子
在这里有从rk2到rk7的所有例子
alljoyland
http://google.com/codesearch?hl=en&q=lang:matlab+rk7+runge+kutta+show:AboMCR3FyR4:iTO94XLDbY4:9grk-9bP3u4&sa=N&cd=1&ct=rc&cs_p=http://www.soe.ucsc.edu/~hongwang/Codes/ODE_RK/ODE_RK.zip&cs_f=ODE_RK/calc_RK7.m#l-2
能否贴出fortan77的7-8阶定步长Rk method, 感激中...
3Q3Q
其实这种经典的算法
这种经典的算法网上有很多成熟的软件包，这里给出两个：【RKSE_tar.gz】【RKsuite_f90.zip】。
这两个包里都是各阶的都有，第一个也许更明了一些，希望会有帮助。
如果自用的话不如自己编一个，比那些通用的包会简明很多。基本算法可以参考上面我补上去的课件链接。 定步长的更简单。
能否贴出7-8阶变步长Ru
能否贴出7-8阶变步长Runge-Kutta方法的c语言版 这里有急用！ 或者由 duanjf88@yahoo.com.cn发给我
C语言版看这里
//h-步长，t、x-自变量和函数，err-积分的tolerance，n-要积分的函数个数。
double RKF78(double h,double t,double *x,double err,int n)
{
double a[13]={0,2./27,1./9,1./6,5./12,1./2,5./6,1./6,2./3,1./3,1,0,1};
double c1[11]={41./840,0,0,0,0,34./105,9./35,9./35,9./280,9./280,41./840};
double b[13][12]
={{0,0,0,0,0,0,0,0,0,0,0,0,},
{2./27,0,0,0,0,0,0,0,0,0,0,0},
{1./36,1./12,0,0,0,0,0,0,0,0,0,0,},
{1./24,0,1./8,0,0,0,0,0,0,0,0,0,},
{5./12,0,-25./16,25./16,0,0,0,0,0,0,0,0},
{1./20,0,0,1./4,1./5,0,0,0,0,0,0,0},
{-25./108,0,0,125./108,-65./27,125./54,0,0,0,0,0,0},
{31./300,0,0,0,61./225,-2./9,13./900,0,0,0,0,0},
{2,0,0,-53./6,704./45,-107./9,67./90,3,0,0,0,0},
{-91./108,0,0,23./108,-976./135,311./54,-19./60,17./6,-1./12,0,0,0},
{2383./4100,0,0,-341./164,4496./1025,-301./82,2133./4100,45./82,45./164,18./41,0,0},
{3./205,0,0,0,0,-6./41,-3./205,-3./41,3./41,6./41,0,0},
{-1777./4100,0,0,-341./164,4496./1025,-289./82,2193./4100,51./82,33./164,12./41,0,1}};
double c2[13]={0,0,0,0,0,34./105,9./35,9./35,9./280,9./280,0,41./840,41./840};
int i,j,m;
double k[6][13];
double x1[6],x2[6];
double delta;
double sum=0;
double temp[6];
for (m=0;m<6;m++)
{
x1[m]=x[m];
x2[m]=x[m];
}
for (m=0;m<6;m++)
{
temp[m]=x[m];
}
for (m=0;m<6;m++)
{
f(t,temp);
k[m][0]=temp[m];
}
for (m=0;m<6;m++)
{
for (i=0;i<13;i++)
{
if(i>0)
{
for (j=0;j<i;j++)
{
sum+=b[i][j]*k[m][j];
}
for (j=0;j<6;j++)
{
temp[j]=x[j]+h*sum;
}
f(t+a[i]*h,temp);
k[m][i]=temp[m];
sum=0;
for (j=0;j<6;j++)
{
temp[j]=x[j];
}
}
}
}
for (i=0;i<6;i++)
{
for (j=0;j<11;j++)
x1[i]+=h*c1[j]*k[i][j];
}
for (i=0;i<6;i++)
{
for (j=0;j<13;j++)
{x2[i]+=(h*c2[j]*k[i][j]);}
}
for (i=0;i<6;i++)
{
delta=fabs(x2[i]-x1[i]);
if (delta>err)
{
h=h*0.999;
flag=-1;
}
}
if (flag==1)
{
for (m=0;m<6;m++)
x[m]=x2[m];
}
return h;
}