解线性方程组——(Jacobi)雅克比迭代法

解线性方程组——(Jacobi)雅克比迭代法 | 北太天元

一、Jacobi迭代法

$n = 3$ , 阶数为 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公式为
$\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迭代可以用于并行计算。

或等价的，将 $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{gathered} (D-L-U)x = b\\ Dx = b+(L+U)x\\ x=D^{-1}(b+(L+U)x)\\ \end{gathered}$
得Jacobi公式 $x^{(k+1)}=D^{-1}(b+(L+U)x^{(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_{12} &a_{13}\\ a_{21} & &a_{23}\\ a_{31} & a_{32} &\\ \end{bmatrix}\begin{bmatrix} x_1^{(k)}\\x_2^{(k)}\\ x_3^{(k)} \end{bmatrix}\right)$
下面是其一般形式下的算法

二、算法

💖 Jacobi迭代法

主要思路

输入: $A,b,x^{(0)}$ ，输出 $x$
$\begin{aligned} x^{(0)} &= \text{ initial vector } \\ x^{(k+1)} &= D^{-1}(b+(L+U)x^{(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^{-1}(b+(L+U)x^{(0)})$
步骤 4： $x 0 = x$ ， $k = k + 1$ ，转步骤 2；
步骤 5: 输出 $x$ 。

在这里插入图片描述

三、北太天元源程序

Jacobi迭代

function [x,k,r] = myJacobi(A,b,x0,e_tol,N)
% Jacobi迭代法解线性方程组
% Input: A, b(列向量), x0(初始值)
%        e_tol: error tolerant  
%        N: 限制迭代次数小于 N 次
% Output: x , k(迭代次数),r:残差
%   Version:            1.0
%   last modified:      01/27/2024n = 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 < Nx(:,k+2) = inv(D)*(b+(L+U)*x(:,k+1));r = norm(b - A*x(:,k+2),2); % 残差k = k+1;endx = x(:,2:k+1); % x取迭代时的结果if k>Nfprintf('迭代超出最大迭代次数');elsefprintf('迭代次数=%i\n',k);end
end

保存为 myJacobi.m文件

四、数值算例

例1

$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.0000000000000002.000000000000000-1.0000000000000001.000000000000000

下面我们用Jacobi迭代法进行求解

设 N = 100 , e_tol = 1e-8, x0=[0; 0; 0; 0]

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];
t1 =tic;[x11,k1] = myJacobi(A,b,x0,e_tol,N)
toc(t1);
t2 = tic;[x12] = gsem_column(A,b)
toc(t2);
% 作图查看误差变化x_exact=[1;2;-1;1]; %真解n = length(b);error=zeros(n,k1);% 每个分量的误差error = abs(x_exact - x11) res =zeros(1,k1); % 残差for i=1:1:k1res(i) = norm(b-A*x11(:,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('每个数值解的变化')% 残差变化figure(2);plot(1:k1,res,'-*r');legend('残差');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('误差变化')

得

迭代次数=26

在这里插入图片描述

例2

$\begin{bmatrix}3&-1&0&0&0&\frac{1}{2}\\-1&3&-1&0&\frac{1}{2}&0\\0&-1&3&-1&0&0\\0&0&-1&3&-1&0\\0&\dfrac{1}{2}&0&-1&3&-1\\\frac{1}{2}&0&0&0&-1&3\end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\\x_4\\x_4\\x_5\\x_6\end{bmatrix}=\begin{bmatrix}\frac{5}{2}\\ \frac{3}{2}\\1\\1\\\frac{3}{2}\\\frac{5}{2}\end{bmatrix}.$
使用Gauss消去法与Jacobi迭代法分别对其进行求解，并观察有何差异。
试着得到在更高阶数下的结果。

$x$ 的真解为 $1,1,1,1,1,1]^T$

稀疏矩阵的构造：

function [A,b,x_sp] = setup_Sparse1(n)% 定义一个 n 阶的稀疏矩阵， 真解全为1% Input: n% Output: A,b%%   Version:            1.0%   last modified:      09/28/2023A = zeros(n,n);b = ones(n,1) * 1.5;b([1,n],1) = 2.5; b([n/2,n/2 + 1],1) = 1;x_sp = ones(n,1); % 真解for i = 1:1:nA(i,n+1-i) = 1/2;A(i,i) = 3;endfor i =1:1:n-1A(i,i+1)= -1;A(i+1,i)= -1;end
end

保存为 setup_Sparse1.m 文件

实现代码

%% 测试 消去法 与 迭代法 在处理稀疏矩阵问题上的差距
clc;clear all,format long;
N = 100; e_tol = 1e-8;
n = 6; %  n 为 6 50 100 500 1000
[A,b,x]=setup_Sparse1(n);
x0 = zeros(n,1); t1 =tic;
[x11,k1,r1] = myJacobi(A,b,x0,e_tol,N);
toc(t1);t2 = tic;
[x12] = gsem_column(A,b);
toc(t2);
r2 = norm(b-A*x12);