example2_14.java

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

import java.applet.*;
import java.awt.*;
public class Example2_14 extends Applet
{
 public static void  MSymTri(double[][] a,double[][] Q,int n)
  {int i,j,k,l;
   double[] u=new double[n];
   double[] q=new double[n];
   double alpha,beta,temp;
   for(i=0;i<n;i++)
   {Q[i][i]=1;for(j=i+1;j<n;j++){Q[j][i]=0;Q[i][j]=0;}}

  for(k=0;k<n-2;k++)
  { l=n-k-1;
	alpha=0;
    for(j=0;j<l;j++)alpha+=a[l][j]*a[l][j];
    if(alpha==0)continue;
    alpha=Math.sqrt(alpha);
    if(a[l][l-1]<0)alpha*=-1;
    beta=(alpha+a[l][l-1])*alpha;
    for(i=0;i<l-1;i++)u[i]=a[l][i];
    u[l-1]=a[l][l-1]+alpha;
    for(i=l;i<n;i++)u[i]=0;
    for(i=0;i<=l;i++)
    {q[i]=0;
     for(j=0;j<l;j++)q[i]+=a[i][j]*u[j];
     q[i]/=beta;}
     temp=0;
     for(j=0;j<l;j++)temp+=u[j]*q[j];
     temp/=2*beta;
    for(i=0;i<l;i++)q[i]-=temp*u[i];
    for(i=0;i<l;i++)
      for(j=0;j<l;j++)
     {a[i][j]-=q[i]*u[j]+q[j]*u[i];}
    for(i=0;i<l-1;i++)a[l][i]=a[i][l]=0;
    a[l][l-1]=a[l-1][l]=-alpha;

     for(i=0;i<n;i++)
     {q[i]=0;
      for(j=0;j<l;j++)q[i]+=Q[i][j]*u[j];}

    for(i=0;i<n;i++)
      for(j=0;j<l;j++)
      Q[i][j]-=q[i]*u[j]/beta;
   }
}

public static void  QR(double[][] a,double[][] Q,int m,int n )
  {int i,j,k,l=Math.max(m,n);
   double[] t=new double[l];
   double alpha,beta;
   for(i=0;i<m;i++)
    for(j=0;j<m;j++)
     {if (i==j) Q[i][j]=1;
      else Q[i][j]=0; }

  for(k=0;k<Math.min(n,m)-1;k++)
  { alpha=0;
    for(l=k;l<m;l++)alpha+=a[l][k]*a[l][k];
    alpha=Math.sqrt(alpha);
    if(alpha==0)continue;
    if(a[k][k]<0)alpha*=-1;
    beta=(alpha+a[k][k])*alpha;
    for(i=0;i<m;i++)
     {t[i]=Q[i][k]*alpha;
      for(l=k;l<m;l++)t[i]+=Q[i][l]*a[l][k];
     }

    for(i=0;i<m;i++)
    { Q[i][k]-=t[i]/alpha;
	  for(j=k+1;j<m;j++)Q[i][j]-=(t[i]*a[j][k]/beta);}
    for(j=0;j<n;j++)
     {t[j]=alpha*a[k][j];
      for(l=k;l<m;l++)t[j]+=a[l][k]*a[l][j];
     }

   for(j=k+1;j<n;j++)a[k][j]-=t[j]/alpha;
   for(i=k+1;i<m;i++)
     for(j=k+1;j<n;j++)a[i][j]-=(a[i][k]*t[j]/beta);
   for(l=k+1;l<m;l++)a[l][k]=0;
   a[k][k]=-alpha;
   }
 }


 public void Mrcheng(double[][] a,double[][] b,double[][]c,int m,int n,int l)
   {double[][] d=new double[m][l];
    int i,j,k;
    for(i=0;i<m;i++)
      for(j=0;j<l;j++)
       {d[i][j]=0;
        for(k=0;k<n;k++)d[i][j]+=a[i][k]*b[k][j];}
   for(i=0;i<m;i++)
    for(j=0;j<l;j++)
     c[i][j]=d[i][j];
}

 public boolean MSymTriEigenv(double[][] a,double[][] Z,int n )
 {int Max=1000,i,j,k=n-1;
  double d,u;
  double[][] Q=new double[n][n];
  for(i=0;i<Max;i++)
  {if(a[k][k-1]==0)k--;
   if(k==0)break;
   d=(a[k-1][k-1]-a[k][k])/2;
   u=Math.sqrt(d*d+a[k-1][k]*a[k-1][k]);
   if(d<0)u=-u;
   u=a[k][k]+d-u;
   for(j=0;j<n;j++)a[j][j]-=u;
   QR(a,Q,n,n);
   Mrcheng(Z,Q,Z,n,n,n);
   Mrcheng(a,Q,a,n,n,n);
   for(j=0;j<n;j++)a[j][j]+=u;
  }
  if(i==Max+1)return false;
  else return true;
  }

public void paint(Graphics g)
{double[][] a={{10,1,2,3,4},{1,9,-1,2,-3},{2,-1,7,3,-5},{3,2,3,12,-1},{4,-3,-5,-1,15}};
 double[][] Q=new double[5][5];
 int i;
 MSymTri(a,Q,5);
 MSymTriEigenv(a,Q,5);
 g.drawString("特征值为：",10,10);
 for(i=0;i<5;i++)
 g.drawString(""+i+":"+a[i][i],10,i*10+20);
 g.drawString("特征向量为：",10,70);
 for(i=0;i<5;i++)
 g.drawString(""+i+":"+Q[0][i]+"  "+Q[1][i]+"  "+Q[2][i]+"  "+Q[3][i]+"  "+Q[4][i],10,80+10*i);
 }
 }