example4_10.java

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

import java.applet.*;
import java.awt.*;
import corejava.*;

public class Example4_10 extends Applet
{public static void F(double[] x,double[] y)
 {y[0]=x[0]*(x[0]+7*x[1])+3*x[1]*x[1]+0.5;
  y[1]=(x[0]-x[1])*(x[0]-x[1])-1;
  y[2]=x[0]+x[1]+1;
 }

 public static void FJacobi(double [] x,double[] [] a,int m,int n)
 {a[0][0]=2*x[0]+7*x[1];
  a[0][1]=7*x[0]+6*x[1];
  a[1][0]=2*x[0]-2*x[1];
  a[1][1]=-2*x[0]+2*x[1];
  a[2][0]=a[2][1]=1;
 }

 public static void NLGIN(double[] x,int m,int n)
 {int i,j,k;
  double alpha,H1=0,H2=0,z=0,y1=0,y2=0;
  double[][] a=new double[m][n];
  double [] b=new double[8];
  double [] y=new double[8];
  double [] v=new double[m];
  double [] u=new double[n];
  double [] x1=new double[n];
  alpha=1.0;
  while(true)
  {F(x,v);FJacobi(x,a,m,n);LEinv(a,v,u,m,n);
   j=-1;
   while(++j<8)
   { if(j<=2)z=alpha+j*0.01;
     else z=H2;
     for(i=0;i<n;i++)x1[i]=x[i]-z*u[i];
     F(x1,v);y1=0;
     for(i=0;i<m;i++)y1+=v[i]*v[i];
     for(i=0;i<n;i++)x1[i]-=0.00001*u[i];
     F(x1,v);y2=0;
     for(i=0;i<m;i++)y2+=v[i]*v[i];
     if(Math.abs(y2-y1)<1e-15)break;
     H1=(y2-y1)*1e5;H2=z;
     if(j==0){y[0]=H1;b[0]=H2;}
     else { y[j]=H1;k=0;
            for(i=0;i<j;i++)
            {if(k==0)
               {if(Math.abs(H2-b[i])+1.0==1.0)k=1;
                else H2=(H1-y[i])/(H2-b[i]);}
	        }
	        b[j]=H2;if(k!=0)b[j]=1e+35;
	        H2=0.0;
	        for(i=j-1;i>=0;i--)H2=-y[i]/(H2+b[i+1]);
	        H2+=b[0];
	       }
	}
    alpha=z;y1=y2=0;
    for(i=0;i<n;i++)
    {u[i]*=alpha;x[i]-=u[i];y1+=Math.abs(x[i]);y2+=Math.abs(u[i]);}
     if(y2/y1<1e-6)break;
   }

  }


  public void paint(Graphics g)
  {double[] x={0.5,-1.0};
   double[] y=new double[3];
   NLGIN(x,3,2);
   g.drawString("X0="+x[0],10,20);
   g.drawString("X1="+x[1],10,40);
   F(x,y);
  }




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

public static void  svd(double[][] a,double[][] U,double[][] V, int m,int n )
{int i,j,k,l,loop=0;
 double c,s,t,u,v,w,x,y;
 double[] b=new double[m]; double[] d=new double[m];
 boolean flag= true;
 for(i=0;i<m;i++)
    {for(j=0;j<m;j++)U[i][j]=0;
     U[i][i]=1;}
 for(i=0;i<n;i++)
	{for(j=0;j<n;j++)V[i][j]=0;
     V[i][i]=1;}
  //化为双对角部分
  for(k=0;k<m-1;k++)
  { //下三角部分
	c=0;
    for(l=k+1;l<m;l++)c+=a[l][k]*a[l][k];
    if(c!=0)
    {c=Math.sqrt(c+a[k][k]*a[k][k]);
    if(a[k][k]<0)c=-c;
    s=(c+a[k][k])*c;
    for(i=0;i<m;i++)
     {b[i]=U[i][k]*c;
      for(l=k;l<m;l++)b[i]+=U[i][l]*a[l][k];}
    for(i=0;i<m;i++)
    {U[i][k]-=b[i]/c;
	 for(j=k+1;j<m;j++)U[i][j]-=b[i]*a[j][k]/s;}
     for(j=0;j<n;j++)
     {b[j]=c*a[k][j];
      for(l=k;l<m;l++)b[j]+=a[l][k]*a[l][j];}
    for(j=k+1;j<n;j++)
     {a[k][j]-=b[j]/c;
      for(i=k+1;i<m;i++)a[i][j]-=a[i][k]*b[j]/s;}
    for(l=k+1;l<m;l++)a[l][k]=0;
   a[k][k]=-c;}
   //上三角次对角线外部分
   if(k>=n-2)continue;
   c=0;
   for(l=k+2;l<n;l++)c+=a[k][l]*a[k][l];
   if(c!=0)
   {c=Math.sqrt(c+a[k][k+1]*a[k][k+1]);
    if(a[k][k+1]<0)c=-c;
    s=(c+a[k][k+1])*c;
    for(j=0;j<n;j++)
    {b[j]=V[k+1][j]*c;
     for(l=k+1;l<n;l++)b[j]+=V[l][j]*a[k][l];}
    for(j=0;j<n;j++)
    {V[k+1][j]-=b[j]/c;
	 for(i=k+2;i<n;i++)V[i][j]-=a[k][i]*b[j]/s;}
    for(i=0;i<m;i++)
    {b[i]=a[i][k+1]*c;
     for(l=k+1;l<n;l++)b[i]+=a[i][l]*a[k][l];}
     for(i=k+1;i<m;i++)
     {a[i][k+1]-=b[i]/c;
     for(j=k+2;j<n;j++)a[i][j]-=b[i]*a[k][j]/s;}
     for(j=k+2;j<n;j++)a[k][j]=0;
    a[k][k+1]=-c;}
   }
  //对角化部分
   while(loop++<100)
   {for(k=0;k<n-1;k++)
    {u=Math.abs(a[k][k+1]);v=Math.abs(a[k][k])+Math.abs(a[k+1][k+1]);
     if(u+v==v)a[k][k+1]=0;
     else break;}
     if(k==n-1)break;
      t=a[k][k+1]/a[k][k];c=1/Math.sqrt(1+t*t);s=c*t;
      u=a[k][k]*c+a[k][k+1]*s;v=s*a[k+1][k+1];w=-c*a[k+1][k+1];
      a[k][k]=u;a[k][k+1]=0;a[k+1][k]=v;a[k+1][k+1]=w;
      for(j=0;j<n;j++){b[j]=c*V[k][j]+s*V[k+1][j];d[j]=s*V[k][j]-c*V[k+1][j];}
      for(j=0;j<n;j++){V[k][j]=b[j];V[k+1][j]=d[j];}
     for(i=k;i<n-2;i++)
     { if(a[i+1][i]!=0)
      {t=a[i+1][i]/a[i][i];c=1/Math.sqrt(1+t*t);s=c*t;
       u=c*a[i][i]+s*a[i+1][i];v=c*a[i][i+1]+s*a[i+1][i+1];w=s*a[i+1][i+2];
       x=s*a[i][i+1]-c*a[i+1][i+1];y=-c*a[i+1][i+2];
       a[i][i]=u;a[i][i+1]=v;a[i][i+2]=w;a[i+1][i]=0;a[i+1][i+1]=x;a[i+1][i+2]=y;
       for(j=0;j<m;j++){b[j]=U[j][i]*c+U[j][i+1]*s;d[j]=U[j][i]*s-U[j][i+1]*c;}
       for(j=0;j<m;j++){U[j][i]=b[j];U[j][i+1]=d[j];}}
      if(a[i][i+2]!=0)
     {t=a[i][i+2]/a[i][i+1];c=1/Math.sqrt(1+t*t);s=c*t;
      u=a[i][i+1]*c+a[i][i+2]*s;v=a[i+1][i+1]*c+a[i+1][i+2]*s;w=s*a[i+2][i+2];
      x=a[i+1][i+1]*s-c*a[i+1][i+2];y=-c*a[i+2][i+2];
      a[i][i+1]=u;a[i][i+2]=0;a[i+1][i+1]=v;a[i+1][i+2]=x;a[i+2][i+1]=w;a[i+2][i+2]=y;
      for(j=0;j<n;j++){b[j]=V[i+1][j]*c+s*V[i+2][j];d[j]=V[i+1][j]*s-c*V[i+2][j];}
      for(j=0;j<n;j++){V[i+1][j]=b[j];V[i+2][j]=d[j];}}
    }
   if(a[n-1][n-2]!=0)
   {t=a[n-1][n-2]/a[n-2][n-2];c=1/Math.sqrt(1+t*t);s=c*t;
    u=c*a[n-2][n-2]+s*a[n-1][n-2];v=c*a[n-2][n-1]+s*a[n-1][n-1];w=s*a[n-2][n-1]-c*a[n-1][n-1];
    a[n-2][n-2]=u;a[n-2][n-1]=v;a[n-1][n-2]=0;a[n-1][n-1]=w;
    for(j=0;j<m;j++){b[j]=U[j][n-2]*c+U[j][n-1]*s;d[j]=U[j][n-2]*s-c*U[j][n-1];}
    for(j=0;j<m;j++){U[j][n-2]=b[j];U[j][n-1]=d[j];}}
  }

   for(i=0;i<n;i++)
   {if(a[i][i]<0)
    {a[i][i]=-a[i][i];
     for(j=0;j<n;j++)V[i][j]=-V[i][j];}
   }
   for(i=0;i<n-1;i++)
   {k=i;
    for(j=i;j<n;j++){if(a[k][k]<a[j][j])k=j;}
    if(k!=i)
    {t=a[k][k];a[k][k]=a[i][i];a[i][i]=t;
     for(j=0;j<n;j++){t=V[k][j];V[k][j]=V[i][j];V[i][j]=t;}
     for(j=0;j<m;j++){t=U[j][k];U[j][k]=U[j][i];U[j][i]=t;}}
   }
}

 public static  void MInv(double[][] a,double[][] aa,int m,int n)
 {double[][] U=new double[m][m];
  double[][] V=new double[n][n];
  double[][] bb=new double[n][n];
  int i,j,k,p;double t;
  svd(a,U,V,m,n);
  for(k=n-1;k>=0;k--){if(a[k][k]!=0)break;}
  p=k+1;
  for(k=0;k<p;k++){bb[k][k]=1/a[k][k];}
  for(i=0;i<m;i++)
   for(j=i+1;j<m;j++)
   {t=U[i][j];U[i][j]=U[j][i];U[j][i]=t;}
  for(i=0;i<n;i++)
      for(j=i+1;j<n;j++)
   {t=V[i][j];V[i][j]=V[j][i];V[j][i]=t;}
  Mrcheng(bb,U,aa,p,p,m);Mrcheng(V,aa,aa,n,p,m);
 }

public static  void LEinv(double[][] a,double[] b,double[] x,int m,int n)
{double[][] aa=new double[n][m];
 int i,j;
 MInv(a,aa,m,n);
 for(i=0;i<n;i++)
 {x[i]=0;
  for(j=0;j<m;j++)x[i]+=aa[i][j]*b[j];}
  }

}