example3_16.java

清华大学2002年出版的《科学与工程数值计算算法Java》配套源码

import java.applet.*;
import java.awt.*;
public class Example3_16 extends Applet
{
  public static void  Matrixequal(double[][] ta,double[][] a,int n,int m)
  {int i,j;
   for(i=0;i<n;i++)
    for(j=0;j<m;j++)
     ta[i][j]=a[i][j];
   }


 public static double GaussSolve(double[][] a,double[][] b,double[][] c,int n,int m)
  { int i,j,k,row,colunm;
    double temp,max,d=1;
    int[] p=new int[n];
    for(i=0;i<n;i++)p[i]=i;
      for(i=0;i<n-1;i++)
      { //找主元
		row=i;colunm=i;
        max=Math.abs(a[i][i]);
	    for(j=i;j<n;j++)
         for(k=i;k<n;k++)
          {temp=Math.abs(a[j][k]);
           if(temp>max){max=temp;row=j;colunm=k;}
	      }
	    if (max==0){return 0;}
	    //第row行、第i行互换
	    if(row!=i)
	    { for(k=i;k<n;k++)
	      {temp=a[i][k];a[i][k]=a[row][k];a[row][k]=temp;}
	      for(k=0;k<m;k++)
	      {temp=b[i][k];b[i][k]=b[row][k];b[row][k]=temp;}
	    }
	    //第column列、第i列互换
	    if(colunm!=i)
	    { for(j=0;j<n;j++)
	      {temp=a[j][colunm];a[j][colunm]=a[j][i];a[j][i]=temp;}
	      k=p[i];p[i]=p[colunm];p[colunm]=k;
	      }
	    d*=a[i][i];
	    //第i行归一化
	     for(k=i+1;k<n;k++)a[i][k]/=a[i][i];
	     for(k=0;k<m;k++)b[i][k]/=a[i][i];
	     a[i][i]=1;
	    //消元
	    for(j=i+1;j<n;j++)
	     { for(k=i+1;k<n;k++)a[j][k]-=a[j][i]*a[i][k];
		   for(k=0;k<m;k++)b[j][k]-=a[j][i]*b[i][k];
		   a[j][i]=0;
	     }
	   }
       d*=a[n-1][n-1];
       //回代
       for(k=0;k<m;k++)
       {b[n-1][k]/=a[n-1][n-1];
        for(i=n-2;i>=0;i--)
          for(j=i+1;j<n;j++)
           b[i][k]-=a[i][j]*b[j][k];
       }

       //调整次序
       for(i=0;i<n;i++)
         for(j=0;j<m;j++)
         c[p[i]][j]=b[i][j];

       return d;
  }


public static void Morbid(double[][] a,double[] b,int n,double eps)
 {double[][] ta=new double[n][n];
  double[][] tb=new double[n][1];
  double[][] x=new double[n][1];
  double[][] y=new double[n][1];
  double[][] r=new double[n][1];
  int i,j;
  Matrixequal(ta,a,n,n);
  for(i=0;i<n;i++)tb[i][0]=b[i];
  GaussSolve(ta,tb,x,n,1);
  while(true)
  { for(i=0;i<n;i++)
    {r[i][0]=b[i];
     for(j=0;j<n;j++)r[i][0]-=a[i][j]*x[j][0];}
    Matrixequal(ta,a,n,n);
    GaussSolve(ta,r,y,n,1);
    for(i=0;i<n;i++)
    {if(Math.abs(y[i][0])/(1+Math.abs(x[i][0]+y[i][0]))>=eps)break;}
    if(i==n)break;
    for(i=0;i<n;i++)x[i][0]+=y[i][0];
    Matrixequal(ta,a,n,n);
}
for(i=0;i<n;i++)b[i]=x[i][0]+y[i][0];
}


 public void paint(Graphics g)
 {double[][] a=new double[5][5];
   double[] b={1,0,1,0,1};
   double[] c={1.001,-0.001,0.999,0.001,1.001};
   int i,j;
   double eps=1e-7;
   for(i=0;i<5;i++)
     for(j=0;j<5;j++)
      a[i][j]=1.0/(i+j+1);
 Morbid(a,b,5,eps);
 Morbid(a,c,5,eps);
 for(i=0;i<5;i++)g.drawString("X1"+i+"="+b[i]+"  X2"+i+"="+c[i],10,i*20+20);
}
}