example4_5.java

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

import java.applet.*;
import java.awt.*;
public class Example4_5 extends Applet
{
 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 static void QRroot(int n,double[] b,double[] wr,double[] wi )
  { double[][] a=new double[n][n];
    int i;
    for(i=1;i<n;i++)a[i][i-1]=1;
    for(i=0;i<n;i++)a[0][i]=-b[n-1-i]/b[n];
    francis(n,a,wr,wi);
  }

  public void paint(Graphics g)
  {double[] b={-60,21,-21,3,9,-15,3};
   double[] r1=new double[6];
   double[] r2=new double[6];
   QRroot(6,b,r1,r2);
   int i;
   g.drawString("方程的6个根为",0,20);
   for(i=0;i<6;i++)
   {if(r2[i]==0)g.drawString(""+r1[i],10,40+20*i);
    else if(r2[i]<0)g.drawString(""+r1[i]+""+r2[i]+"  I",10,40+20*i);
    else g.drawString(""+r1[i]+"+"+r2[i]+"  I",10,40+20*i);
   }
  }
  }