1004ac.cpp

数值分析的c语言几种算法实现

#include <stdio.h>

#include <math.h>



#define MAX_SIZE 100

#define bound pow(2, 127)



int SOR(int n, double a[][MAX_SIZE], double b[], double x[],

       double w, double TOL, int MAXN);



int main()

{

    int n, MAXN, m, i, j, k;

    double a[MAX_SIZE][MAX_SIZE], b[MAX_SIZE], x[MAX_SIZE];

    double w, TOL;



    while (scanf("%d", &n) != EOF) {

       for (i=0; i<n; i++) {

           for (j=0; j<n; j++)

              scanf("%lf", &a[i][j]);

           scanf("%lf", &b[i]);

       }

       scanf("%lf %d", &TOL, &MAXN);

       scanf("%d", &m);

       for (i=0; i<m; i++) {

           scanf("%lf", &w);

           switch (SOR(n, a, b, x, w, TOL, MAXN)) {

              case -2:

                  printf("No convergence.\n");

                  break;

              case -1:

                  printf("Matrix has a zero column.  No unique solution exists.\n");

                  break;

              case 0:

                  printf("Maximum number of iterations exceeded.\n");

                  break;

              case 1:

                  for (k=0; k<n; k++)

                     printf("%12.8f\n", x[k]);

                  break;

           }

       }
       printf("\n");

    }

    return 0;

}


int SOR(int n, double a[][MAX_SIZE], double b[], double x[], double w, double TOL, int MAXN)
{
	
	
   double S, S1;
   int i,j,k;
   double dX;
   double min,max;
   for(i=0;i<n;i++) x[i]=0.0;
   for(i=0;i<n;i++)
   {
      min=fabs(a[i][i]);
      k=i;     
	  for(j=i;j<n;j++)
      {
			if(fabs(a[j][i])>min) 
			{
				k=j;
				min=fabs(a[j][i]);
			 }
      }
      if(min>1e-10)
	  {
		    for(j=0;j<n;j++)
			{ max=a[i][j];a[i][j]=a[k][j];a[k][j]=max;}
			max=b[i];b[i]=b[k];b[k]=max;//b[i]+=b[k];
      }
     else if(min<1e-10)
	 {
		 min=0;
         k=i;
         for(j=i-1;j>=0;j--)
			if(fabs(a[j][i])>min) 
			{
				k=j;
				min=fabs(a[j][i]);
			 }

			if(min>1e-10)
			{
			  for(j=0;j<n;j++)
					a[i][j]+=a[k][j];
			  b[i]+=b[k];
			}
			else 
				return -1;
    }
  }
  /*for(i=0;i<n;i++)
	{  
		for(j=0;j<n;j++)	  
		  printf("%lf ",a[i][j]);
		  printf("%lf\n",b[i]);
    }*/
   for (k=0; k<MAXN; k++)
   {
     S = 0.0;
     for (i=0; i<n; i++)
     {
       double sum = 0.0;       
	   for ( j=0; j<n; j++)
			 sum += a[i][j]*x[j];
       
	   if (fabs(a[i][i]) < 1e-10)
			return -1;

       dX = w * (b[i] - sum ) /a[i][i];
       S1 = dX;
	   // printf("%lf\n",dX[i]);
       if (fabs(S1) > fabs(S)) S = S1;
	   //printf("%lf\n",S);
       x[i] = x[i] + S1;
	  // printf("%lf\n",pow((double)2,(double)127));
       if(fabs(x[i])>pow((double)2,(double)127)) 
	   {
		 //  printf("%lf\n",S);
		 //  printf("%d\n",k);
		  // cout << k << endl;
		   return -2;
	   }
     }
     if (fabs(S)<=TOL)
		 return 1;
   }
   return 0;
}