gpt4 book ai didi

matlab - 雅可比迭代没有结束

转载 作者:太空宇宙 更新时间:2023-11-03 20:14:18 25 4
gpt4 key购买 nike

我试图在 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,它都会爆炸。这可能是我忽略的一个小错误,但如果有人能解释问题所在,我将不胜感激,因为这应该是正确的,但在实践中并非如此。

最佳答案

您的代码是正确的。它似乎不起作用的原因是因为您指定的系统在使用雅可比迭代时可能不会收敛。

具体来说(感谢@Saraubh),如果您的矩阵 A严格 diagonally dominant,此方法将收敛.换句话说,对于矩阵中的每一行 ii 行所有 j 列的绝对总和,不包括对角线系数i 必须小于 对角线本身。换句话说:

blah

但是,有些系统会与 Jacobi 收敛,即使不满足此条件,但在尝试将 Jacobi 用于您的系统之前,您应该将此作为一般规则使用。如果你使用 Gauss-Seidel,它实际上更稳定。唯一的区别是,您正在重新使用 x 的解决方案,并在您沿着行前进时将其提供给其他变量。要制作这个高斯赛德尔,您所要做的就是更改 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**

下面是我将向您展示的两个示例:

  1. 我们指定一个系统,该系统不按雅可比收敛,但有一个解决方案。该系统对角占优势。
  2. 我们在这里指定一个通过雅可比收敛的系统。具体来说,这个系统是对角线占优的。

示例#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

现在,如果我使用 Gauss-Seidel 解,这就是我得到的:

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 求解的系统。对于围绕这种求解的迭代方法,最好使用 Gauss-Seidel。原因是因为您立即使用当前迭代中的信息并将其传播到其余变量。雅可比不这样做,这就是它发散更快的原因。对于 Jacobi,您可以看到示例 #1 无法收敛,而示例 #2 可以。 Gauss-Seidel 收敛于两者。事实上,当它们都收敛时,它们非常接近真正的解决方案。

同样,您需要确保您的系统是对角线优势的,这样您才能保证收敛。不执行此规则...好吧...您将承担风险,因为它可能收敛也可能不收敛。

祝你好运!

关于matlab - 雅可比迭代没有结束,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/24730993/

25 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com