matlab & numeric analysis basics_01

写在前面：刚刚开始使用这个平台，第一篇记录，还不是很习惯，如有问题，还望指出。本身是学生，专业不是学计算机，方向上属于金融。

在数值分析课程中，解线性代数方程组主要有迭代法和直接法。直接法具体有LU分解，Gauss消去法，主元素消去法等。这篇主要说一下Gauss消去法和列主元素消去法的matlab实现。

其实，对于解方程而言，matlab是可以直接进行方程组的解的，即利用矩阵除的命令。

假设系数矩阵是A，右端项矩阵是b，原方程组可以写成Ax = b，直接输入x = A\b命令可以得到线性代数方程组的解。当原方程组具有唯一解时输出唯一解，输出的x是唯一解；若原方程组“无解”，即超定方程组时，这个命令返回方程组最小二乘意义下的最优解；若原方程组有无穷多解时（欠定方程组），返回无穷多解的其中一组解。

Gauss消去法解线性方程组

Gauss消去法是直接法解线性代数方程组的一种原理较为容易的方法。在数值分析的入门课程中大都有所描述，具体原理如下：

step1:取方程组的系数矩阵A和右端项，组成增广矩阵[A b]，对增广矩阵进行初行变换使之变成阶梯形矩阵。
例如，以三阶系数矩阵为例：

$\left[\begin{array}{cccc}{a}_{11}& a_{12}& a_{13}& b_1\\ a_{21}& a_{22}& a_{23}& b_2\\ a_{31}& a_{32}& a_{33}& b_3\end{array}\right]\stackrel{初等行变换}{\Rightarrow }\left[\begin{array}{cccc}{a}_{11}^{\prime }& a_{12}^{\prime }& a_{13}^{\prime }& b_1^{\prime }\\ & a_{22}^{\prime }& a_{23}^{\prime }& b_2^{\prime }\\ & & a_{33}^{\prime }& b_3^{\prime }\end{array}\right]$

step2:通过阶梯形矩阵的最后一行可以算出最后一个未知数的值。

$x_3=b_3^{\prime }/a_{33}^{\prime }$

step3:最后一个值回代，可以算出每一个未知数的值。

$x_2=(b_2^{\prime }-a_{23}^{\prime }\ast x_3)/a_{22}^{\prime }$ 其中，$x_3$的值已知。
$x_1=(b_1^{\prime }-a_{13}^{\prime }\ast x_3-a_{12}^{\prime }\ast x_2)/a_{11}^{\prime }$ 其中 $x_2$，$x_3$的值已经解出。

Gauss消去法核心在于化成上三角形矩阵（阶梯形矩阵）和回代。

具体的Matlab算法实现：

function S = SolveEquationG(A, b)
% create the augmented matrix combined coefficient matrix and right items
B = [A b];
% the num of dimensions of the coefficient matrix (the number of rows)
n = size(A, 1);
% transfer the coefficient matrix to an upper triangular matrix
for k = 1 : n-1
    if B(k, k) ~= 0 %only when this element is not zero, the equations system could be solved in this method.
        for i = k+1 : n
            B(i, k) = B(i, k)/B(k, k);
            for j = k+1 : n+1
                B(i, j) = B(i, j) - B(i, k) * B(k, j);
            end
        end
    else break
    end
end
% reprocess the matrix and solve out the solution vector.
S = rand(n, 1);% create a vector to store the solutions.
S(n, 1) = B(n, n+1)/B(n, n); % calculate the 'xn'
%calculate rest elements of the solution vector.
for t1 = n-1:-1:1
    for t2 = n:-1:t1+1
        B(t1, n+1) = B(t1, n+1) - S(t2, 1)*B(t1, t2);
    end
    S(t1, 1) = B(t1, n+1)/ B(t1, t1);
end

列主元素消去法解方程组

Gauss消去法确实提供了一个解线性代数方程组的一个逻辑较为简单的方法，但是Gauss消去法却有着弊端——有可能会产生较大的舍入误差。如果消元过程中，第k个系数矩阵第k行第k列的元素是一个极小的数（eg. $0.4\ast 10^{-11}$) 。这时，这一行下一个元素的消元运算就有可能被“吃掉”，从而使最后计算出来的解向量产生一个较大的舍入误差。

为了避免舍入误差，需要对Gauss消去法进行改进，便引出列主元素消去法。（事实上，使用主元素消去法也可以达到同样的效果，但是显然算法上，列主元素消去法更容易实现）

列主元素消去法的算法过程大致和Gauss消去法一致，核心思想都是先处理系数矩阵与右端项列向量组成的增广矩阵，将增广矩阵化成上三角形矩阵（阶梯型矩阵）。

改进在于，每一列元素消元的过程中，将最大的元素换到第k行，第k列，然后消去下面的几行的第k个元素。

以进行第二列的消元过程为例：

$\left[\begin{array}{cccccc}{a}_{11}& a_{12}& ...& & a_{1n}& b_1\\ & a_{22}& ...& & a_{2n}& b_2\\ & a_{32}& ...& & a_{3n}& b_3\\ & ...& ...& & ...& ...\\ & a_{m2}& ...& & a_{mn}& b_m\end{array}\right]\underset{exchange{a}_{i2}with{a}_{22}}{\overset{assume{a}_{i2}isthebiggestelementofcolumn2}{\Rightarrow }}\left[\begin{array}{cccccc}{a}_{11}& a_{12}& ...& & a_{1n}& b_1\\ & a_{i2}& ...& & a_{in}& b_i\\ & a_{32}& ...& & a_{3n}& b_3\\ & ...& ...& & ...& ...\\ & a_{22}& ...& & a_{2n}& b_2\\ & ...& ...& & ...& ...\\ & a_{m2}& ...& & a_{mn}& b_m\end{array}\right]\stackrel{eliminatetheelementsofcolumn2}{\Rightarrow }\left[\begin{array}{cccccc}{a}_{11}& a_{12}& ...& & a_{1n}& b_1\\ & a_{i2}& ...& & a_{in}& b_2\\ & & a_{33}^{\prime }& ...& a_{3n}^{\prime }& b_3^{\prime }\\ & & a_{43}^{\prime }& ...& a_{4n}^{\prime }& b_4^{\prime }\\ & & ...& ...& ...& ...\\ & & a_{m3}^{\prime }& ...& a_{mn}^{\prime }& b_m^{\prime }\end{array}\right]$

按照这个过程进行消去可以有效减少最后算出的方程组解矩阵的舍入误差。

其余的过程类似于未改进的Gauss消去法的过程。
具体算法如下：

function S = SolveEquationM(A, b)
% create the augmented matrix combined the coefficiet matrix and right items
B = [A b];
% the num of dimensions of the coefficient matrix;
n = size(A, 1);
% transfer the coefficient matrix to an upper triangular matrix
for k = 1 : n
    IM = B(:, k);
% select the max value of this column as main element
    [max_IM, index] = max(IM, [], 1);
% exchange the value of rows
    B([k, index], :) = B([index k], :);
    if B(k, k) ~= 0
        for i = k+1 : n
            B(i, k) = B(i, k)/B(k, k);
            for j = k+1 : n+1
                B(i, j) = B(i, j) - B(i, k) * B(k, j);
            end
        end
    else break
    end
end
% reprocess the matrix (same as Gauss-elimination)
S = rand(n,1);
S(n, 1) = B(n, n+1)/B(n, n);
for t1 = n-1 : -1 : 1
    for t2 = n : -1 : t1 + 1
        B(t1, n+1) = B(t1, n+1) - S(t2, 1) * B(t1, t2);
    end
    S(t1, 1) = B(t1, n+1)/ B(t1, t1);
end