controlsystemdigitalsimuation_1.cpp

数值积分方法类型Model,1=欧拉法,2=二阶龙格库塔,4=四阶龙格库塔

	// ControlSystemDigitalSimuation_1.cpp : 定义控制台应用程序的入口点。
	//

	#include "stdafx.h"
	#include "math.h"
	#include "iostream"

	double u=1.0;
	double f[3];
	using namespace std;

	void Model(double f[3],double x[3]);
	void ModelRgkt_4(double x[], int k, double h);
	void ModelRgkt_2(double x[],int k,double h);
	void Modelola(double x[],int k,double h);
	int _tmain(int argc, _TCHAR* argv[])
	{ 
		double y;
		double u=1.0;
		double h;
		double x[]={0,0,0};
		int k=3;
		int Model;
		int count=60;
		double t=0.0;
		cout<<"输入数值积分的步长h"<< endl;
		cin>>h;
		cout<<"输入数值积分方法类型Model,1=欧拉法,2=二阶龙格库塔,4=四阶龙格库塔"<<endl;
		cin>>Model;
		FILE *fp;

		if(Model==1)

		{

			if(!(fp=fopen("Modelola_data.txt","w")))
				printf("creat file failed!");   
			for(int i=0;i<count;i++) 
			{
				t=t+h;
				Modelola(x,k,h);                      
				y=333.3333*x[1]+66.6667*x[2];
				fprintf(fp,"%f\t%f\t\n",i*h,y);

			}
		}

		if(Model==2)
		{

			if(!(fp=fopen("D:\\Java\\federate\\ControlSystemDigitalSimuation_1\\Rgkt2_data.txt","w")))
				printf("creat file failed!");   
			for(int i=0;i<count;i++) 
			{
				t=t+h;
				ModelRgkt_4(x,k,h);                      
				y=333.3333*x[1]+66.6667*x[2];
				fprintf(fp,"%f\t%f\t\n",i*h,y);

			}
		}


		if (Model==4)
		{


			if(!(fp=fopen("D:\\Java\\federate\\homework\\ControlSystemDigitalSimuation_1\\Rgkt_4.txt","w")))
				printf("creat file failed!");   
			for(int i=0;i<count;i++) 
			{
				t=t+h;
				ModelRgkt_4(x,k,h);                      
				y=333.3333*x[1]+66.6667*x[2];
				fprintf(fp,"%f\t%f\t\n",i*h,y);

			}
		}


		return 0;
	}                                           

	void Model(double f[3],double x[3])

	{   double x1,x2,x3;
	x1=x[0]; x2=x[1];x3=x[2];
	f[0]=-7.7667*x1-340.7667*x2-67.3333*x3+u;
	f[1]=x1;
	f[2]=x2;
	}


	void Modelola(double x[],int k,double h)
	{
		Model(f,x);
		for(int i=0;i<k;i++)
		{ x[i]+=h*f[i];}
		Model(f,x);
	}
	void ModelRgkt_4(double x[], int k, double h)
	{
		double a[4], *b, *f;

		b = new double[k];
		f = new double[k];
		a[0] = h/2.0; a[1] = a[0]; a[2] = h; a[3] = h;
		Model(f, x);
		for (int i = 0; i <=k-1; i++)  b[i] = x[i];
		for (int j = 0; j <=2; j++)
		{
			for (i = 0; i <=k-1; i++)
			{
				x[i] += a[j]*f[i];
				b[i] += a[j+1]*f[i]/3.0;
			}

			Model(f, x);
		}
		for (i = 0; i <=k-1; i++) x[i] = b[i]+h*f[i]/6.0;//解出了上面3个微分方程组
		delete(b);  delete(f);   //
	}
	void ModelRgkt_2(double x[],int k,double h)
	{  double *b;
	b=new double[k];
	Model(f,x);
	for (int i = 0; i <=k-1; i++)b[i] = x[i];
	for(i=0;i<=k-1;i++)
	{

		x[i]+=h*f[i];
	}
	Model(f,x);
	delete(b);
	}