danbai01.c

本程序以自然界最平常的单摆现象为案例

#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include <graphics.h>
#include<dos.h>
#include<conio.h>
void inisub(void);
void runsub(void);
void intsub(void);
//void stosub(void);
void outsub(void);
void modsub(void);
void repsub(void);
void blot(void);
float t,t1,t2,t3,t4,t7,l,s;
int l3,j8,j9,n1,n2,n3,n4,n5;
float y[31],m[31],e[31],g[31],p[21],z[21],k[21][5],w[802][6],a[5];
int h,j,q,gdrive=VGA,gmode=VGAHI;
void main()
{
	initgraph(&gdrive,&gmode,".\\bgi");
	delay(2000);
	inisub();
	t1=50;t2=0.01;n1=3;n2=0;n3=5000;y[1]=0.6;y[2]=1;y[3]=0;
	l=y[2]/y[1];
	cleardevice();
	runsub();
	closegraph();
       }
      void runsub(void)
      {
	int i;
	for(i=1;i<=n1;i++)
	  z[i]=y[i];
	  n4=1;
	  t=0.0;
	  t4=t1/(n3-1);
	  if(t2>=t4) t2=t4;
	  t3=t4;
	  t7=0.5*t2;
	modsub();
       do
       {
	 do intsub();while(t<t3);
	 n4++;
	 t3=n4*t4;
	 s=kbhit();
	}while(!s);
	for(i=1;i<=n1;i++)
	   y[i]=z[i];
	}
	void inisub(void)
	{
	   int i,j;
	   for(i=0;i<=20;i++)
	   {
	       p[i]=0.0;
	       g[i]=0.0;
	       z[i]=0.0;
	       for(j=0;j<=4;j++)
	       {
		      k[i][j]=0.0;
		      a[j]=0.0;
		}
	    }
		t=0.0;
		l3=0.0;
		for(i=0;i<=30;i++)
		{
		   y[i]=0.0;
		   m[i]=0.0;
		   e[i]=0.0;
		 }
		 j9=0;
		 for(i=1;i<=801;i++)
		 {
		    for(j=0;j<=5;j++)
			w[i][j]=0.0;
		  }
	 }
	 void intsub(void)

	 {
	     int i;
	     for(i=1;i<=n1;i++)
	     {
		k[i][2]=g[i]*t7;
		k[i][1]=y[i];
		y[i]=k[i][1]+k[i][2];
	      }
	      t=t+t7;
	      modsub();
	      for(i=1;i<=n1;i++)
	      {
		  k[i][3]=g[i]*t7;
		  y[i]=k[i][1]+k[i][3];
	      }
	       modsub();
	       for(i=1;i<=n1;i++)
	       {
		  k[i][4]=g[i]*t2;
		  y[i]=k[i][1]+k[i][4];
		}
		  t=t+t7;
		  modsub();
		for(i=1;i<=n1;i++)
		     y[i]=k[i][1]+(2*k[i][2]+2*k[i][4]+4*k[i][3]+g[i]*t2)/6;
			blot();
		modsub();
	  }
	  void stosub(void)
	  {
	      int i;
	      n5=n4-1;
	      w[n5][1]=t;
	      repsub();
	      for(i=1;i<=30;i++)
	      {
		if(y[i]>m[i]) m[i]=y[i];
		if(y[i]<e[i]) e[i]=y[i];
	       }
	       for(i=2;i<=5;i++)
		    w[n5][i]=a[i-1];
	  }
	  void outsub(void)
	  {
	     int i,j;
	     printf("Y[I]                 Y[I]max             Y[I]min\n");
	     for(i=1;i<=30;i++)
		if(m[i]!=0||e[i]!=0)
		     printf("%2d          %15.5f              %15.5f\n",i,m[i],e[i]);
	     printf("SIMULATION STEP T2=%15.5f\n",t2);
	     printf("    T     A[1]       A[2]       A[3]        A[4]\n");
	     for(i=0;i<=n5;i++)
	     {
		for(j=1;j<=5;j++)
		    printf("%13.5f",w[i][j]);
		printf("\n");
	      }
	   }
	   void modsub(void)
	   {
		g[1]=y[3];
		g[2]=l*y[3];
		g[3]=-9.8/l*sin(y[1]);
		y[4]=y[2]/y[1];
	    }
	    void repsub(void)
	   {
	       a[1]=y[1];
	       a[2]=y[2];
	       a[3]=y[3];
	       a[4]=y[4];
	   }

	void blot()
	{
	    double j=1;q=0;
	    j=(int)(100*l*sin(y[1]));
	    q=(int)(100*l*cos(y[1]));
	    
	    setcolor(WHITE);
	    moveto(300,100);
	    lineto(300+j,100+q+1);
	    circle(300+j,100+q,4);
	    delay(10);
	    setcolor(BLACK);
	    moveto(300,100);
	    lineto(300+j,100+q+1);
	    setcolor(BLACK);
	    circle(300+j,100+q,4);
  }