bijin.cpp

数值分析中实现对给定数表进行插值、曲面拟合、逼近

#include "BiJin.h"

//矩阵B和G初始化
void IniMatrixBG(POINT fxy[11][21],double **B,double **G,int k)
{
	int i,j;
	for (i=0; i<11; i++)
	{
		B[i] = (double*)malloc(k * sizeof(double));
		for (j=0; j<k; j++)
		{
			B[i][j] = pow(fxy[i][0].a1,j);
			//cout<<B[i][j]<<" ";
		}
	}
	for (i=0; i<21; i++)
	{
		G[i] = (double*)malloc(k * sizeof(double));
		for (j=0; j<k; j++)
		{
			G[i][j] = pow(fxy[0][i].a2,j);
		}
	}
}

//矩阵A(m,n)转置为B(n,m)
void TransMatrix(double **a, double **b, int m, int n)
{
	int i, j;
	for (i=0; i<n; i++)
	{		
		for (j=0; j<m; j++)
		{
			b[i][j] = a[j][i];
		}
	}
}

//求矩阵A(m*k),B(k*n)的乘积C(m*n)
void MultiMatrix(double **a, double **b, double **c, int m, int n, int k)
{
	int i, j, t;
	for (i=0; i<m; i++)
	{
		for (j=0; j<n; j++)
		{			
			c[i][j] = 0.0;
			for (t=0; t<k; t++)
			{
				c[i][j] = c[i][j] + a[i][t] * b[t][j];
			}
		}
	}
}

//求拟合的系数矩阵C
void ComputeMatrixC(POINT fxy[11][21], double **B, double **G,
					double **C,int m, int n, int k)
{
	int i, j;
	double **Bt, **Gt;	//用来存储矩阵B的转置(k行m列)和G的转置(k行n列)
	double **BtB, **GtG;//分别用来存储矩阵BT与B的乘积(k行k列)及矩阵GT与G的乘积(k行k列)
	double **u;		//U用来存储f(xi,yi),共m行n列
	double **BtuG;	//用来存储乘积BT*U*G，btug为k行k列
	double *x, *y;	//求解线性方程组Ax=y时所用的向量x，y
    int *M;			//符号数组，存储主元LU分解时的符号
	//建立动态矩阵空间
	Bt = (double**)malloc(k * sizeof(double));
	Gt = (double**)malloc(k * sizeof(double));
	BtB = (double**)malloc(k * sizeof(double));
	GtG = (double**)malloc(k * sizeof(double));
	u = (double **)malloc(m * sizeof(double));
	BtuG =(double **)malloc(k * sizeof(double));
	x = (double *)malloc(k * sizeof(double));
	y = (double *)malloc(k * sizeof(double));
    M = (int *)malloc(k * sizeof(int));
	for (i=0; i<k; i++)	
	{
		Bt[i] = (double*)malloc(m * sizeof(double));
		Gt[i] = (double*)malloc(n * sizeof(double));
		BtB[i] = (double*)malloc(k * sizeof(double));
		GtG[i] = (double*)malloc(k * sizeof(double));
		BtuG[i] = (double*)malloc(k * sizeof(double));
	}
	for (i=0; i<m; i++)
	{
		u[i] = (double *)malloc(n * sizeof(double));
		for (j=0; j<n; j++)
		{
			u[i][j] = fxy[i][j].a3;
		}
	}

	//求矩阵B和G的转置Bt和Gt
	TransMatrix(B,Bt,m,k);
	
	TransMatrix(G,Gt,n,k);
	
	//求BtB,GtG
	MultiMatrix(Bt,B,BtB,k,k,m);
	
	MultiMatrix(Gt,G,GtG,k,k,n);

    //下面求BT*U*G
	//先求BT*U，结果放在GT(k行n列)中
	MultiMatrix(Bt,u,Gt,k,n,m);

	//再求GT*G，结果放在BTUG中
	MultiMatrix(Gt,G,BtuG,k,k,n);
//		for (i=0;i<k;i++)
//	{
//		for (j=0;j<k;j++)
//		{
//			cout<<"btb"<<BtuG[i][j]<<" ";
//		}
//		cout<<endl;
//	}

	//求矩阵C，先对BTB作选主元LU分解
	LU(BtB,M,k);
	//做Doolittle分解，解出CGtG,付给C(k*k)
	for (j=0; j<k; j++)
	{
		//初始化b
		for (i=0; i<k; i++)
		{
			y[i] = BtuG[i][j];
		}
		
		Doolittle(BtB,y,M,x,k);
//		for (i=0;i<k;i++)
//		{
//			cout<<"x"<<x[i]<<" ";
//		}
//		
		for (i=0; i<k; i++) 
			C[i][j] = x[i];
		
	}
//	for (i=0;i<k;i++)
//	{
//		for (j=0;j<k;j++)
//		{
//			cout<<"c"<<C[i][j]<<" ";
//		}
//		cout<<endl;
//	}
	//把求(BTB(-1)*BTUG)的转置，赋给Btb
	TransMatrix(C,BtB,k,k);
//		for (i=0;i<k;i++)
//	{
//		for (j=0;j<k;j++)
//		{
//			cout<<"btb"<<BtB[i][j]<<" ";
//		}
//		cout<<endl;
//	}

	///对GTG作选主元LU分解
	LU(GtG,M,k);
	//做Doolittle分解，解出最终结果的转置,付给BTUG
	for (j=0; j<k; j++)
	{
		//初始化b
		for (i=0; i<k; i++) 
			y[i] = BtB[i][j];	
		Doolittle(GtG,y,M,x,k);
//		for (i=0;i<k;i++)
//		{
//			cout<<"x"<<x[i]<<" ";
//		}
		for (i=0; i<k; i++) 
			BtuG[i][j] = x[i];					
	}
	//最后把BTUG的转置付给C
	TransMatrix(BtuG,C,k,k);
//			for (i=0;i<k;i++)
//	{
//		for (j=0;j<k;j++)
//		{
//			cout<<"c"<<C[i][j]<<" ";
//		}
//		cout<<endl;
//			}


}