问题
求解线性方程组
$$ \left{ \begin{array}{c} a_{11}x_1+a_{12}x_2+…+a_{1n}x_n=b_1 \ a_{21}x_1+a_{22}x_2+…+a_{2n}x_n=b_2 \ … \ a_{n1}x_1+a_{n2}x_2+…+a_{nn}x_n=b_n \end{array} \right. $$
其中方程组系数 $a_{ij}(i,j=1,2,3,…,n)$ 与右端项 $b_i(i=1,2,3,…,n)$ 均为实数，且 $b_1,b_2,…,b_n$ 不全为零，上述方程组简记为 $$ Ax=b $$ 其中 $$ A=\begin{bmatrix} a_{11} & a_{12} & … & a_{1n} \ a_{21} & a_{22} & … & a_{2n} \ …&…&…&… \ a_{n1} & a_{n2} & … & a_{nn} \ \end{bmatrix}, x=\begin{bmatrix} x_{1}\ x_{2}\ … \ x_{n}\ \end{bmatrix}, b=\begin{bmatrix} b_{1}\ b_{2}\ … \ b_{n}\ \end{bmatrix}. $$ 解线性方程组直接法（也称精确解法）的定义为：
若计算过程中没有舍入误差，则经过有限次的算术运算可以求得方程组的准确解的方法。
其主要有两类方法：
- Gauss 消去法
- Gauss 消去法
- Gauss 列主元消去法
- Gauss 按比例列主元消去法
- Gauss-Jordan 消去法
- 直接三角分解法
- Corut 方法
- $L^TDL$ 分解
- 追赶法
本篇只用 Matlab 实现上述算法，不涉及算法的介绍与证明，如果对证明与原理比较感兴趣，详情参考：
并且事实上这些算法在 Matlab 中都有既成的函数或者工具包，虽然 CS 里面有一句话：“如果有现成的轮子那么绝对不自己造。”很有道理，但是对于算法的理解还是需要亲自尝试，这样才能比较好的理解。
在最初编写这些算法时，我的 Matlab 代码风格还没有确定，以至相比于后续算法的源码，这些源码会显得比较的不规范，我会进行一部分的修正，但是毕竟时间比较久，修正程度可能比较有限，如果有一些问题，欢迎联系。
正文
Gauss 消去法
Gauss 消去法
function x = Gaussian_elimination( A, b )
% 高斯消去法解线性方程组
% Input: A: 系数矩阵; b: 方程组右端项;
% Output: x: 线性方程组解;
n = size(A, 1); % 矩阵 A 阶数
if size(A,1)~=size(A,2) || size(A,1)~=n
error('Matrix does not match')
end
for k = 1 : n - 1
if A(k, k) == 0 % 检测 A 是否为正定矩阵
for i = k + 1 : n
if A(i, k) ~= 0
for j = k : n
t = A(k, j); A(k, j) = A(i, j); A(i, j) = t;
end % for
else
error('A is singular.');
end % if
end % for
end % if
for i = k + 1 : n
l(i, k) = A(i, k) / A(k, k);
for j = k : n
A(i, j) = A(i, j) - l(i, k) * A(k, j);
end % for
b(i) = b(i) - A(i, k) * b(k);
end % for
end % for
x(n) = b(n) / A(n, n);
for k = n - 1 : -1 : 1
sum = 0;
for i = k + 1 : n
sum = sum + A(k, j) * x(j);
end % for
x(k) = (b(k) - sum) / A(k, k);
end % for
end % function
Gauss 消去法就非常的基础了，不过是平时我们解方程的固有过程，用来练习一下 Matlab 基本语言的使用。
当初在做这些课题的时候对代码的健壮性考虑的较少，主要是来理解算法的。如果有相应的需求，需要进一步的修改。
Gauss 列主元消去法
function x = Gaussian_principal_element_elimination( A, b )
% 高斯列主元消去法解线性方程组
% Input: A: 系数矩阵; b: 方程组右端项;
% Output: x: 线性方程组解;
n = size(A, 1); % 矩阵 A 阶数
if size(A, 1) ~= size(A, 2) || size(A, 1) ~= n
error('Matrix does not match')
end % if
for i = 1 : n
[ap, p] = max(abs(A(i : n, i)));
p = p + i - 1;
if ap == 0
error('divide by zero!')
end % if
if p > i
t = A(i, :);
A(i, :) = A(p, :);
A(p, :) = t;
t = b(i);
b(i) = b(p);
b(p) = t;
end % if
%消元计算
m = A(i+1:n, i) ./ A(i, i);
A(i+1:n, i+1:n) = A(i+1:n, i+1:n) - m * A(i, i+1:n);
b(i+1:n) = b(i+1:n) - m * b(i);
A(i+1:n, i) = zeros(n-i, 1);
end % for
x = zeros(n, 1);
x(n) = b(n) / A(n, n);
for i = (n - 1) : -1 : 1
x(i) = (b(i) - A(i, i+1:n) * x(i+1:n)) / A(i, i);
end % for
end % function
Gauss 按比例列主元消去法
function x = Gauss_proportional_column_principal_elimination( A, b )
% 高斯消去法解线性方程组
% Input: A: 系数矩阵; b: 方程组右端项;
% Output: x: 线性方程组解;
n = size(A, 1); % 矩阵 A 阶数
if size(A, 1) ~= size(A, 2) || size(A, 1) ~= n
error('Matrix does not match')
end
A_max = max(A, [], 2); % A 行最大值
for i = 1 : n
s(i) = A_max(i);
if s(i) == 0
error('A is singular.')
end % if
p(i) = i;
end % for
for k = 1 : n - 1
%%%%%%%%%% 待完成
r = 1; % 解 abs(a(p(r), k)) / s(p(r)) = max_k<=1<=n (abs(a(p(i), k)) / s(p(i));
%%%%%%%%%%
if A(p(r), k) == 0
error('A is singular.')
end % if
if k ~= r
temp = p(r);
p(k) = p(r);
p(r) = temp;
end % if
for i = k + 1 : n
A(p(i), k) = A(p(i), k) / A(p(k), k);
for j = k + 1 : n
A(p(i), j) = A(p(i), j) - A(p(i), k) * A(p(i), j);
end % for
end % for
end % for
b2(p(2)) = b(P(1));
for i = 2 : n
sum = 0;
for j = 1 : i - 1
sum = sum + A(p(i), j) * b2(p(j));
end % for
b2(p(i)) = b(p(i)) - sum;
end % for
if A(p(n), n) == 0
error('A is singular.')
end % if
x(n) = b2(p(n)) / A(p(n), n);
for i = n - 1 : -1 : 1
sum = 0;
for j = i + 1 : n
sum = sum + a(p(i), j) * x(j);
end % for
x(i) = (b2(p(i)) - sum) / a(p(i), i);
end % for
end % function
Gauss-Jordan 消去法
Gauss-Jordan 消去法就是在上述任一方法之中，消去步骤中对主元行外的所有行进行操作，其他没有任何变化，代码依方法选择而变，在此不进一步考虑。
直接三角形分解法
Crout 方法：
function x = Court( A, b )
% Court 方法解线性方程组
% Input: A: 系数矩阵; b: 方程组右端项;
% Output: x: 线性方程组解;
n = size(A);
L = zeros(n, n);
U = zeros(n, n);
for i = 1 : n
L(i, i) = 1;
end % for
for k = 1 : n
for j = k : n
U(k, j) = A(k, j) - sum(L(k, 1 : k - 1) .* U(1 : k - 1, j)');
end % for
for i = k + 1 : n
L(i, k) = (A(i, k) - sum(L(i, 1 : k - 1) .* U(1 : k - 1 , k)')) / U(k, k);
end % for
end % for
[L,U] = myLU(A);
[n,m] = size(A);
y(1) = b(1);
for i = 2 : n
for j = 1 : i - 1
b(i) = b(i) - L(i, j) * y(j);
end
y(i) = b(i);
end % for
x(n) = y(n) / U(n, n);
for i = (n - 1) : -1 : 1
for j = n : -1 : i + 1
y(i) = y(i) - U(i, j) * x(j);
end
x(i) = y(i) / U(i, i);
end % for
end % function
Court 方法是三角分解法中最简单的一种，计算量会大一些，但该方法适用于任意情况的方程组。
$L^TDL$ 分解法
$L^TDL$ 分解法适用面比较窄：必须是实对称矩阵。此时使用这种方法可以比较好的减少计算量，增加计算速度，但实际原理与 LU 分解法类似。
若是追求比较通用的方法，LDU 分解法是 LU 分解法的优化，但整体上大同小异，在此不占用篇幅专门书写它们的算法了。
追赶法
追赶法的适用范围也相当的窄：三对角实线性方程组。其本质就是 $LU$ 分解法在这里的运用，解这样的方程非常的迅速，需求的计算量也非常的小。
对于追赶法有一篇比较好的文章可以参考：Thomas, L.H., Elliptic Problems in Linear Differential Equations over a Network. Watson Science Computer Laboratory Report, 1949.
其中追赶法在国外常称为“Thomas 方法”，在热传导方程等物理方面有着非常广泛的用处。
结尾
对于实线性方程组，通过恰当的变换，都可以得到任意精确度的数值解或者直接获得解析解，在数值处理上属于最简单的一个类型。
在数值分析中，这些过程需要交给计算机计算，所以就需要着重考虑其中会由计算机计算产生的误差。因此诞生出了诸多的方法并不断优化，避免了 7种情况 带来的计算误差，使得数据处理更加精确可信。
数值分析也研究了常微分方程组、函数逼近、超越方程寻根等问题。在这些问题的算法中更加体现出了数值分析的精妙逻辑与魅力，我们使用精妙的方法通过计算机去一步步逼近真实结果，以此用于实际问题中的复杂方程求解。
这是非常有现实价值与研究价值的，之后我也会分别总结这些方面涉及到的一些算法，本次的线性方程组就是一个简单的入门。