我正在为科学计算做家庭作业,特别是matlab中的迭代方法Gauss-Seidel和SOR,问题是对于一个矩阵给了我意想不到的结果(解决方案不收敛)并且对于另一个矩阵收敛。
这是sor的代码,其中:
- A:系统矩阵 A * x = b
- xini:初始迭代数组
- b:与系统无关的数组 A * x = b
- maxiter:最大迭代次数
- tol:公差;
- 特别是,SOR 方法将接收称为 w 的第六个参数,它对应于松弛参数。
这是sor方法的代码:
function [x2,iter] = sor(A,xIni, b, maxIter, tol,w)
x1 = xIni;
x2 = x1;
iter = 0;
i = 0;
j = 0;
n = size(A, 1);
for iter = 1:maxIter,
for i = 1:n
a = w / A(i,i);
x = 0;
for j = 1:i-1
x = x + (A(i,j) * x2(j));
end
for j = i+1:n
x = x + (A(i,j) * x1(j));
end
x2(i) = (a * (b(i) - x)) + ((1 - w) * x1(i));
end
x1 = x2;
if (norm(b - A * x2) < tol);
break;
end
end
这是 Gauss-seidel 方法的代码:
function [x, iter] = Gauss(A, xIni, b, maxIter, tol)
x = xIni;
xnew = x;
iter = 0;
i = 0;
j = 0;
n = size(A,1);
for iter = 1:maxIter,
for i = 1:n
a = 1 / A(i,i);
x1 = 0;
x2 = 0;
for j = 1:i-1
x1 = x1 + (A(i,j) * xnew(j));
end
for j = i+1:n
x2 = x2 + (A(i,j) * x(j));
end
xnew(i) = a * (b(i) - x1 - x2);
end
x= xnew;
if ((norm(A*xnew-b)) <= tol);
break;
end
end
对于此输入:
A = [1 2 -2; 1 1 1; 2 2 1];
b = [1; 2; 5];
当调用函数 Gauss-Seidel 或 sor 时:
[x, iter] = gauss(A, [0; 0; 0], b, 1000, eps)
[x, iter] = sor(A, [0; 0; 0], b, 1000, eps, 1.5)
高斯的输出是:
x =
1.0e+304 *
1.6024
-1.6030
0.0011
iter =
1000
对于 sor 是:
x =
NaN
NaN
NaN
iter =
1000
但是对于以下系统能够找到解决方案:
A = [ 4 -1 0 -1 0 0;
-1 4 -1 0 -1 0;
0 -1 4 0 0 -1;
-1 0 0 4 -1 0;
0 -1 0 -1 4 -1;
0 0 -1 0 -1 4 ]
b = [1 0 0 0 0 0]'
解决方案:
[x, iter] = sor(A, [0; 0; 0], b, 1000, eps, 1.5)
x =
0.2948
0.0932
0.0282
0.0861
0.0497
0.0195
iter =
52
方法的行为取决于两个矩阵的条件?因为我注意到第二个矩阵比第一个矩阵的条件更好。有什么建议么?