nonlinear.cpp

迭代计算一个四阶黎卡提方程

/***********************************************************/
/*            This is the simulation program CSS02.c               */
/*   Remark: z[i] represents a nonlinear block.                       */
/*           If z[i] is odd, a nonlinear block is located before the block No.i.              */
/*           If z[i] is even, a nonlinear block is located behind the block No.i.             */
/*         If z[i] is equal zero, there isn't a nonlinear block              */
/*                before or behind the block No.i.                 */
/**********************************************************/
#include "stdio.h"
#include "math.h"
#include "conio.h"
#define lpn stdout
#define Kp 23
#define Tp 0.27
/*-------------------- The description of variables ------------------------ */
typedef double REAL[31];
typedef int   INT[31];
REAL a,b,c,d,e,f,g,x,y;
REAL w,k,s,o,q,u,v;
REAL yout[300];
double t1[300],p[30][31];
double t,step,rm;
INT h,z;
int n0,n1,n2,n3,n4,n;
int l1,l2;
/*-------------------------- The description of functions -------------------*/
void subinput_s();
void subinput_b();
void subinput();
void subcoe();
void subhead();
void subdeu();
void subaux();
void subcrt();
void subout();
void sub1(int m , double c1);
void sub2(int m , double c1);
void sub3(int m , double c1);
void sub4(int m );
void sub5(int m );
/*----------------------------- Main function -----------------------------*/
void main()
{
	int iq1,iq2;
	int i,j;
	subinput();
	subcoe();
	t=0.0;
	subhead();
	t1[0]=0.0;
	for(i=1;i<=n;i++)
		yout[0][i]=y[i];
	y[0]=rm;
	for(iq1=1;iq1<=l2;iq1++)
	{ 
		for(iq2=1;iq2<=l1;iq2++)
		{
			for(i=1;i<=n;i++)
  			{
				u[i]=0.0;
				for(j=0;j<=n;j++)
  					u[i]+=p[i][j]*y[j];
			}
			subdeu();
			for(i=1;i<=n;i++)
			{
				v[i]=(u[i]-w[i])/step;
				w[i]=u[i];
				x[i]=e[i]*x[i]+f[i]*u[i]+g[i]*v[i];
				switch(h[i])
				{  
	   	 			case 0:
					case 2:
					case 4:
						y[i]=x[i];
						break;
					case 1:
						y[i]=x[i]+b[i]*k[i]*u[i]+b[i]*k[i]*v[i]*step;
						break;
					case 3:
						y[i]=(b[i]-a[i])*x[i]+k[i]*u[i]+k[i]*v[i]*step;
						break;
					default:
						printf("** h[%d]=%d incorrect!\n",i,h[i]);
				}
			}
			t+=step;
			subaux();
		}
	subcrt();
	for(i=1;i<=n;i++)
    	yout[iq1][i]=y[i];
	t1[iq1]=t;
	}
	subout();
	getch();
}

/*-------------------------------- Function: subinput_s ----------------------*/
void subinput_s()
{ 
  printf("The order of system N=");
  scanf("%d",&n);
  printf("Input N0,N2,N3,N4=");
  scanf("%d,%d,%d,%d",&n0,&n2,&n3,&n4);
  printf("Enter the input R=");
  scanf("%f",&rm);
  printf("Enter the STEP=");
  scanf("%f",&step);
  printf("Enter L1,L2=");
  scanf("%d,%d",&l1,&l2);
}
/*------------------------------- Function: subinput_b ----------------------*/
void subinput_b()
{
  int i,j;
  int i1,j1;
  double p1;
  printf("Enter:\n     A, B, C, D, X, Y, Z, S \n");
  for(i=1;i<=n;i++)
  { 
	 printf("No.%d: ",i);
     scanf("%f,%f,%f,%f,%f,%f,%d,%f",&a[i],&b[i],&c[i],&d[i],&x[i],&y[i], &z[i],&s[i]);
     for(j=0;j<=n;j++) 
		p[i][j]=0.0;
  }
  printf("Enter link matrix:\nI, J , P  \n");
  while(1)
  {
      scanf("%d,%d,%f",&i1,&j1,&p1);
      if((i1>0&&i1<=n)&&(j1>=0&&j1<=n)) 	  
		  p[i1][j1]=p1;
      else  if(i1>n||j1>n)  
		  printf("p[%d][%d] : i or j over range! retry.\n",i1,j1);
      else  break;
   }
}
/*----------------------------- Function: subcoe-----------------------------*/
void subcoe()
{
	int i;
	for(i=1;i<=n;i++)
   	{ 
		if(a[i]==0.0)
        {
          if(d[i]==0.0)
          { 
			h[i]=0;
			k[i]=c[i]/b[i];
			a[i]=0.0;
			b[i]=0.0;
          }
		  else
		  {
			h[i]=1;
			k[i]=c[i]/b[i];
			a[i]=0.0;
			b[i]=d[i]/c[i];
		  }
		}
     else
     {
         if(d[i]==0.0)
         {
             if(b[i]==0.0)
	         {
				h[i]=4;
				k[i]=c[i]/a[i];
				a[i]=0.0;
				b[i]=0.0;
             }
             else
	         {
				h[i]=2;
				k[i]=c[i]/b[i];
				a[i]=a[i]/b[i];
				b[i]=0.0;
	         }
         }
         else
         {
			h[i]=3;
			k[i]=d[i]/b[i];
			a[i]=a[i]/b[i];
			b[i]=c[i]/d[i];
         }
     }
     w[i]=0.0;
	 if(z[i]==6)    
		 q[i]=y[i];
	 else          
		 q[i]=0.0;
	 o[i]=0.0;
	 switch(h[i])
      { 
		case 0:
		case 1:
			e[i]=1.0;
			f[i]=k[i]*step;
			g[i]=k[i]*step*step/2.0;
			break;
		case 2:
		case 3:
			e[i]=exp(-a[i]*step);
			f[i]=k[i]*(1.0-e[i])/a[i];
			g[i]=k[i]*(e[i]-1.0)/(a[i]*a[i])+k[i]*step/a[i];
			break;
		case 4:
			e[i]=0.0;
			f[i]=k[i];
			g[i]=0.0;
      }

	 printf("i=%d\t a=%f\t,b=%f\t\n",i,a[i],b[i]);
   }
}
/*----------------------------- Function: subhead----------------------------*/
void subhead()
{
	switch(n0)
	{ 
		case 1:
			printf("  T                Y(%d)  \n",n2);
			printf("  %15.5f    %15.5f    \n",t,y[n2]);
			break;
		case 2:
			printf("  T                Y(%d)            Y(%d)    \n",n2,n3);
			printf("  %15.5f    %15.5f    %15.5f    \n",t,y[n2],y[n3]);
			break;
		case 3:
			printf("  T                Y(%d)            Y(%d)            Y(%d)  \n",n2,n3,n4);
			printf("  %15.5f    %15.5f    %15.5f    %15.5f   \n",t,y[n2],y[n3],y[n4]);
			break;
		default:
			printf(" cannot output more than 3 varibles.\n");
   }
}
/*----------------------------- Function: subcrt-----------------------------*/
void subcrt()
{
	switch(n0)
	{ 
		case 1:
			printf("  %15.5f    %15.5f    \n",t,y[n2]);
			break;
		case 2:
			printf("  %15.5f    %15.5f    %15.5f    \n",t,y[n2],y[n3]);
			break;
		case 3:
			printf("  %15.5f    %15.5f    %15.5f    %15.5f   \n",t,y[n2],y[n3],y[n4]);
			break;
		default:
			printf(" cannot output more than 3 varibles.\n");
   }
}
/*----------------------------- Function: subdeu-----------------------------*/
void subdeu()
{ 
	int i;
	for(i=1;i<=n;i++)
	{
		switch(z[i])
		{
			case 1: sub1(i,s[i]); break;
			case 3: sub2(i,s[i]); break;
			case 5: sub3(i,s[i]); break;
			case 7: sub4(i)		; break;
			case 9: sub5(i)		; break;
			case 0:
			case 2:
			case 4:
			case 6:
			case 8:
			case 10:break;
			default: printf("z[%d] value incorrect.\n",i);
		}
   }
}
/*--------------------------- Function: subaux-----------------------------*/
void subaux()
{
	int i;
	for(i=1;i<=n;i++)
	{
		u[i]=y[i];
		switch(z[i])
		{
			case 2: sub1(i,s[i]);y[i]=u[i]; break;
			case 4: sub2(i,s[i]);y[i]=u[i]; break;
			case 6: sub3(i,s[i]);y[i]=u[i]; break;
			case 8:	sub4(i);	 y[i]=u[i]; break;
			case 10:sub5(i);	 y[i]=u[i]; break;
			case 0:
			case 1:
			case 3:
			case 5:
			case 7:
			case 9: break;
			default: printf("z[%d] value incorrect.\n",i);
		}
	}
}
/*---------------------------- Function: sub1-----------------------------*/
void sub1(int m,double c1)
{
	if(fabs(u[m])>=c1)
	{
		if(u[m]>=c1)
    		 u[m]=c1;
		else
			 u[m]=-c1;
	}
	else
   		u[m]=u[m];
}
void sub2(int m,double c1)
{
	if(fabs(u[m])>=c1)
	{
		if(u[m]>=c1)
			u[m]-=c1;
        else
			u[m]+=c1;
	}
	else
  		u[m]=0.0;
}
void sub3(int m,double c1)
{
	if(u[m]<=o[m])
	{
		if(u[m]<o[m])
		{
			if(q[m]>=(u[m]+c1))
     		{ 
				o[m]=u[m];
				u[m]+=c1;
				q[m]=u[m];
			}
			else
			{ 
				o[m]=u[m];
				u[m]=q[m];
				q[m]=u[m];
			}
		}
		else
		{
			o[m]=u[m];
			u[m]=q[m];
			q[m]=u[m];
		}
	}
	else
	{
		o[m]=u[m];
		if(q[m]<=(u[m]-c1))
		{
			u[m]-=c1;
			q[m]=u[m];
		}
		else
		{
			u[m]=q[m];
			q[m]=u[m];
		}
	}
}

/* ---------------------Function: subout -------------------------------*/
void subout()
{
	FILE *file=fopen("data.txt","w");
	int i;
	switch(n0)
    {
       	  case 1:
			printf("	  T                Y(%d)  \n",n2);
			break;
		  case 2:
			printf("	  T                Y(%d)             Y(%d)    \n",n2,n3);
			break;
          case 3:
			printf("	  T                Y(%d)             Y(%d)             Y(%d)  \n",n2,n3,n4);
			break;
          default:
			printf(" cannot output more than 3 varibles.\n");
	}
	for(i=0;i<=l2;i++)
	{  
		switch(n0)
		{
			case 1:
				printf("  %15.5f    %15.5f    \n",t1[i],yout[i][n2]);
				break;
			case 2:
				printf("  %15.5f    %15.5f    %15.5f    \n",t1[i],yout[i][n2] ,yout[i][n3]);  
				break;
			case 3:
				printf("  %15.5f    %15.5f    %15.5f    %15.5f   \n",t1[i],yout[i][n2],yout[i][n3],yout[i][n4]);
				fprintf(file,"%f\t%f\n",t1[i],yout[i][n4]);
				break;
            default:
	           printf(" cannot output more than 3 varibles.\n");
        }
    }
	fclose(file);
 }
/* -------------------- The end of program CSS02.c-----------------------*/

/* --------------------"以下是修改的输入部分子程序"----------------------*/
void sub4(int m)
{
	if(fabs(u[m]>=6.01))
		u[m]=0;
	else if(fabs(u[m]<=3.01))
		u[m]=1.46;
	else
	{
		if(u[m]>0)
			u[m]=1.46-(u[m]-3.01)*3*1.46;
		else
			u[m]=1.46-(fabs(u[m])-3.01)*3*1.46;
	}
}

void sub5(int m)
{
	double temp;
	temp=fabs(u[m]);
	if(temp>=113.7)
		temp=(temp-113.7)*0.00046+0.2145;
	else if(temp>=27.97)
		temp=(temp-27.97)*0.0009+0.1375;
	else if(temp>=8.53)
		temp=(temp-8.53)*0.0025+0.088;
	else 
		temp=temp*0.013;
	if(u[m]>=0)
		u[m]=temp;
	else
		u[m]=-1*temp;
}

void subinput()
{
	int i,j;
	
	n=24;
	n0=3;n2=1;n3=18;n4=19;
	rm=1.0;
	step=0.001;
	l1=20;l2=200;

	for(i=1;i<=n;i++)
	{ 
		x[i]=0.0;y[i]=0.0;
		z[i]=0;s[i]=0.0;
		for(j=0;j<=n;j++)
		p[i][j]=0.0;
	}
	
	i=1;a[i]=0.0;b[i]=1.0;c[i]=Kp;d[i]=Kp*Tp;
	i=2;a[i]=27.9;b[i]=0.0;c[i]=100.0;d[i]=0.0;z[i]=2;s[i]=10.0;
	i=3;a[i]=1.0;b[i]=0.0;c[i]=1.0;d[i]=0.0;z[i]=1;s[i]=10.0;
	i=4;a[i]=1.0;b[i]=0.0;c[i]=1.0;d[i]=0.0;z[i]=1;s[i]=10.0;
	i=5;a[i]=1.0;b[i]=0.0033;c[i]=27.9;d[i]=0.0;z[i]=3;s[i]=1.2;
	i=6;a[i]=1.0;b[i]=0.0033;c[i]=27.9;d[i]=0.0;z[i]=3;s[i]=1.2;
	i=7;a[i]=2.0;b[i]=0.25;c[i]=1.0;d[i]=0.0;
	i=8;a[i]=2.0;b[i]=0.25;c[i]=1.0;d[i]=0.0;
	i=9;a[i]=1.0;b[i]=0.0;c[i]=0.055;d[i]=0.0;
	i=10;a[i]=1.0;b[i]=0.0;c[i]=0.055;d[i]=0.0;
	i=11;a[i]=0.0;b[i]=0.000566;c[i]=1.0;d[i]=0.0;
	i=12;a[i]=0.0;b[i]=0.000566;c[i]=1.0;d[i]=0.0;
	i=13;a[i]=1.0;b[i]=0.0;c[i]=1.0;d[i]=0.0;z[i]=10;s[i]=0;
	i=14;a[i]=1.0;b[i]=0.0;c[i]=1.0;d[i]=0.0;z[i]=10;s[i]=0;
	i=15;a[i]=0.0;b[i]=1.0;c[i]=0.517;d[i]=0.0;
	i=16;a[i]=0.0;b[i]=1.0;c[i]=0.517;d[i]=0.0;
	i=17;a[i]=0.0011856;b[i]=0.0;c[i]=1.0;d[i]=0.0;
	i=18;a[i]=0.0;b[i]=1.0;c[i]=1.0;d[i]=0.0;
	i=19;a[i]=0.0;b[i]=1.0;c[i]=1.0;d[i]=0.0;
	i=20;a[i]=1.0;b[i]=0.0034;c[i]=1.0;d[i]=0.0;
	i=21;a[i]=1.0;b[i]=0.0034;c[i]=1.0;d[i]=0.0;
	i=22;a[i]=1.0;b[i]=0.0034;c[i]=1.0;d[i]=0.0;
	i=23;a[i]=1.0;b[i]=0.0;c[i]=1.0;d[i]=0.0;z[i]=3;s[i]=0.2517;
	i=24;a[i]=1.0;b[i]=0.0;c[i]=1.0;d[i]=0.0;z[i]=3;s[i]=0.2517;
	i=25;a[i]=1.0;b[i]=0.0;c[i]=1.0;d[i]=0.0;z[i]=7;s[i]=0;

	 p[1][0]=1.0;p[1][19]=-1.0;
	 p[2][1]=1.0;p[2][7]=-0.2512201;p[2][8]=-0.2512201;
	 p[2][17]=-0.0092166;p[2][20]=-0.114088;
	 p[2][21]=-0.1505775;p[2][22]=-0.1505775;

	 p[3][2]=1.0;p[3][25]=1.0;  
	 p[4][2]=1.0;p[4][25]=1.0;
	 p[5][3]=1.0;  p[6][4]=1.0;
	 p[7][5]=1.0; p[7][11]=-0.544;
	 p[8][6]=1.0; p[8][12]=-0.544;
	 p[9][7]=1.0; p[10][8]=1.0;
	 p[11][9]=1.0; p[11][13]=-1.0;p[11][23]=-1.0;
	 p[12][10]=1.0; p[12][14]=-1.0;p[12][24]=-1.0;

	 p[13][11]=1.0;  p[14][12]=1.0;
	 p[15][11]=1.0;  p[16][12]=1.0;
	 p[17][0]=-0.2; p[17][23]=1.0; p[17][24]=1.0;
	 p[18][17]=1.0;
	 p[19][18]=1.0;

	 p[20][18]=1.0;
	 p[21][11]=1.0;p[22][12]=1.0;
	 p[23][15]=1.0;p[23][19]=-0.517;
	 p[24][16]=1.0;p[24][19]=-0.517;

	 p[25][7]=0.2512201;p[25][8]=0.2512201;
}