高斯消元法

最近的一个大作业中需要用到最小二乘拟合曲线。最小二乘法实际上就是求解一个线性方程组得到拟合多项式的系数。线性代数中常用的策略就是高斯消元，即满足:

1. 系数矩阵每一行第一个非0元素是主对角线元素
2. 系数矩阵主对角线元素均为1.

在这个条件下，高斯消元的求解过程可以描述为:
设系数矩阵A为$n\times n$矩阵，常数矩阵b为 $n\times 1$矩阵，增广矩阵为 $(A|b)$

//逐行得到主值为1的对角线系数矩阵,最后一行元素特判（运算过程中操作的是增广矩阵)
for i=0 to n-2 do begin
    //寻找i to n-1行中所有第i列绝对值最大的元素所在行号
	pivot=elementseek(i);
	//如果行号与当前行不同，则交换这两行的全部元素
	if (pivot<>i) then
		swap(pivot,i);
	if (A[i][i]==0) then begin
		writeln(系数矩阵的秩r(A)<n，方程可能无解或有无穷多解);
		break;
	end
	simple(i); //对第i行元素进行归一化使主对角线上的系数值为1
	for j=i+1 to n-1 do begin
		//对后面的每一行元素消元，使得后面每一行的第i列元素为0。
		Elimation(i,j);
	end
end
//消元结束后最后一行仅最后一列的系数非0，可以直接运算得到一个解，重新赋给A[n-1][n-1]是为了回代前面的行运算方便。
A[n-1][n-1]=b[n-1]/A[n-1][n-1];
for i=n-2 downto 0 do begin
	for j=i+1 to n-1 do begin
		tot=tot+A[i][j]*A[j][j]; //求出已知部分的解与方程系数相乘的累加和。
	end
	A[i][i]=b[i]-tot; //将累加和右移到等号右边得到第i个元素的解。
end

for i=0 to n-1 do begin
	writeln(A[i][i]); //输出每个方程的解。
end

#include<iostream>
#define maxn 10005
using namespace std;
double a[101*maxn],b[maxn];
int n;

double abss(double x){
	if (x<0)
		return -x;
	return x;
}

int SelectPivotalElement(int c){
	double tmp=a[c],ans=0;
	for (int i=0; i<n; i++){
		if (abss(a[i*n+c])>tmp){
			tmp=abss(a[i*n+c]);
			ans=i;
		}
	}
	return ans;
}

void SwapRow(int row1,int row2){
	double tmp;
	for (int i=0; i<n; i++){
		tmp=a[row1*n+i];
		a[row1*n+i]=a[row2*n+i];
		a[row2*n+i]=tmp;
	}
	tmp=b[row1]; b[row1]=b[row2]; b[row2]=tmp;
	return;
}

void SimplePivotalRow(int r,int c){
	double x=a[r*n+c];
	for (int i=0; i<n; i++){
		a[r*n+i]/=x;
	}
	b[r]/=x;
	return;
}

void RowElimination(int r1,int r2){
	double x=a[r2*n+r1];
	for (int i=r1; i<n; i++){
		a[r2*n+i]-=a[r1*n+i]*x;
	}
	b[r2]-=b[r1]*x;
	return;
}

int main(){
	cin>>n;
	for (int i=0; i<n; i++){
		for (int j=0; j<n; j++){
			scanf("%lf",&a[i*n+j]);
		}
		scanf("%lf",&b[i]);
	}
	
	for (int i=0; i<n-1; i++){
		int pivotRow=SelectPivotalElement(i);
		if (pivotRow!=i){
			SwapRow(i,pivotRow);
		}
		if (a[i*n+i]==0){
			cout<<"No Solution"<<endl;
			return 0;
		}
		SimplePivotalRow(i,i);
		for (int j=i+1; j<n; j++){
			RowElimination(i,j);
		}
	}

	if (abss(a[(n-1)*n+n-1])<1e-17){
		cout<<"No Solution"<<endl;
		return 0;
	}
	
	a[(n-1)*n+n-1]=b[n-1]/a[(n-1)*n+n-1];
	for (int i=n-2; i>=0; i--){
		double totcof=0.0;
		for (int j=i+1; j<n; j++){
			totcof+=a[i*n+j]*a[j*n+j];
		}
		a[i*n+i]=(b[i]-totcof);///a[i*n+i];
	}
	
	for (int i=0; i<n; i++){
		b[i]=a[i*n+i];
		printf("%.2lf\n",b[i]);
	}
	return 0;
}