Jacobi（雅可比）迭代原理与matlab代码_雅各比迭代的理论基础和算法原理

Jacobi（雅可比）迭代原理与matlab，C代码

Jacobi迭代分量形式:
$$x_i^{(k+1)}=1/a_{ii}(b_i-\sum _{\stackrel{i=1}{j\ne i}}^n{a}_{ij}{x}_j^{(k)})$$
矩阵形式：
$$X^{(k+1)}=D^{-1}(L+U)X^{(k)}+D^{-1}b$$

1.理论分析

比较gauss seidel迭代法与J迭代法解方程组差异见图1（以3*3方程组为例），gauss seidel迭代法每次计算使用解x的最新值，而J迭代法使用上一次x的解，如果解是收敛的，那么gauss seidel迭代过程使用的最好的估计值，性能应该比J迭代法好：

2.算法描述

Jacobi迭代法：
根据A=D-L-U将A拆分成对角矩阵D，下三角矩阵L，上三角矩阵U；
循环迭代操作迭代矩阵$$x^{(k+1)}=B_J{x}^{(k)}+g,其中B_J=D^{-1}\ast （L+U），g=D^{-1}b$$；
使用为无穷范数小于等于E-6终止条件；

3.计算程序（matlab）

3.1分量形式用循环不用matlab自带写迭代过程

% Jacobi method A*X=b
 clear
 clc
 n=input('输入问题维度 n:  ');
 A = zeros(n,n);                         %生成矩阵需要的储存单元
 b = zeros(n,1);                         %生成矩阵需要的储存单元
 Xnow = zeros(n);                          %生成解向量需要的储存单元
 Xafter = zeros(n);                          %生成解向量需要的储存单元
 tol = input('输入你题目中的误差: ');     %生成解向量误差x(k+1)-x(k)需要的储存单元
 m = 300;                                %最多迭代300次。超过可能有问题？
 
 A=[4 2 3; 3 -5 2; -2 3 8];
 b=[8 -14 27 ];
 Xnow=[0 0 0];
 
 k = 1;                                 %记录迭代次数
 while k <= m 
   err = 0;
   for i = 1 : n 
      s = 0;
      for j = 1 : n 
        s = s-A(i,j)*Xnow(j);
      end
      s = (b(i)+s)/A(i,i);
      if abs(s) > err 
        err = abs(s);
      end
      Xafter(i) = Xnow(i)+s;
  end

  if err <= tol 
     break;
  else
     k = k+1;
     for i = 1 : n 
       Xnow(i) = Xafter(i);
     end
  end
 end
 
 fprintf('Jacobi方法迭代了 %d 次 :\n', k-1);
 for i = 1 : n 
    fprintf(' %11.8f \n', Xafter(i));
 end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

3.2矩阵迭代用matlab自带的函数求解

clear
clc
 n = input('Enter number of equations, n:  ');          %根据提示输入题目要求的n
 %构造矩阵A
 for m=n
    A=zeros(m,m);
    for m=1:m
        A(m,m)=20;
    end
    for m=2:m
        A(m,m-1)=-8;
        A(m-1,m)=-8;
    end
    for m=3:m
        A(m,m-2)=1;
        A(m-2,m)=1;
    end
 end

x = ones(n,1);
b = zeros(n,1);
D = diag(diag(A));                              %求A的对角矩阵
L = tril(A,-1);                                %求A的下三角矩阵
U = triu(A,1);                                 %求A的上三角矩阵
err = 1E-6;
%Jacobil迭代
I=eye(n);										%生成单位矩阵
BJ=I-D\A;     									%生成迭代矩阵BJ
pJ=vrho(BJ);                                    % 矩阵的谱半径
R_j = -log(pJ);                                  %J迭代渐进收敛速度
g=D\b;
jacobilk = 1;                                   %jacobil迭代次数
xJ = ones(n,1);								%初始迭代向量（1，1，1，1,..）T
it_max = 200;                                   %设定上限，万一无限循环
while jacobilk <= it_max
    xjafter=BJ*xJ+g; 
	xJ=xjafter;
    if norm(xJ,inf)<err
            break;
    end
	jacobilk = jacobilk +1;
end 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42