工程计算——奇异值分解法&迭代法

本次整理主要参考工程数学计算方法(第二版)_张诚坚_何南忠_覃婷婷_高教社及MATLAB数值计算实战_占海明_机械工业出版社两本书籍。

由于迭代法具有：

• 存储单元少
• 程序简单
• 计算速度快
  的优点，因此在实际大规模科学工程计算中迭代法用的更多。

例题引入

给定线性方程组
$$\left(\begin{array}{ccccc}1& 2& 3& 4& 5\\ -2& 3& 4& 5& 6\\ -3& -4& 5& 6& 7\\ -4& -5& -6& 7& 8\\ -5& -6& -7& -8& 9\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\end{array}\right)=\left(\begin{array}{c}55\\ 66\\ 63\\ 36\\ -25\end{array}\right)$$
试以$\Vert \begin{array}{c}{X}^{(k+1)}-X^{(k)}\end{array}\Vert$ $<$ $10^{-5}$作为终止准则，分别利用

• Jacobi迭代法
• Gauss-Seidel迭代法
• 具有最佳松弛因子的JOR迭代法
• SOR迭代法

求解该方程组，并比较这些方法的计算时间。

以上呈现的是书本上的原始题目，但通过求解该线性方程组的条件数$\text{cond}(A)$可以得到$\text{cond}(A)\gg 1$,也即：该线性方程组是病态的，这是不争的事实。因此，如果真的按照题目要求解题，那么求出的结果一定是发散的，这是没有太大的实际意义的，也更无必要通过这个例子去比较这些方法的计算时间了。对于求解病态的线性方程组我们可以使用——奇异值分解法（SVD）
关于条件数的求解可见工程计算——求解线性方程组
我们先来介绍SVD分解法：

SVD分解法

对于$n$阶实矩阵$A$,存在$n$阶正交阵$U、V$和对角阵$S$，使得$$A=US{V}^T$$
结合：
$$\begin{gathered} A=US{V}^T\Rightarrow A^{-1}=V{S}^{-1}{U}^T \\ Ax=b \end{gathered}$$
得到：
$$x=A^{-1}b=V{S}^{-1}{U}^{T}b(1)$$
假设：
$$\begin{gathered} U=[u_1,u_2,u_3,...,u_n] \\ V=[v_1,v_2,v_3,...,v_n] \\ S=diag({\sigma }_1,{\sigma }_2,...,{\sigma }_n) \end{gathered}$$
结合$(1)$得到：
$$x=\sum _{i=1}^n\frac{u_i^{T}b}{{\sigma }_i}{v}_i$$
因为，当${\sigma }_i$接近于$0$时，一方面会对计算$x$产生较大的误差；另一方面，由$A=US{V}^T=\sum _{i=1}^n{\sigma }_i{u}_i{v}_i^T$可知，当 ${\sigma }_i$ 接近于0时，后面的项对矩阵$A$的贡献率非常小。取一正的阈值（假设为 $ϵ$）当 $|{\sigma }_i|<ϵ$ 时，其相应的项就不再计算，即$$A\approx \sum _{|{\sigma }_i|\ge ϵ}{\sigma }_i{u}_i{v}_i^T$$
由此我们可以进一步得到：$$x\approx \sum _{|{\sigma }_i|\ge ϵ}\frac{u_i^{T}b}{{\sigma }_i}{v}_i$$
根据上式编写MATLAB程序：

function varargout = svd_solve(A,b,ep)
%SVD_SOLVE 奇异值分解法求解病态方程组
n = length(b);                  % 线性方程组的阶数
if nargin < 3 
    ep = 1e-10;                 % 默认精度
end
[U,S,V] = svd(A);               %奇异值分解
%MATLAB中提供了svd函数实现矩阵的奇异值分解，该函数的调用格式即为
%[U,S,V] = svd(A)
%其中U、V为n阶正交阵，V为对角阵
sigma = diag(S);                %提取对角线元素
x = zeros(n,1);                 %预置零向量存储迭代解
k = 1;                          %计数器
while 1
    if abs(sigma(k))<ep         %若不满足条件
        A = A-sigma(k)*U(:,k)*V(:,k)';%从A中剔除
    else 
        x = x+(U(:,k)'*b)/sigma(k)*V(:,k);% 累加
    end
    k = k+1;                    %计数器累加
    if k > n
        break                   % 判断次数大于方程组的阶数退出循环
    end
end
[varargout{1:2}] = deal(x,...
    A);

对于该题求解：
主函数可为：

A = [1 2 3 4 5 ;-2 3 4 5 6;-3 -4 5 6 7;-4 -5 -6 7 8;-5 -6 -7 -8 9];
b = [55;66;63;36;-25];
n = length(b);
x0 = zeros(n,1);
ep = 1e-5;
maxiter = 100;
[x,exitflag,iter,xs,err]= jacobic(A,b,x0,ep,maxiter);
if exitflag == 1
    disp("该线性方程组迭代收敛，其对应的解为：");
    x
    disp("迭代次数为：");
    iter
    disp("迭代序列为：");
    xs
    disp("近似相对误差为：")
else
    disp("该线性方程组迭代发散，改为使用奇异值分解法分解求解。");
    disp("奇异值分解法下：")
    [x,A]= svd_solve(A,b,ep);
    disp("利用奇异值分解法得到的该线性方程组的解为:");
    x
    disp("利用奇异值分解法得到的该线性方程组对应的近似系数矩阵为:");
    A
end

结果：

该线性方程组迭代发散，改为使用奇异值分解法分解求解。
奇异值分解法下：
利用奇异值分解法得到的该线性方程组的解为:

x =

    1.0000
    2.0000
    3.0000
    4.0000
    5.0000

利用奇异值分解法得到的该线性方程组对应的近似系数矩阵为:

A =

     1     2     3     4     5
    -2     3     4     5     6
    -3    -4     5     6     7
    -4    -5    -6     7     8
    -5    -6    -7    -8     9

迭代法

接下来介绍迭代法（iteration method）

Jacobi迭代法

$$\begin{gathered} B_J=-D^{-1}(L+U) \\ {F}_J=D^{-1}b \end{gathered}$$

function varargout = jacobic(A,b,x0,ep,maxiter)
%JACOBI 分量形式Jacobi迭代法求解线性方程组
n = length(b);                %线性方程组的阶数
if nargin < 5
    maxiter = 500;            %默认最大迭代次数
end
if nargin < 4
    ep = 1e-8;                %默认精度
end
if nargin < 3 || isempty(x0)
    x0 = zeros(n,1);           %默认迭代初值
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上为初始化过程相当于：
% maxiter = 500;
% ep = 1e-8;
% x0 = zeros(n,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x0 = x0(:);                    %将初始迭代解向量变成列向量
x = zeros(n,1);                %预置零向量存储迭代解
iter = 1 ;                      %迭代次数初始值
exitflag = 1;                  %迭代收敛标志： 1-收敛；0-发散
xs(iter,:) = x0;               %存储每次迭代的值
err(iter,1) = nan;             %初始误差置为nan
while exitflag
    for k=1:n                  %Jacobi迭代
        x(k) = (b(k)-A(k,[1:k-1,k+1:n])*x0([1:k-1,k+1:n]))/A(k,k);
    end
    xs(iter+1,:) = x;          %将迭代解存储在xs矩阵中
    err(iter+1,1) = norm(x-x0);%计算误差并存储
    if err(iter+1,1)<=ep       %前后迭代值误差在允许范围内
        break;                 %退出循环
    end
    x0 = x;                    %更新迭代初始值
    iter = iter + 1;           %迭代次数累加
    if iter > maxiter  %若迭代次数大于最大迭代次数，则认为迭代发散，即根不可靠
        exitflag = 0;          %置exitflag值为0，即发散
    end
end
[varargout{1:5}] = deal(x,...  %第1个输出参数为线性方程组的解
    exitflag,...               %第2个输出参数为迭代收敛的标志
    iter,...                   %第3个输出参数为迭代次数
    xs,...                     %第4个输出参数为迭代序列
    err(2:end));               %第5个输出参数为近似相对误差

Gauss-Seidel迭代法

$$\begin{gathered} B_{GS}=-(D+L{)}^{-1}U \\ {F}_{GS}=(D+L{)}^{-1}b \end{gathered}$$

function varargout = seidelc(A,b,x0,ep,maxiter)
%SEIDEL 分量形式Gauss-Seidel迭代法求解线性方程组
n = length(b);                %线性方程组的阶数
if nargin < 5
    maxiter = 500;            %默认最大迭代次数
end
if nargin < 4
    ep = 1e-8;                %默认精度
end
if nargin < 3 || isempty(x0)
    x0 = zeros(n,1);           %默认迭代初值
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上为初始化过程相当于：
% maxiter = 500;
% ep = 1e-8;
% x0 = zeros(n,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x0 = x0(:);                    %将初始迭代解向量变成列向量
x = zeros(n,1);                %预置零向量存储迭代解
iter = 1 ;                      %迭代次数初始值
exitflag = 1;                  %迭代收敛标志： 1-收敛；0-发散
xs(iter,:) = x0;               %存储每次迭代的值
err(iter,1) = nan;             %初始误差置为nan
while exitflag
    for k=1:n                  %Gauss-Seidel迭代
        x(k) = (b(k)-A(k,[1:k-1,k+1:n])*x1([1:k-1,k+1:n]))/A(k,k);
        x1(k) = x(k);
    end
    xs(iter+1,:) = x;          %将迭代解存储在xs矩阵中
    err(iter+1,1) = norm(x-x0);%计算误差并存储
    if err(iter+1,1)<=ep       %前后迭代值误差在允许范围内
        break                  %退出循环
    end
    x0 = x;                    %更新迭代初始值
    iter = iter + 1;           %迭代次数累加
    if iter > maxiter  %若迭代次数大于最大迭代次数，则认为迭代发散，即根不可靠
        exitflag = 0;          %置exitflag值为0，即发散
    end
end
[varargout{1:5}] = deal(x,...  %第1个输出参数为线性方程组的解
    exitflag,...               %第2个输出参数为迭代收敛的标志
    iter,...                   %第3个输出参数为迭代次数
    xs,...                     %第4个输出参数为迭代序列
    err(2:end));               %第5个输出参数为近似相对误差

SOR迭代法

$$\begin{gathered} B_{SOR}=(D+\omega L{)}^{-1}[(1-\omega )D-\omega U] \\ {F}_{SOR}=\omega (D+\omega L{)}^{-1}b \end{gathered}$$

function varargout = sorc(A,b,w,x0,ep,maxiter)
%SORC 分量形式SOR法求解线性方程组
n = length(b);                %线性方程组的阶数
if nargin < 5
    maxiter = 500;            %默认最大迭代次数
end
if nargin < 4
    ep = 1e-8;                %默认精度
end
if nargin < 3 || isempty(x0)
    x0 = zeros(n,1);           %默认迭代初值
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上为初始化过程相当于：
% maxiter = 500;
% ep = 1e-8;
% x0 = zeros(n,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x0 = x0(:);                    %将初始迭代解向量变成列向量
x = zeros(n,1);                %预置零向量存储迭代解
iter = 1 ;                      %迭代次数初始值
exitflag = 1;                  %迭代收敛标志： 1-收敛；0-发散
xs(iter,:) = x0;               %存储每次迭代的值
err(iter,1) = nan;             %初始误差置为nan
while exitflag
    x1 = x0(:);
    for k=1:n                  %Gauss-Seidel迭代
        x(k) = (1-w)*x0(k)+...
            w*(b(k)-A(k,[1:k-1,k+1:n])*x1([1:k-1,k+1]))/A(k,k);
        x1(k) = x(k);
    end
    xs(iter+1,:) = x;          %将迭代解存储在xs矩阵中
    err(iter+1,1) = norm(x-x0);%计算误差并存储
    if err(iter+1,1)<=ep       %前后迭代值误差在允许范围内
        break                  %退出循环
    end
    x0 = x;                    %更新迭代初始值
    iter = iter + 1;           %迭代次数累加
    if iter > maxiter  %若迭代次数大于最大迭代次数，则认为迭代发散，即根不可靠
        exitflag = 0;          %置exitflag值为0，即发散
    end
end
[varargout{1:5}] = deal(x,...  %第1个输出参数为线性方程组的解
    exitflag,...               %第2个输出参数为迭代收敛的标志
    iter,...                   %第3个输出参数为迭代次数
    xs,...                     %第4个输出参数为迭代序列
    err(2:end));               %第5个输出参数为近似相对误差