example4_9.java

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

import java.applet.*;
import java.awt.*;
public class Example4_9 extends Applet
{
  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 void F(double[] x,double[] y)
  {y[0]=x[0]*x[0]+x[1]*x[1]+x[2]*x[2]-1;
   y[1]=2*x[0]*x[0]+x[1]*x[1]-4*x[2];
   y[2]=3*x[0]*x[0]-4*x[1]+x[2]*x[2];
  }

  public void dnewton(double[] x,double h,double eps)
  {int i,j,n=x.length;
   double[] y=new double[n];
   double[] x1=new double[n];
   double[][] a=new double[n][n];
   double[][] z=new double[n][1];
   double[][] z1=new double[n][1];
    double max,temp;
   while(true)
   {F(x,y);max=0;
    for(i=0;i<n;i++)z1[i][0]=y[i];
    for(i=0;i<n;i++)
    {temp=Math.abs(y[i]);
    if(temp>max)max=temp;}
    if(max<eps)break;
    for(j=0;j<n;j++)
     {for(i=0;i<n;i++)x1[i]=x[i];
      x1[j]+=h;F(x1,y);
      for(i=0;i<n;i++)a[i][j]=y[i];
     }
   GaussSolve(a,z1,z,n,1);
   temp=1;
   for(i=0;i<n;i++)temp-=z[i][0];
   for(i=0;i<n;i++)x[i]-=h/temp*z[i][0];
   h*=0.1;
   }
 }

 public void paint(Graphics g)
 {double[] x={1,1,1};
  double h=0.1,eps=1e-7;
  dnewton(x,h,eps);
  for(int i=0;i<x.length;i++)
  {g.drawString("X"+i+"="+x[i],10,20+20*i);}
 }
}