adamsin.h

常微分方程数值解法

#ifndef Adamsin_h
#define Adamsin_h
#include "Adamsout.h"
#define Adamsinf "adamsinitial.txt"
template<class Type> class Adamsin:public Adamsout<Type>{
	private:
		Type b[6][6]; Type *u2;
	public:
		Adamsin(Type (*pf)(Type,Type)):Adamsout<Type>(pf){}
		/*{u1=new Type[N+1];e1=new Type[N+1];}*/
		Adamsin(Type (*pu)(Type),Type (*pf)(Type,Type));
		~Adamsin(){delete []u2;}//{delete []u1;delete []e1;}
		void Calculate();
		void Solution();
		void Error();
};
template<class Type> Adamsin<Type>::Adamsin(Type (*pu)(Type),Type (*pf)(Type,Type))
:Adamsout<Type>(pu,pf){
	u2=new Type[N+1];
	int i,j; 
	ifstream fin(Adamsinf,ios::in|ios::nocreate);
	if(!fin){
		cout<<Adamsinf<<" cannot be openned!"<<endl;exit(1);
	}
	for(i=0;i<6;i++)
		for(j=0;j<i+1;j++)
		{fin>>b[i][j];}
		for(j=0;j<6;j++){
			b[1][j]=b[1][j]/2;
			b[2][j]=b[2][j]/12;
			b[3][j]=b[3][j]/24;
			b[4][j]=b[4][j]/720;
			b[5][j]=b[5][j]/1440;
		}
		fin.close();
}
template<class Type> void Adamsin<Type>::Calculate(){
	Adamsout<Type>::Calculate();
	int i,l;Type temp=0;
	for(i=0;i<N+1;i++)
		u2[i]=u1[i];
	/*for( i=1;i<=k;i++)
		u1[i]=u1[i-1]+h*F(t[i-1],u1[i-1],f);*/
	for(i=k+2;i<N+1;i++){
		for(l=1;l<=k+1;l++)
		{temp=temp+b[k+1][l]*F(t[i-l],u2[i-l],f);
		}//cout<<temp<<endl;
		//for(long x=0;x<10000;x++)
		u2[i]=u2[i-1]+h*b[k+1][0]*F(t[i],u2[i],f)+h*temp;//cout<<u1[i]<<endl;
		temp=0;
	}
}
template<class Type> void Adamsin<Type>::Error(){
	Type x=0;
	Exact();
    cout<<setw(6)<<"t"<<setw(14)<<"Adamsin"<<setw(14)<<"Precision"
		<<setw(14)<<"Error"<<endl;
	for(int i=0;i<N+1;i++)
	{
		x=fabs(u0[i]-u2[i]);
		cout<<setw(6)<<t[i]<<setw(14)<<u2[i]<<setw(14)<<u0[i]<<
			setw(14)<<x<<endl;
	}
}
#endif