test1_2.java.bak

数值分析之计算方法与实验5：曲线拟合 by java

//指数曲线拟合程序：
class Test1_2{
	static int n=0,p=1;//p表示拟合次数
	static int row=0,line=0;//解方程时ab矩阵的行数和列数
	public static void main(String[] args){
		double[] x={1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8}; 
		double[] y={33.4,79.5,122.65,159.05,189.15,214.15,238.65,252.5,267.55,280.50,296.65,301.40,310.40,318.15,325.15};
		
		//double[] x={1,2,3,4};double[] y={7,11,17,27};
		//double[] x={2.2,2.7,3.5,4.1,4.8};double[] y={65,60,53,50,46};
		n=x.length;
		for(int i=0;i<n;i++)y[i]=Math.log(y[i]);
		double[][] fei=new double[p+1][n];
		fei=qiufei(x);//将fei矩阵求出
		//将s矩阵求出
		double[][] s=new double[p+1][p+1];
		for(int i=0;i<=p;i++){
			for(int j=0;j<=i&&j<=p;j++){
				s[i][j]=neiji(fei[i],fei[j]);
				s[j][i]=s[i][j];
			}
		}
		/*for(int i=0;i<=p;i++){
			for(int j=0;j<=p;j++){
				System.out.print(s[i][j]+"   ");
			}
			System.out.println();
		}*/
		
		//将t矩阵求出
		double[] t=new double[p+1];
		for(int i=0;i<=p;i++){
			t[i]=neiji(y,fei[i]);
			//System.out.print(t[i]+"   ");
		}
		
		//将s矩阵和t矩阵合并成ab矩阵
		double[][] ab=new double[p+1][p+2];
		for(int i=0;i<=p;i++){
			for(int j=0;j<=p;j++)ab[i][j]=s[i][j];
		}
		for(int i=0;i<=p;i++)ab[i][p+1]=t[i];
		row=p+1;line=p+2;
		/*for(int i=0;i<=p;i++){
			for(int j=0;j<=p+1;j++){
				System.out.print(ab[i][j]+"   ");
			}
			System.out.println();
		}*/
		
		
		double[] a=directLU(ab);
		/*for(int i=0;i<=p;i++){
			System.out.println("a["+i+"]="+a[i]);
		}*/
		double aa=Math.exp(a[0]),bb=a[1];
		System.out.println("结果a="+aa+"   b="+bb);		
	}
	
	
	
	static double neiji(double[] x,double[] y){
		double s=0;
		for(int i=0;i<x.length;i++)s+=x[i]*y[i];
		return s;
	}	
	static double[][] qiufei(double[] x){
		double[][] fei=new double[p+1][n];
		for(int k=0;k<=p;k++){
			for(int j=0;j<n;j++)fei[k][j]=Math.pow(x[j],k);
		}
		return fei;
	}
	//改进平方根法解方程
	static double[] directLU(double[][] ab){
		double[][] u=new double[row][line];
		double[][] l=new double[row][row];
		for(int j=0;j<line;j++)u[0][j]=ab[0][j];
		for(int i=1;i<row;i++)l[i][0]=ab[i][0]/ab[0][0];
		for(int k=1;k<row;k++){
			//将U矩阵算出
			for(int j=k;j<line;j++){
				u[k][j]=ab[k][j];
				for(int q=0;q<=k-1;q++){
					u[k][j]=u[k][j]-u[q][j]*l[k][q];
				}
			}
			//算出矩阵L
			for(int i=k+1;i<row;i++){
				l[i][k]=u[k][i]/u[k][k];
			}
		}
		double[] x=new double[row];
		x[row-1]=u[row-1][line-1]/u[row-1][line-2];
		for(int i=row-2;i>=0;i--){
			x[i]=u[i][line-1];
			
			for(int j=i+1;j<=line-2;j++){
				x[i]=x[i]-u[i][j]*x[j];
			}
			x[i]=x[i]/u[i][i];
		}
		
		return x;	
	}
}