chengxu.txt

用vc++求解五阶矩阵的算法代码

#include "stdio.h"
#include "math.h"

//输出函数
void print(double *a,int n)
{
	int i,j;
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			if(j%n==0)printf("\n");
			printf("%10.4f",a[i*n+j]);
		}
	printf("\n");
}

//gauss消去法求解线性方程组 ax=b
int agaus(double a[],double b[],int n)
  { int *js,l,k,i,j,is,p,q;
    double d,t;
    //js=malloc(n*sizeof(int));
	js=new int[n];
    l=1;
    for (k=0;k<=n-2;k++)
      { d=0.0;
        for (i=k;i<=n-1;i++)
          for (j=k;j<=n-1;j++)
            { t=fabs(a[i*n+j]);
              if (t>d) { d=t; js[k]=j; is=i;}
            }
        if (d+1.0==1.0) l=0;
        else
          { if (js[k]!=k)
              for (i=0;i<=n-1;i++)
                { p=i*n+k; q=i*n+js[k];
                  t=a[p]; a[p]=a[q]; a[q]=t;
                }
            if (is!=k)
              { for (j=k;j<=n-1;j++)
                  { p=k*n+j; q=is*n+j;
                    t=a[p]; a[p]=a[q]; a[q]=t;
                  }
                t=b[k]; b[k]=b[is]; b[is]=t;
              }
          }
        if (l==0)
          { //free(js);
			delete js;
			printf("fail\n");
            return(0);
          }
        d=a[k*n+k];
        for (j=k+1;j<=n-1;j++)
          { p=k*n+j; a[p]=a[p]/d;}
        b[k]=b[k]/d;
        for (i=k+1;i<=n-1;i++)
          { for (j=k+1;j<=n-1;j++)
              { p=i*n+j;
                a[p]=a[p]-a[i*n+k]*a[k*n+j];
              }
            b[i]=b[i]-a[i*n+k]*b[k];
          }
      }
    d=a[(n-1)*n+n-1];
    if (fabs(d)+1.0==1.0)
      { //free(js);
		delete js;
		printf("fail\n");
        return(0);
      }
    b[n-1]=b[n-1]/d;
    for (i=n-2;i>=0;i--)
      { t=0.0;
        for (j=i+1;j<=n-1;j++)
          t=t+a[i*n+j]*b[j];
        b[i]=b[i]-t;
      }
    js[n-1]=n-1;
    for (k=n-1;k>=0;k--)
      if (js[k]!=k)
        { t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;}
    //free(js);
	  delete js;
    return(1);
  }
/*****************************************************/

 void main()
  { 
    //static double a[4][4]=    
    //       { {0.2368,0.2471,0.2568,1.2671},
    //         {0.1968,0.2071,1.2168,0.2271},
    //         {0.1581,1.1675,0.1768,0.1871},
    //         {1.1161,0.1254,0.1397,0.1490} };
    //static double b[4]={1.8471,1.7471,1.6471,1.5471};

	double a[5][5],b[5];
	/*  double a[5][5]={{0.9501,0.7621,0.6154,0.4057,0.0579},
				{0.2311,0.4565,0.7919,0.9355,0.3529},
				{0.6068,0.0185,0.9218,0.9169,0.8132},
				{0.4860,0.8214,0.7382,0.4103,0.0099},
				{0.8913,0.4447,0.1763,0.8936,0.1389}};
	  double b[]={0.2028,0.1987,0.6038,0.2722,0.1988};*/
	double x0,y0,theta,aa,bb,A,B,C,D,E,t1,t2;
	int i,j,k,n=9;
	double x[]={0,0.5,1,1.5,2,2.5,3,3.5,4};
	double y[]={5,1,1.5,2,2.5,3,3.5,4,5};


	//printf("\ninput x[]=\n");
	//for(k=0;k<n;k++)scanf("%f",&x[k]);
	//printf("\ninput y[]=\n");
	//for(k=0;k<n;k++)scanf("%f",&y[k]);
	
	for(i=0;i<5;i++)
	{   b[i]=0;
		for(j=0;j<5;j++)
			a[i][j]=0;
	}
	print(*a,5);
	for(k=0;k<n;k++)
	{
		a[0][0]=a[0][0]+x[k]*x[k]*y[k]*y[k];
		a[0][1]=a[0][1]+x[k]*y[k]*y[k]*y[k];
		a[0][2]=a[0][2]+x[k]*x[k]*y[k];
		a[0][3]=a[0][3]+x[k]*y[k]*y[k];
		a[0][4]=a[0][4]+x[k]*y[k];

		a[1][1]=a[1][1]+y[k]*y[k]*y[k]*y[k];
		a[1][2]=a[1][2]+x[k]*y[k]*y[k];
		a[1][3]=a[1][3]+y[k]*y[k]*y[k];
		a[1][4]=a[1][4]+y[k]*y[k];

		a[2][2]=a[2][2]+x[k]*x[k];
		a[2][3]=a[2][3]+x[k]*y[k];
		a[2][4]=a[2][4]+x[k];

		a[3][3]=a[3][3]+y[k]*y[k];
		a[3][4]=a[3][4]+y[k];
		a[4][4]=n;

		b[0]=b[0]+x[k]*x[k]*x[k]*y[k];
		b[1]=b[1]+x[k]*x[k]*y[k]*y[k];
		b[2]=b[2]+x[k]*x[k]*x[k];
		b[3]=b[3]+x[k]*x[k]*y[k];
		b[4]=b[4]+x[k]*x[k];
	}
	
	for(i=1;i<5;i++)
		for(j=0;j<i;j++)
			a[i][j]=a[j][i];
	print(*a,5);
	for(k=0;k<5;k++)
		printf("%10.4f",b[k]);
	printf("\n");
	


    if (agaus(*a,b,5)!=0)
      for (i=0;i<5;i++)
        printf("x(%d)=%10.4f\n",i,b[i]);
	print(*a,5);
	A=b[0];B=b[1];C=b[2];D=b[3];E=b[4];
	x0=(2*B*C-A*D)/(A*A-4*B);
	y0=D*(2-A)/(A*A-4*B);
	t1=2*(A*C*D-B*C*C-D*D+4*B*E-A*A*E);
	t2=(A*A-4*B)*(B+1-sqrt(A*A+(1-B)*(1-B)));
	aa=sqrt(t1/t2);
	t2=(A*A-4*B)*(B+1+sqrt(A*A+(1-B)*(1-B)));
	bb=sqrt(t1/t2);
	theta=atan(sqrt((aa*aa-bb*bb*B)/(aa*aa*B-bb*bb)));
  }
/*****************************************************/