1.      非齊次線性方程組有解的條件

如下非齊次線性方程組：

由係數矩陣和常數列向量構成的增廣矩陣如下：

無解情況：

唯一解情況：

無窮解情況：

2.      高斯消元法求解

步驟：

1）    消元法

通過矩陣的初等變換，將增廣矩陣變換為上三角矩陣

2）    回代法

採用回代法求解上三角矩陣對應的非齊次線性方程組，即從第n個方程開始求解，解出

然後再將解出之代入第n-1個方程解出，…，直到解出所有來。

       注意到在消元和回代的過程中均需使用矩陣A的主對角線元素(稱為主元素)作除數，因此如果原方程組的某個主對角線元素等於0，就會導致求解失敗。而且從數值計算的特點可知，即使某個主元素並未等於0，而只是其絕對值很小，也會導致較大的計算誤差。因此，我們可以通過調換方程組中某些方程的相對位置來保證主對角線上的元素取較大的絕對值，這就是列主元素法。

3.     演算法實現

考慮簡單的情況，求解如下：

方程。

C語言程式碼：

/*
*
* Copyright (c) 2011,武漢大學國家多媒體中心
* All rights reserved.
*
* 檔名稱：SolutionLinearEquation.cpp
* 摘要：列主元素高斯消元法解n元一次方程組
*
* 當前版本：1.0
* 作者：王汪
* 完成日期：2011年9月10日
*
*/

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

/*
 * 功能:選取列主元素
 * 輸入:A[]		係數矩陣A；
 *		B[]		常數列向量B；
 *		A_Rows	係數矩陣A的行數；
 *		Kst_Row	待求的第k行列主元素
 * 輸出:void
 */
void ColumnPrimaryElement(double A[], double B[], int A_Rows, int Kst_Row)
{
	// 用於存放列主元素的值
	double main_element = 0.0;
	// 用於存放列主元素所在行
	int main_line = 0;

	// 中間變數，用於交換
	double temp = 0.0;
	// 迴圈變數
	int i, j;

	// 暫定A[k, k]為列主元素
	main_element = A[Kst_Row * A_Rows + Kst_Row];
	main_line = Kst_Row;

	for (i = Kst_Row + 1; i < A_Rows; ++i)
	{
		if (fabs(A[i * A_Rows + Kst_Row]) > fabs(main_element))
		{
			main_element = A[i * A_Rows + Kst_Row];
			main_line = i;
		}

		// 如果第k列元素中絕對值最大的不是a[k，k]，則交換兩個方程
		if (main_line != Kst_Row)
		{
			for (j = Kst_Row; j < A_Rows; ++j)
			{
				temp = A[Kst_Row * A_Rows + j];
				A[Kst_Row * A_Rows + j] = A[main_line * A_Rows + j];
				A[main_line * A_Rows + j] = temp;
			}

			temp = B[Kst_Row];
			B[Kst_Row] = B[main_line];
			B[main_line] = temp;
		}
	}
}

/*
 * 功能:列主元素高斯消元法解n元一次方程組
 * 輸入:A[]		係數矩陣A；
 *		B[]		常數列向量B；
 *		A_Rows	係數矩陣A的行數；
 * 輸出:X[]		結果向量
 */

void Gauss_ColumnPrimaryElement(double A[], double B[], int A_Rows, double X[])
{
	// 迴圈變數
	int i, j, k;
	// 高斯消元比例因子
	double c;

	// 消元
	for (k = 0; k < A_Rows; ++k)
	{
		// 選取列主元素
		ColumnPrimaryElement(A, B, A_Rows, k);
		for (i = k + 1; i < A_Rows; ++i)
		{
			// 求得高斯消元比例因子
			c = A[i * A_Rows + k] / A[k * A_Rows + k];

			for (j = k + 1; j < A_Rows; ++j)
			{
				A[i * A_Rows + j] = A[i * A_Rows + j] - A[k * A_Rows + j] * c;
			}
			B[i] = B[i] - B[k] * c; 
		}
	}

	// 有解條件判斷
	// 係數矩陣A的秩等於A的維數n(即行數或者列數)
	if (fabs(A[(A_Rows - 1) * A_Rows + (A_Rows - 1)]) < 10e-6)
	{
		// 不存在唯一解
		printf("不存在唯一解!\n");

		for (i = 0; i < A_Rows; ++i)
		{
			X[i] = 0.0;
		}
		return;
	}

	// 回代求解
	for (i = A_Rows -1; i >= 0; --i)
	{
		X[i] = B[i];
		for (j = i + 1; j < A_Rows; ++j)
		{
			X[i] = X[i] - A[i * A_Rows + j] * X[j];
		}
		X[i] = X[i] / A[i * A_Rows + i];
	}
}

void main()
{
	double A[4][4] = 
	{
		{0.2368, 0.2471, 0.2568, 1.2671},
		{0.1968, 0.2071, 1.2168, 0.2271},
		{0.1581, 1.1675, 0.1768, 0.1871},
		{1.1161, 0.1254, 0.1397, 0.1490},
	};

	double B[4] = 
	{
		1.8471, 1.7471, 1.6471, 1.5471
	};

	double X[4];

	int i = 0;

	Gauss_ColumnPrimaryElement(A[0], B, 4, X);

	printf("The result is x = \n");
	for (i = 0; i < 4; ++i)
	{
		printf("%10.6lf\n", X[i]);
	}
}

// MATLAB測試
// A = [[0.2368, 0.2471, 0.2568, 1.2671];[0.1968, 0.2071, 1.2168, 0.2271];[0.1581, 1.1675, 0.1768, 0.1871];[1.1161, 0.1254, 0.1397, 0.1490]]
// B = [1.8471; 1.7471; 1.6471; 1.5471]
// A \ B