解线性方程组——(Gauss-Seidel)高斯-赛德尔迭代法 | 北太天元

最新推荐文章于 2024-05-17 14:15:08 发布

清水折木

最新推荐文章于 2024-05-17 14:15:08 发布

阅读量811

点赞数 9

分类专栏：数值计算方法 | 北太天元文章标签：线性代数数值计算 matlab 算法解线性方程组矩阵

本文链接：https://blog.csdn.net/math_boy_w/article/details/138171722

版权

数值计算方法 | 北太天元专栏收录该内容

16 篇文章 2 订阅

订阅专栏

一、Gauss-Seidel迭代法

$n = 3$ 时
$A=\begin{pmatrix} a_{11} & a_{12} &a_{13}\\ a_{21} & a_{22} &a_{23}\\ a_{31} & a_{32} &a_{33}\\ \end{pmatrix} ,\quad b=\begin{pmatrix} b_1\\b_2\\ b_3 \end{pmatrix},$

先看一下Jacobi公式的特点

Jacobi公式为
$\begin{gathered} x_1^{(k+1)} =\frac{b_{1}-a_{12}x_{2}^{(k)}-a_{13}x_{3}^{(k)}}{a_{11}} \\ x_2^{(k+1)} =\frac{b_{2}-a_{21}x_{1}^{(k)}-a_{23}x_{3}^{(k)}}{a_{22}} \\ x_3^{(k+1)} =\frac{b_{3}-a_{31}x_{1}^{(k)}-a_{32}x_{2}^{(k)}}{a_{33}} \end{gathered}$

可以看出Jacobi方法每次迭代使用的是上一步的结果，每行全部独立计算完后才进入下一轮迭代。

而实际上计算 $x_2$ 时, 本次迭代产生的 $x_1$ 已经更新, 使用 $x_3$ 时， $x_1,x_2$ 已经更新。由此想法（尽可能使用最新产生的结果），得到Gauss-Seidel方法

Gauss-Seidel公式为
$\begin{aligned} x_1^{(k+1)} &=\frac{b_{1}-a_{12}x_{2}^{(k)}-a_{13}x_{3}^{(k)}}{a_{11}} \\ x_2^{(k+1)} &=\frac{b_{2}-a_{21}x_{1}^{(k+1)}-a_{23}x_{3}^{(k)}}{a_{22}} \\ x_3^{(k+1)} &=\frac{b_{3}-a_{31}x_{1}^{(k+1)}-a_{32}x_{2}^{(k+1)}}{a_{33}} \end{aligned}$

或等价的，将 $A$ 分解为 $A = D - L - U$ ，其中
$D=diag(a_{11},a_{22},a_{33}),$ $L=-\begin{bmatrix} 0 & 0 &0\\ a_{21} & 0&0\\ a_{31} & a_{32} &0\\ \end{bmatrix},\quad U=-\begin{bmatrix} 0 & a_{12} &a_{13}\\ 0 & 0&a_{23}\\ 0 & 0 &0\\ \end{bmatrix}.$ $\begin{aligned} (D-L-U)x &= b\\ Dx &= b+(L+U)x\\ Dx^{(k+1)}&=b+Lx^{(k+1)}+Ux^{(k)}\\ (D-L)x^{(k+1)}&=b+Ux^{(k)}\\ x^{(k+1)}&=(D-L)^{-1}(b+Ux^{(k)}) \end{aligned}$

其中 $x^{(k+1)}=D^{-1}(b+Lx^{(k+1)}+Ux^{(k)})$
写成矩阵的形式
$\begin{bmatrix} x_1^{(k+1)}\\x_2^{(k+1)}\\ x_3^{(k+1)} \end{bmatrix} = \begin{bmatrix} \frac{1}{a_{11}} & &\\ & \frac{1}{a_{22}} &\\ & &\frac{1}{a_{33}}\\ \end{bmatrix}\left( \begin{bmatrix} b_1\\b_2\\ b_3 \end{bmatrix}-\begin{bmatrix} & &\\ a_{21} & &\\ a_{31} & a_{32} &\ \\ \end{bmatrix}\begin{bmatrix} x_1^{(k+1)}\\x_2^{(k+1)}\\ x_3^{(k+1)} \end{bmatrix} -\begin{bmatrix} \ & a_{12} &a_{13}\\ & &a_{23}\\ & &\\ \end{bmatrix}\begin{bmatrix} x_1^{(k)}\\x_2^{(k)}\\x_3^{(k)}\end{bmatrix}\right)$
下面是其一般形式下的算法

二、算法

$\heartsuit$ Gauss-Seidel迭代法

主要思路

输入 : $A,b,x^{(0)}$
输出 : $x$
$\begin{aligned} x^{(0)} &=\text{ initial vector } \\ x^{(k+1)}&=(D-L)^{-1}(b+Ux^{(k)}) \end{aligned}$

添加一些限制

容许误差 e_tol,
最大迭代步 $N$ .

当残差 <e_tol 或迭代步数 $\geq N$ 时，都会停止迭代，输出结果

实现步骤

步骤 $1$ : $k = 0$ , $x=x^{(0)}$ ；
步骤 $2$ : 计算残差 $r=\|b-Ax\|$ ，
- 如果残差 $r$ >e_tol 且 $k < N$ ，转步骤 $3$ ;
- 否则, 转步骤 $5$ ；
步骤 $3$ : 更新解向量 $x=(D-L)^{-1}(b+Ux^{(0)})$
步骤 $4$ ： $x 0 = x$ , $k = k + 1$ ，转步骤 $2$ ;
步骤 $5$ : 输出 $x$ .

在这里插入图片描述

三、北太天元源程序

Gauss-Seidel

function [x,k,r] = myGS(A,b,x0,e_tol,N)
% Gauss-Seidel迭代法解线性方程组
% Input: A, b(列向量), x0(初始值)
%        e_tol: error tolerant  
%        N: 限制迭代次数小于 N 次
% Output: x , k(迭代次数),r:残差
%   Version:            1.0
%   last modified:      01/29/2024
    n = length(b); k = 0; 
    x=zeros(n,N); % 记录每一次迭代的结果，方便后续作误差分析
    x(:,1)=x0; 
    L = -tril(A,-1); U = -triu(A,1); D = diag(diag(A));
    r = norm(b - A*x,2);
    while r > e_tol && k < N
        x(:,k+2) = inv(D-L)*(b+U*x(:,k+1)); % 不同之处
        r = norm(b - A*x(:,k+2),2); % 残差
        k = k+1;
    end
    x = x(:,2:k+1); % x取迭代时的结果
    
    if k>N
        fprintf('迭代超出最大迭代次数');
    else
        fprintf('迭代次数=%i\n',k);
    end
end

保存为 myGS.m 文件

四、数值算例

（此例子与 Jaboci迭代法文章中的例子相同，可以对比着来看）

$A x = b$ ,求 $x$
$\begin{pmatrix} 10 & -1 & 2 & 0 \\ -1 & 11 & -1 & 3 \\ 2 & -1 & 10 & -1 \\ 0 & 3 & -1 & 8 \\ \end{pmatrix}\quad b = \begin{pmatrix} 6 \\ 25 \\ -11 \\ 15 \\ \end{pmatrix}$

用Gauss列主元消去法,得

x = 
   1.000000000000000
   2.000000000000000
  -1.000000000000000
   1.000000000000000

下面我们用Gauss-Seidel 迭代法进行求解

%% Gauss-Seidel test
% time :      4/24/2024

%% example 1
clc;clear all,format long;
N = 100; e_tol = 1e-8;
A=[10 -1 2 0; -1 11 -1 3; 2 -1 10 -1; 0 3 -1 8];
b=[6; 25; -11; 15];
x0=[0; 0; 0; 0];
    [x11,k1] = myGS(A,b,x0,e_tol,N)
    [x12,k2] = myJacobi(A,b,x0,e_tol,N)
% 作图查看误差变化
    x_exact=[1;2;-1;1]; %真解
    n = length(b);
    error=zeros(n,k1);% 每个分量的误差
        error = abs(x_exact - x11) 
    res =zeros(1,k1); % 残差
    res1 = res;
    res2 = res;
    for i=1:1:k1
        res1(i) = norm(b-A*x11(:,i),2);
        
    end 
    for i=1:1:k2
        res2(i) = norm(b-A*x12(:,i),2);
    end
    % 数值解
        figure(1);
        plot(1:k1,x11(1,:),'-*r',1:k1,x11(2,:),'-og', 1:k1,x11(3,:),'-+b',1:k1,x11(4,:),'-dk');
        legend('x_1','x_2','x_3','x_4');
        title('G-S下每个数值解的变化')
    % 残差变化
        figure(2);
        plot(1:k1,res1,'-*r');
        hold on
        plot(1:k2,res2,'-*b');
        legend('G-S','Jacobi');
        title('残差变化')
    % 误差
        figure(3);
        plot(1:k1,error(1,:),'-*r',1:k1,error(2,:),'-og', 1:k1,error(3,:),'-+b',1:k1,error(4,:),'-dk');
        legend('x_1','x_2','x_3','x_4');
        title('G-S下每个数值解的误差变化')