ex6.cpp

方程的数值计算方法

#include<iostream.h>
#include<math.h>
#define N 100
double exactF(double x)
{
	return 3/(1+(pow(x,3)));
}
double Fnctn(double x,double y)
{
	double f=-1*(pow(x,2)*pow(y,2));
	return f;
}
void Runge(double x[N],double y[N],int Number,double h)
{
	double k1;
	double k2;
	double k3;
	double k4;
	for(int i=0;i<Number;i++)
	{
		k1=Fnctn(x[i],y[i]);
		k2=Fnctn((x[i]+h/2),(y[i]+1/2*h*k1));
		k3=Fnctn((x[i]+h/2),(y[i]+1/2*h*k2));
		k4=Fnctn((x[i]+h),(y[i]+h*k3));
		y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6;
	}
}
//AB4公式
void CalAB(double x[N],double y[N],int Number,double h)
{
	for(int i=2;i<Number;i++)
	{
		y[i+1]=y[i]+h*(55*Fnctn(x[i],y[i])-59*Fnctn(x[i-1],y[i-1])+37*Fnctn(x[i-2],y[i-2])-9*Fnctn(x[i-3],y[i-3]))/24;
	}
}
//AB4-AM4公式
void CalAb4Am4(double x[N],double y[N],int Number,double h)
{
	double k;
	for(int i=2;i<Number;i++)
	{
		k		=y[i]+h*(55*Fnctn(x[i],y[i])-59*Fnctn(x[i-1],y[i-1])+37*Fnctn(x[i-2],y[i-2])-9*Fnctn(x[i-3],y[i-3]))/24;
		//AB4-AM4公式
		y[i+1]	=y[i]+h*(9*Fnctn(x[i+1],k)+19*Fnctn(x[i],y[i])-5*Fnctn(x[i-1],y[i-1])+Fnctn(x[i-2],y[i-2]))/24;
	}
}
void CalAdvancedAb4Am4(double x[N],double y[N],int Number,double h)
{
	double k;
	double q;
	for(int i=2;i<Number;i++)
	{
		k		=y[i]+h*(55*Fnctn(x[i],y[i])-59*Fnctn(x[i-1],y[i-1])+37*Fnctn(x[i-2],y[i-2])-9*Fnctn(x[i-3],y[i-3]))/24;
		//AB4-AM4公式
		q	=y[i]+h*(9*Fnctn(x[i+1],k)+19*Fnctn(x[i],y[i])-5*Fnctn(x[i-1],y[i-1])+Fnctn(x[i-2],y[i-2]))/24;
		y[i+1]	=(k*251+q*19)/270;
	}

}

void main()
{
	double hLength;
	double a;
	double b;
	double x[N]				={0};	//自变量
	double RungeKutta[N]	={0};	//Runge=Kutta方法数值解
	double AB4[N]			={0};	//AB4方法数值解
	double ABAM4[N]			={0};	//AB4-AM4方法数值解
	double AdvancedAbAm4[N]	={0};	//改进的AB4-Am4数值解
	double y[N];					//精确解，用来与以上数值解比较
	int Number;
	cout<<"********本程序实现常微分方程的数值解法********\n";
	cout<<"请输入自变量的最小值：";
	cin>>a;
	cout<<"请输入自变量的最大值：";
	cin>>b;
	cout<<"请输入步长h：";
	cin>>hLength;
	cout<<"请输入初值y("<<a<<"):";
	cin>>RungeKutta[0];
	Number=(b-a)/hLength;
	for(int i=0;i<=Number;i++)
	{
		x[i]=a+i*hLength;
	}
	Runge(x,RungeKutta,Number,hLength);
	for(int j=0;j<4;j++)
	{
		AB4[j]=RungeKutta[j];
		ABAM4[j]=RungeKutta[j];
		AdvancedAbAm4[j]=RungeKutta[j];
	}
	CalAB(x,AB4,Number,hLength);
	CalAb4Am4(x,ABAM4,Number,hLength);
	CalAdvancedAb4Am4(x,AdvancedAbAm4,Number,hLength);
	for(i=0;i<=Number;i++)
		y[i]=exactF(x[i]);
	//输出结果
	cout<<'\t'<<"精确解"<<'\t'<<'\t'<<"数值解"<<'\t'<<'\t'<<"\t与精确值之差\n";
	//Runge-Kutta方法的结果及其误差
	cout<<"Runge-Kutta方法:\n";
	for(i=0;i<=Number;i++)
	{
		cout<<"y("<<x[i]<<")="<<'\t'<<y[i]<<'\t'<<'\t';
		if(i<14)
			cout<<'\t';						//调整格式
		cout<<RungeKutta[i]<<'\t'<<'\t';
		cout<<fabs(y[i]-RungeKutta[i])<<'\n';
	}
	cout<<"AB4方法:\n";
	for(i=0;i<=Number;i++)
	{
		cout<<"y("<<x[i]<<")="<<'\t'<<y[i]<<'\t'<<'\t';
		if(i<14)
			cout<<'\t';						//调整格式
		cout<<AB4[i]<<'\t'<<'\t';
		cout<<fabs(y[i]-AB4[i])<<'\n';
	}
	cout<<"AB4-Am4方法:\n";
	for(i=0;i<=Number;i++)
	{
		cout<<"y("<<x[i]<<")="<<'\t'<<y[i]<<'\t'<<'\t';
		if(i<14)
			cout<<'\t';						//调整格式
		cout<<ABAM4[i]<<'\t'<<'\t';
		cout<<fabs(y[i]-ABAM4[i])<<'\n';
	}
	cout<<"改进的AB4-Am4方法:\n";
	for(i=0;i<=Number;i++)
	{
		cout<<"y("<<x[i]<<")="<<'\t'<<y[i]<<'\t';
		if(i<14)
			cout<<'\t';						//调整格式
		cout<<AdvancedAbAm4[i]<<'\t'<<'\t';
		cout<<fabs(y[i]-AdvancedAbAm4[i])<<'\n';
	}

}