r-k-f.cpp

较著名的解初值微分方程的数值方法——自适应Runge-Kutta-Fehlberg算法

#include<iostream.h>
#include<math.h>
double a,b,y,Y0,t,TOL,h,Hmax,Hmin,K1,K2,K3,K4,K5,K6,R,delta;
double f(double t,double y)
{return -y+t*t+1;}
void main()
{cout<<"输入端点a,b的值:"<<endl;
cin>>a>>b;
cout<<"输入误差容限TOL:"<<endl;
cin>>TOL;
cout<<"分别输入最大步长Hmax，最小步长Hmin:"<<endl;
cin>>Hmax>>Hmin;
cout<<"输入初值Y0:"<<endl;
cin>>Y0;
t=a;
y=Y0;
h=Hmax;
cout<<"("<<t<<","<<y<<")"<<'\n';
for(;h>=Hmin;){
if(t<b)
{K1=h*f(t,y);
K2=h*f(t+h/4,y+K1/4);
K3=h*f(t+3*h/8,y+3*K1/32+9*K2/32);
K4=h*f(t+12*h/13,y+1932*K1/2197-7200*K2/2197+7296*K3/2197);
K5=h*f(t+h,y+439*K1/216-8*K2+3680*K3/513-845*K4/4104);
K6=h*f(t+h/2,y-8*K1/27+2*K2-3544*K3/2565+1859*K4/4104-11*K5/40);
R=fabs(K1/360-128*K3/4275-2197*K4/75240+K5/50+2*K6/55)/h;
delta=0.84*pow(TOL/R,0.25);
if(R<=TOL)
{t=t+h;
y=y+25*K1/216+1408*K3/2565+2197*K4/4104-K5/5;
cout<<"("<<t<<","<<y<<","<<h<<")"<<'\n';
}
if(delta<=0.1)h=0.1*h;
else if(delta>=4.0)h=4*h;
else h=delta*h;
if(h>=Hmax)h=Hmax;
else cout<<"Minimum h exceeded!"<<'\n';
}}
cout<<"The procedure is complete"<<endl;
}