Jacobi(雅可比)迭代原理与matlab,C代码
Jacobi迭代分量形式:
x
i
(
k
+
1
)
=
1
/
a
i
i
(
b
i
−
∑
j
≠
i
i
=
1
n
a
i
j
x
j
(
k
)
)
\ x_i^{(k+1)}=1/a_{ii}(b_i-\sum_{\stackrel{i=1}{ j\neq i} }^na_{ij}x_j^{(k)})
xi(k+1)=1/aii(bi−j=ii=1∑naijxj(k))
矩阵形式:
X
(
k
+
1
)
=
D
−
1
(
L
+
U
)
X
(
k
)
+
D
−
1
b
\ X^{(k+1)}=D^{-1}(L+U)X^{(k)}+D^{-1}b
X(k+1)=D−1(L+U)X(k)+D−1b
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
∗
(
L
+
U
)
,
g
=
D
−
1
b
\ x^{(k+1)}=B_Jx^{(k)}+g,其中B_J=D^{-1}*(L+U),g=D^{-1}b
x(k+1)=BJx(k)+g,其中BJ=D−1∗(L+U),g=D−1b;
使用为无穷范数小于等于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
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