rungekutta.h

常微分方程数值解法

#ifndef RungeKutta_h
#define RungeKutta_h
#include "Differential.h"
#define RungeKuttaf "rungeainitial.txt"
template<class Type> class RungeKutta:public Differential<Type>{
	private:
		int m;Type a[6],b[6][6],c[6],k[6];
	protected:
		Type *u1,*e1;
	public:
		RungeKutta(Type (*pf)(Type,Type)):Differential<Type>
		(pf){}/*u1=new Type[N+1];e1=new Type[N+1];*/
		RungeKutta(Type (*pu)(Type),Type (*pf)(Type,Type));
		~RungeKutta(){delete []u1;delete []e1;}
		void Calculate();
		void Solution();
		void Error();
};
template<class Type> RungeKutta<Type>::RungeKutta(Type (*pu)(Type),Type (*pf)(Type,Type))
:Differential<Type>(pf){
	int i,j;Type d;
	u1=new Type[N+1];e1=new Type[N+1];	u=pu;
	cout<<"Input m:";
    //n=N-k+1;
	ifstream fin(RungeKuttaf,ios::in|ios::nocreate);
	if(!fin){
		cout<<RungeKuttaf<<" cannot be openned!"<<endl;exit(1);
	}
	fin>>m;
	for(i=1;i<m+1;i++)
	{fin>>a[i];}//cout<<a[i]<<endl;
	for(i=2;i<m+1;i++)
		for(j=1;j<i;j++)
		{fin>>b[i][j];}//cout<<b[i][j]<<endl;
		for(j=1;j<m+1;j++)
		{fin>>c[j];}//cout<<c[j]<<endl;
		fin>>d;//cout<<d;
		for(i=1;i<m+1;i++)
			c[i]=c[i]/d;
		fin.close();
}
template<class Type> void RungeKutta<Type>::Calculate(){
	int l,i,j;Type temp=0;
    cout<<"Input u("<<low<<"):";cin>>u1[0];
    for(l=1;l<N+1;l++){
	for(i=1;i<m+1;i++){
		for(j=1;j<i;j++)temp=temp+b[i][j]*k[j];
		k[i]=F(t[l-1]+h*a[i],u1[l-1]+h*temp,f);//cout<<k[i]<<"  ";
		temp=0;
	}
	for( i=1;i<m+1;i++)
		temp+=c[i]*k[i];
	u1[l]=u1[l-1]+h*temp;//cout<<u1[l]<<endl;
	temp=0;}
}
template<class Type> void RungeKutta<Type>::Error(){
	Type x=0;
	Exact();
    cout<<setw(6)<<"t"<<setw(14)<<"RungeKutta"<<setw(14)<<"Precision"
		<<setw(14)<<"Error"<<endl;
	for(int i=0;i<N+1;i++)
	{
		x=fabs(u0[i]-u1[i]);
		cout<<setw(6)<<t[i]<<setw(14)<<u1[i]<<setw(14)<<u0[i]<<
			setw(14)<<x<<endl;
	}
}
#endif