高阶分方程组初值问题的数值解法.txt

工程问题中常用微分方程的形式建立数学模型

#include<iostream.h>
#include<iomanip.h>
#include<math.h>
double fun1(double y,double z)
{
	double r;
	r=z;
	return r;
}
double fun2(double y,double z)
{
	double r;
	r=-2*y+z;
	return r;
}
void Euler(double yn,double zn,double h,int N)
{
	double K,L;
	int i;
	cout<<"Euler公式计算得："<<endl;
	for(i=1;i<=N;i++)
	{
		K=fun1(yn,zn);
		L=fun2(yn,zn);
		yn+=h*K;
		zn+=h*L;
		cout<<"y"<<i<<"="<<yn<<setw(10)<<"z"<<i<<"="<<zn<<endl;
	}
	cout<<endl;
}
void A_Euler(double yn,double zn,double h,int N)
{
	double K1,K2,L1,L2;
	int i;
	cout<<"改进Euler公式计算得："<<endl;
	for(i=1;i<=N;i++)
	{
		K1=fun1(yn,zn);
		L1=fun2(yn,zn);
		K2=fun1(yn+h*K1,zn+h*L1);
		L2=fun2(yn+h*K1,zn+h*L1);
		yn+=h/2*(K1+K2);
		zn+=h/2*(L1+L2);
		cout<<"y"<<i<<"="<<yn<<setw(10)<<"z"<<i<<"="<<zn<<endl;
	}
	cout<<endl;
}	
void jd_R_K(double yn,double zn,double h,int N)//经典公式
{
	double K1,K2,K3,K4,L1,L2,L3,L4;
	int i;
	cout<<"经典公式计算得："<<endl;
	for(i=1;i<=N;i++)
	{
		K1=fun1(yn,zn);
		L1=fun2(yn,zn);
		K2=fun1(yn+h/2*K1,zn+h/2*L1);
		L2=fun2(yn+h/2*K1,zn+h/2*L1);

		K3=fun1(yn+h/2*K2,zn+h/2*L2);
		L3=fun2(yn+h/2*K2,zn+h/2*L2);
		
		K4=fun1(yn+h*K3,zn+h*L3);
		L4=fun2(yn+h*K3,zn+h*L3);
		
		yn+=h/6*(K1+2*K2+2*K3+K4);
		zn+=h/6*(L1+2*L2+2*L3+L4);
		cout<<"y"<<i<<"="<<yn<<setw(10)<<"z"<<i<<"="<<zn<<endl;
	}
	cout<<endl;
}
void main()
{
    double a=0,b=1;//定义域；
	double y0=1;
	double z0=1;
	double h=0.1;
	int N=int((b-a)/h);
	int t=0;
	cout<<"高阶分方程组初值问题的数值解法："<<endl;
	cout<<"y''-3y'+2y=0"<<endl;// x<-["<<a<<","<<b<<"],"<<endl;
	cout<<"y(0)=1,y'(0)=1."<<endl<<endl;
    cout<<"请问使用哪种方法求解？"<<endl;
	cout<<"Euler公式，请输入1"<<endl;
	cout<<"改进Euler公式，请输入2"<<endl;
	cout<<"四阶R-K经典公式，请输入3"<<endl;
	cin>>t;
	if(t==1)
	{
		Euler(y0,z0,h,N);//Euler公式很不精确
		cout<<"太不精确了，还是改用其他方法吧！"<<endl;
	}
	else if(t==2)
	{
		A_Euler(y0,z0,h,N);
	}
	else if(t==3)
	{
		jd_R_K(y0,z0,h,N);
	}
	else 
	{
		cout<<"不进行求解。"<<endl;
	}
}