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

大家好，欢迎来到IT知识分享网。

文章对应视频讲解：(Jacobi)雅克比迭代法

一、Jacobi迭代法

$n=3$ , 阶数为 3 时
$$A=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right),b=\left(\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right),$$

或等价的，将 $A$ 分解为 $A=D-L-U$，其中
$$D=diag(a_{11},a_{22},a_{33}),$$$$L=-\left[\begin{array}{ccc}0& 0& 0\\ a_{21}& 0& 0\\ a_{31}& a_{32}& 0\end{array}\right],U=-\left[\begin{array}{ccc}0& a_{12}& a_{13}\\ 0& 0& a_{23}\\ 0& 0& 0\end{array}\right].$$$$\begin{array}{c}(D-L-U)x=b\\ Dx=b+(L+U)x\\ x=D^{-1}(b+(L+U)x)\end{array}$$
得Jacobi公式 $$x^{(k+1)}=D^{-1}(b+(L+U)x^{(k)})$$
写成矩阵的形式
$$\left[\begin{array}{c}{x}_1^{(k+1)}\\ x_2^{(k+1)}\\ x_3^{(k+1)}\end{array}\right]=\left[\begin{array}{ccc}\frac{1}{a_{11}}& & \\ & \frac{1}{a_{22}}& \\ & & \frac{1}{a_{33}}\end{array}\right]\left(\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right]-\left[\begin{array}{ccc}& a_{12}& a_{13}\\ a_{21}& & a_{23}\\ a_{31}& a_{32}& \end{array}\right]\left[\begin{array}{c}{x}_1^{(k)}\\ x_2^{(k)}\\ x_3^{(k)}\end{array}\right]\right)$$
下面是其一般形式下的算法

二、算法

💖 Jacobi迭代法

主要思路

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

添加一些限制

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

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

实现步骤

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

三、北太天元 or matlab 实现

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/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)*(b+(L+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 

保存为 myJacobi.m文件

四、数值算例

例1

$Ax=b$,求 $x$
$$A=\left(\begin{array}{cccc}10& -1& 2& 0\\ -1& 11& -1& 3\\ 2& -1& 10& -1\\ 0& 3& -1& 8\end{array}\right)b=\left(\begin{array}{c}6\\ 25\\ -11\\ 15\end{array}\right)$$

用Gauss列主元消去法,得

x = 1.000000000000000 2.000000000000000 -1.000000000000000 1.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:k1 res(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

$$\left[\begin{array}{cccccc}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& \frac{1}{2}& 0& -1& 3& -1\\ \frac{1}{2}& 0& 0& 0& -1& 3\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_4\\ x_5\\ x_6\end{array}\right]=\left[\begin{array}{c}\frac{5}{2}\\ \frac{3}{2}\\ 1\\ 1\\ \frac{3}{2}\\ \frac{5}{2}\end{array}\right].$$
使用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/2023 A = 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:n A(i,n+1-i) = 1/2; A(i,i) = 3; end for i =1:1:n-1 A(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); 

• Jacobi迭代次数为33次，耗时 0. 秒。 $r=8.5770e-09$
• Gauss列主元耗时 0. 秒。$r=0$

n=50时

• Jacobi迭代次数为84次，耗时 0. 秒。 $r=8.6777e-09$
• Gauss列主元耗时 0. 秒。$r=5.3577e-15$

n=100时

• Jacobi迭代次数为84次，耗时 1.017717 秒。 $r=9.0032e-09$
• Gauss列主元耗时 1. 秒。$r=7.0848e-15$

n=500时

• Jacobi迭代次数为84次，耗时 20. 秒。 $r=9.3455e-09$
• Gauss列主元耗时 21. 秒。$r=8.5065e-15$

n=1000时

• Jacobi迭代次数为84次，耗时 34. 秒。 $r=9.4769e-09$
• Gauss列主元耗时 89. 秒。$r=1.0300e-14$

而Jacobi迭代法，只需把精度控制在所需范围内即算完成任务，能够节约很多时间。