dlt2.cpp

利用高斯消元法解L系数

#include <iostream.h>
#include <math.h>
#include <stdio.h>

#define  DataType long double
#define N 12
int GaussEliminant(DataType a[N][N],DataType b[N],DataType c[N]);//利用高斯消元法解多元一次线性方程组

int main(int argc,char *argv[])
{ 
	DataType control[7][4]= // 6个控制点的(XYZ)坐标
	{
		{0,0,0,0},
		{0,1.311,1.751,1.951},
		{0,5.121,5.531,5.941},
		{0,10.511,10.751,10.191},
		{0,15.311,15.551,15.971},
		{0,20.711,20.551,20.981},
		{0,25.311,25.541,25.951},
	};
	DataType photo_control[7][3] = //6个控制点在图像内的坐标(XY)
	{
		{0,0,0},
		{0,1.211,1.531},
		{0,5.141,5.551},
		{0,10.161,10.751},
		{0,15.181,15.591},
		{0,20.121,20.541},
		{0,25.151,25.561},
	};
	DataType b[N];
	DataType L[N];
	for(int ii=0; ii<N; ii++)
	{
		b[ii] = 0;//b 是方程右边
		L[ii] = 0;//L 是L系数
	}
	
	DataType a[N][N];//相当于方程Al = -B 左侧的A a[0]行与列不参与运算
	for (ii=1; ii<=11; ii++)//初始化为0
		for(int jj=1;jj<=11;jj++)
			a[ii][jj] =0;

	//以下是给a[11][11]赋值过程
	for ( ii=1; ii<=11; ii++)//ii 行下标
		for( int jj=1;jj<=11;jj++)//jj 列下标
		{
			if (ii%2)//奇数行
			{
				if (jj<=3)//奇数行的前3项与实际控制点坐标相同
				{
					a[ii][jj] = control[ii/2+1][jj];
				}
				else if (jj==4)//奇数行的第4项为1
				{
					a[ii][jj] = 1;
				}
				else if (jj>4&&jj<8)//第4到第8项为0
				{
					a[ii][jj] = 0;
				}
				else if (jj>8&&jj<N)//第8到第11项为一乘积
				{
					a[ii][jj] = control[ii/2+1][jj-8]*photo_control[ii/2+1][1];
				}		
			}
			else if (ii%2 == 0)
			{
				if (jj<=4)
				{
					a[ii][jj] = 0;
				}
				else if (jj>4&&jj<8)
				{
					a[ii][jj] = control[ii/2][jj-4];
				}
				else if (jj==8)
				{
					a[ii][jj] = 1;
				}
				else if (jj>8&&jj<N)
				{
						a[ii][jj] = control[ii/2][jj-8]*photo_control[ii/2][2];
				}
			}				
		}//N维数组A赋值过程结束
		//给B数组赋值
		for ( ii=1;ii<N;ii++)
		{
			if (ii%2)
			{
				b[ii] = photo_control[(ii+1)/2][1];
			}
			else if(ii%2==0)
			{
				b[ii] = photo_control[ii/2][2];
			}
		}
		char szBuff[50];//临时数组用于格式化输出DOUBLE
		cout<<"方程左边矩阵A为：";
		for (ii=1;ii<N;ii++)
		{
			cout<<endl;
			for (int jj=1;jj<N;jj++)
			{
				sprintf(szBuff,"%f",a[ii][jj]);
				cout<<szBuff<<" ";
			}
		}
		cout<<endl<<"方程右边矩阵B为："<<endl;
		for (ii=0;ii<N;ii++)
		{
			cout<<b[ii]<<endl;
		}
		GaussEliminant(a,b,L);
		return 0;
}

int GaussEliminant(DataType a[N][N],DataType b[N],DataType c[N])
{
	DataType d=0,t=0,sum=0;
	int k=1,l=0,i=0,j=0,m,n;
	n=N-1;//n代表方程的组数 
	while(k<=n)
	{
		//*****选主元素*****
		d=a[k][k];
		l=k;
		for(i=k+1;i<=n;i++)
		{
			if(fabs(a[i][k])>fabs(d))
			{
				d=a[i][k];//现在d中存放的是：a[k]列中最大的系数 就是对角的那一列
				l=i;// l 中存放：最大数的行标
			}
		}
		if(l!=k)//如果对角元素在这一列中不是最大的
		{
			for(j=k;j<=n;j++)//则将对角元素所在式 与 最大元素所在式进行 系数 交换 //两式对调
			{					//目的是：使最大元素在对角线上
				t=a[l][j];
				a[l][j]=a[k][j];
				a[k][j]=t;
			}
			//下面三步是交换方程 等号 右面的系数
			t=b[k];
			b[k]=b[l];
			b[l]=t;
		}
		//*****消元*****
		for(j=k+1;j<=n;j++)
		{
			a[k][j]=a[k][j]/a[k][k];
		}
		b[k]=b[k]/a[k][k];
		for(i=k+1;i<=n;i++)
		{
			for(j=k+1;j<=n;j++)
			{
				a[i][j]=a[i][j]-a[i][k]*a[k][j];
			}
			j=1;
			b[i]=b[i]-a[i][k]*b[k];
		}
		k++;
	}
	//*****回代*****
	for(i=n-1;i>=1;i--)
	{
		sum=0;
		for(j=i+1;j<=n;j++)
		{
			sum=sum+a[i][j]*b[j];
		}
		b[i]=b[i]-sum;
	}
	cout<<endl<<endl<<"最后结果:"<<endl;
	for(int loop=1;loop<=n;loop++)
	{
		c[loop] = b[loop];
		cout<<"x"<<loop<<"="<<c[loop]<<endl;
	}
	return 0;
}