我试图在MATLAB中实现Jacobi迭代,但无法使其收敛我已经在网上和其他地方寻找工作代码进行比较,但无法找到任何类似于我的代码,仍然有效以下是我所拥有的:
function x = Jacobi(A,b,tol,maxiter)
n = size(A,1);
xp = zeros(n,1);
x = zeros(n,1);
k=0; % number of steps
while(k<=maxiter)
k=k+1;
for i=1:n
xp(i) = 1/A(i,i)*(b(i) - A(i,1:i-1)*x(1:i-1) - A(i,i+1:n)*x(i+1:n));
end
err = norm(A*xp-b);
if(err<tol)
x=xp;
break;
end
x=xp;
end
不管我用什么A和b都会爆炸这可能是我忽略的一个小错误,但如果有人能解释出什么是错误的,我会非常感激,因为这应该是正确的,但在实践中不是这样。
最佳答案
你的代码是正确的它看起来不起作用的原因是,您指定的系统在使用Jacobi迭代时可能不会收敛。
具体来说(多亏了@Saraubh),如果矩阵A
严格地diagonally dominant,这个方法就会收敛换言之,对于矩阵中的每一行i
,所有列j
在i
处的绝对和(在i
处没有对角线系数)必须小于对角线本身换句话说:
然而,有些系统会与Jacobi收敛,即使这个条件不满足,但是在尝试将Jacobi用于您的系统之前,您应该将其作为一般规则使用如果你用高斯赛德尔,它实际上更稳定唯一的区别是,您正在使用x
的解决方案,并在向下执行行时将其输入到其他变量中要制作这个Gauss-Seidel,只需在for
循环中更改一个字符更改为:
xp(i) = 1/A(i,i)*(b(i) - A(i,1:i-1)*x(1:i-1) - A(i,i+1:n)*x(i+1:n));
对此:
xp(i) = 1/A(i,i)*(b(i) - A(i,1:i-1)*xp(1:i-1) - A(i,i+1:n)*x(i+1:n));
**HERE**
下面是两个例子:
这里我们指定一个系统,它不是由Jacobi收敛的,但是有一个解这个系统不是对角占优的。
这里我们指定一个系统,它确实是由Jacobi收敛的具体来说,这个系统是对角占优的。
示例1
A = [1 2 2 3; -1 4 2 7; 3 1 6 0; 1 0 3 4];
b = [0;1;-1;2];
x = Jacobi(A, b, 0.001, 40)
xtrue = A \ b
x =
1.0e+09 *
4.1567
0.8382
1.2380
1.0983
xtrue =
-0.1979
-0.7187
0.0521
0.5104
现在,如果我用高斯-赛德尔解,这就是我得到的:
x =
-0.1988
-0.7190
0.0526
0.5103
哇哦它收敛于Gauss-Seidel而不是Jacobi,即使系统不是对角占优的,我可能对此有一个解释,稍后将提供。
例2
A = [10 -1 2 0; -1 -11 -1 3; 2 -1 10 -1; 0 3 -1 8];
b = [6;25;-11;15];
x = Jacobi(A, b, 0.001, 40);
xtrue = A \ b
x =
0.6729
-1.5936
-1.1612
2.3275
xtrue =
0.6729
-1.5936
-1.1612
2.3274
这就是我从高斯·赛德尔那里得到的:
x =
0.6729
-1.5936
-1.1612
2.3274
这两种情况都肯定会收敛,而且系统是对角占优的。
因此,您的代码没有任何问题您只是指定了一个无法使用Jacobi解决的系统最好使用高斯-赛德尔迭代法来解决这类问题原因是您正在立即使用当前迭代的信息并将其传播到其他变量雅各比没有这样做,这就是为什么它分化得更快对于Jacobi,可以看到示例1未能收敛,而示例2收敛Gauss-Seidel同时收敛事实上,当它们都收敛时,它们非常接近真实的解。
同样,你需要确保你的系统是对角占优的,这样你就能保证收敛不执行此规则好。。。你会冒这个险,因为它可能会,也可能不会收敛。
祝你好运!