hilbert.java

数值计算（科学计算）实验代码高斯消元楚列斯基分解

package hilbert;
import java.math.*;
public class Hilbert {
	/**
	 * @param args
	 */
	double[][]	hilArrary;
	double[] b ,xx;
	public static void main(String[] args) {
		Hilbert aa=new Hilbert();
		for(int t=1;;t++)
		{   
			double[][] tArrary=aa.hilbert(t);
			double[] bValue=aa.valueofb(t);
			double[][] lValue=aa.cholesky(tArrary,t);
			double[][] ltValue=aa.choleskyT(lValue, t);
		    double[] xx=aa.calculatex(t, lValue, ltValue, bValue);
			double[] xcha=aa.xwucha( xx, t);
			double sssssss=aa.man(xcha, t);
			if(sssssss/1>1) break;
			System.out.println("n="+t);
		    // TODO Auto-generated method stub
			for(int i=0;i<t;i++)
				{
				System.out.print(bValue[i]+"\t");
				}
			System.out.println();
			for(int i=0;i<t;i++)
			{
				System.out.print(xx[i]+"\t");
			}
			System.out.println();
			System.out.println(sssssss/1);
			System.out.println();
		}	
	}	
	public double[][] hilbert(int n)
	{
		hilArrary=new double[n][n];
		for(int i=0;i<n;i++)
		{
			for(int j=0;j<n;j++)
				hilArrary[i][j]=(double)1/(i+j+1);
		}
	return hilArrary;	
	}
	public double[] valueofb(int n)
	{  b=new double[n];
           for(int i=0;i<n;i++)
                b[i]=0;
		for(int i=0;i<n;i++)
			for(int j=0;j<n;j++)
				b[i]=hilArrary[i][j]+b[i];	
		return b;
	}
	public double[] valueofx(int n)
	{
		
		
	return xx;	
	}
	public double[][] cholesky(double[][] b,int n)
	{
		double[][] ll=new double[n][n];
                ll=b;
		for(int k=0;k<n;k++)
		{
			ll[k][k]=Math.sqrt(ll[k][k]);
			for(int i=k+1;i<n;i++)
			{
				ll[i][k]=ll[i][k]/ll[k][k];
			}
			for(int j=k+1;j<n;j++)
			{
				for(int i=k+1;i<n;i++)
					ll[i][j]=ll[i][j]-ll[i][k]*ll[j][k];	
			}
	
		}	
               for(int i=0;i<n;i++)
               {
                 for(int j=i+1;j<n;j++)
                 ll[i][j]=0;
                    }
		return ll;
	}

        public double[][] choleskyT(double[][] b,int n)
        { 
        	double[][] uu=new double[n][n];
        	
             for(int i=0;i<n;i++)
                 uu[i][i]=b[i][i];
             for(int j=0;j<n;j++)
 			{
 				for(int i=j+1;i<n;i++)
 					uu[i][j]=b[j][i];	
 			}
            return uu;
          }
        public double[] calculatex(int n,double[][] l,double[][] u,double[] b)
        {
        	double[] a=new double[n];
        	for(int j=0;j<n;j++)
        	{
        		if(l[j][j]==0) break;
        		a[j]=b[j]/l[j][j];
        		for(int i=j+1;i<n;i++)
        			b[i]=b[i]-l[i][j]*a[j];	
        	}
        	for(int j=n-1;j>=0;j--)
        	{
        		if(u[j][j]==0) break;
        		a[j]=b[j]/u[j][j];
        		for(int i=j+1;i<n;i++)
        			b[i]=b[i]-u[i][j]*a[j];		
        	}
        	
        	return a;
        }
        public double[] xwucha(double[] b,int n)
        {
        	double[] cha=new double[n];
        	for(int i=0;i<n;i++)
        	cha[i]=1-b[i];
        	return cha;
        }
        public double man(double[] a,int n)
        {
        	double cond=0;
        	cond=Math.abs(a[0]);
        	for(int i=0;i<n;i++)
        	{
        		if(cond<Math.abs(a[i])) cond=Math.abs(a[i]);
        	}
        	return cond;
        }
}