example2_16.java

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

import java.applet.*;
import java.awt.*;
public  class Example2_16 extends Applet
{  public static  void MHessenberg(double[][] a,int n)
   {int i=0,j,k; double d,t;
    for(k=0;k<=n-3;k++)
    {d=0;
      for(j=k+1;j<n;j++)
      {t=a[j][k];
       if (Math.abs(t)>Math.abs(d)){d=t;i=j;}}

      if(d!=0)
      {if(i!=(k+1))
       {for(j=k;j<n;j++)
        {t=a[i][j];a[i][j]=a[k+1][j];a[k+1][j]=t;}
        for(j=0;j<n;j++)
        {t=a[j][i];a[j][i]=a[j][k+1];a[j][k+1]=t;}
       }
	   for(i=k+2;i<n;i++)
	   {t=a[i][k]/d;a[i][k]=0;
	    for(j=k+1;j<n;j++)a[i][j]-=t*a[k+1][j];
	    for(j=0;j<n;j++)a[j][k+1]+=t*a[j][i];}
	  }
    }
  }


public static  void francis(int n,double[][]a ,double[] wr,double[] wi)
 {int nn=n-1,i,j,k,loop=0;
  double s,t,x,y,z,alpha,beta,temp;
  boolean falg;
  while(nn>=1&&loop<100)
  {loop++;
   s=a[nn][nn]+a[nn-1][nn-1];
   t=a[nn-1][nn-1]*a[nn][nn]-a[nn][nn-1]*a[nn-1][nn];
   x=a[0][0]*a[0][0]+a[0][1]*a[1][0]-s*a[0][0]+t;
   y=a[1][0]*(a[0][0]+a[1][1]-s);z=a[1][0]*a[2][1];
   for(k=0;k<=nn+1-3;k++)
   {alpha=y*y+z*z;
    if(alpha!=0)
    {alpha=Math.sqrt(x*x+alpha);
     if(x<0)alpha=-alpha;
     beta=alpha*(alpha+x);
     for(j=0;j<=nn;j++)
    {temp=(x+alpha)*a[k][j]+y*a[k+1][j]+z*a[k+2][j];
     a[k][j]-=temp/alpha;a[k+1][j]-=y*temp/beta;a[k+2][j]-=z*temp/beta;}
    for(i=0;i<=nn;i++)
    {t=(x+alpha)*a[i][k]+y*a[i][k+1]+z*a[i][k+2];
     a[i][k]-=t/alpha;a[i][k+1]-=y*t/beta;a[i][k+2]-=z*t/beta;}}
    x=a[k+1][k];y=a[k+2][k];
    if(k<nn+1-3) z=a[k+3][k];
    }
    if(y!=0)
    {alpha=Math.sqrt(x*x+y*y);if(x<0)alpha=-alpha;
     beta=alpha*(alpha+x);
     for(j=0;j<=nn;j++)
      {temp=(x+alpha)*a[nn-1][j]+y*a[nn][j];
       a[nn-1][j]-=temp/alpha;a[nn][j]-=y*temp/beta;}
     for(i=0;i<=nn;i++)
     {temp=(x+alpha)*a[i][nn-1]+y*a[i][nn];
     a[i][nn-1]-=temp/alpha;a[i][nn]-=y*temp/beta;}}
     //降阶处理
     falg=true;
     while(falg)
     {if(nn>=1&&(Math.abs(a[nn][nn-1])+a[nn][nn]==a[nn][nn])){wr[nn]=a[nn][nn];wi[nn]=0;nn--;}
      else if(nn>=2&&(Math.abs(a[nn-1][nn-2])+a[nn-1][nn-1]==a[nn-1][nn-1]))
      {s=a[nn][nn]+a[nn-1][nn-1];t=a[nn][nn]*a[nn-1][nn-1]-a[nn][nn-1]*a[nn-1][nn];
       wr[nn]=wr[nn-1]=s/2;wi[nn]=Math.sqrt(t-s*s/4);wi[nn-1]=-wi[nn];nn-=2;}
      else falg=false;}
     if(nn==0){wr[0]=a[0][0];wi[0]=0;break;}
     else if(nn==1){s=a[0][0]+a[1][1];t=a[0][0]*a[1][1]-a[1][0]*a[0][1];
      if(s*s>=4*t){wr[0]=s/2+Math.sqrt(s*s/4-t);wr[1]=s/2-Math.sqrt(s*s/4-t);wi[0]=wi[1]=0;}
      else {wr[0]=wr[1]=s/2;wi[0]=Math.sqrt(t-s*s/4);wi[1]=-wi[0];break;}}
    }
  }

  public void paint(Graphics g)
  {double[][] a={{1,6,-3,-1,7},{8,-15,18,5,4},{-2,11,9,15,20},{-13,2,21,30,-6},{17,22,-5,3,6}};
   double[] wr=new double[5];
   double[] wi=new double[5];
   int i;
   MHessenberg(a,5);
   francis(5,a,wr,wi);
   for(i=0;i<5;i++)
   {if(wi[i]==0)g.drawString("第"+(i+1)+"特征值为"+wr[i],10,20+20*i);
    else g.drawString("第"+(i+1)+"特征值为("+wr[i]+")+("+wi[i]+")I",10,20+20*i);}
  }
}