MATLAB 中的高斯-赛德尔法
Gauss-Seidel Method in MATLAB
我正在尝试在 MATLAB 中实现 Gauss-Seidel 方法。但是我的代码中有两个主要错误,我无法修复它们:
我的代码在小矩阵上收敛得很好,但在大矩阵上从不收敛。
代码进行了冗余迭代。如何防止重复迭代?
Gauss-Seidel Method on wikipedia.
N=5;
A=rand(N,N);
b=rand(N,1);
x = zeros(N,1);
sum = 0;
xold = x;
tic
for n_iter=1:1000
for i = 1:N
for j = 1:N
if (j ~= i)
sum = sum + (A(i,j)/A(i,i)) * xold(j);
else
continue;
end
end
x(i) = -sum + b(i)/A(i,i);
sum = 0;
end
if(abs(x(i)-xold(j))<0.001)
break;
end
xold = x;
end
gs_time=toc;
prompt1='Gauss-Seidel Method Time';
prompt2='x Matrix';
disp(prompt2);
disp(x);
disp(prompt1);
disp(gs_time);
%Gauss-seidal method for three equations
clc;
x1=0;
x2=0;
x3=0;
m=input('Enter number of iteration');
for i=1:1:m
x1(i+1)=(-0.01-0.52*x2(i)-x3(i))/0.3
x2(i+1)=0.67-1.9*x3(i)-0.5*x1(i+1)
x3(i+1)=(0.44-0.1*x1(i+1)-0.3*x2(i+1))/0.5
er1=abs((x1(i+1)-x1(i))/x1(i+1))*100
er2=abs((x2(i+1)-x2(i))/x2(i+1))*100
er3=abs((x3(i+1)-x3(i))/x3(i+1))*100
if er1<=0.01
er2<=0.01
er3<=0.01
break;
end
end
首先,一般性。 Gauß-Seidel 和 Jacobi 方法 仅适用于 对角占优矩阵 ,不适用于一般随机矩阵。因此,要获得正确的测试示例,您需要实际建设性地确保该条件,例如通过
A = rand(N,N)+N*eye(N)
或类似。
Else the method will diverge towards infinity in some or all components.
现在谈谈您的实施中的一些其他奇怪之处。什么
if(abs(x(i)-xold(j))<0.001)
是什么意思?请注意,此指令 在循环 之外,其中 i
和 j
是迭代变量,因此索引值可能是 undefined。由于惯性,它们会意外地同时具有值 N
,因此该标准至少有一点意义。
你要测试的是向量作为一个整体的差异范数,因此使用sum(abs(x-xold))/N
或max(abs(x-xold))
。在右侧,您可能希望乘以应用于 x
的相同范数构造,以便测试相对误差,同时考虑问题的规模。
根据给定代码中的说明,您正在实施 Jacobi 迭代,首先计算所有更新,然后推进迭代向量。对于 Gauß-Seidel 变体,您需要就地替换单个组件,以便立即使用新计算的值。
此外,您可以 shorten/simplify 内部循环
xold = x;
for i = 1:N
sum = b(i);
for j = 1:N
if (j ~= i)
sum = sum - A(i,j) * x(j);
end
end
x(i) = sum/A(i,i);
end
err = norm(x-xold)
使用 matlab 的语言特性甚至更短
xold = x
for i = 1:N
J = [1:(i-1) (i+1):N];
x(i) = ( b(i) - A(i,J)*x(J) )/A(i,i);
end
err = norm(x-xold)
我正在尝试在 MATLAB 中实现 Gauss-Seidel 方法。但是我的代码中有两个主要错误,我无法修复它们:
我的代码在小矩阵上收敛得很好,但在大矩阵上从不收敛。
代码进行了冗余迭代。如何防止重复迭代?
Gauss-Seidel Method on wikipedia.
N=5;
A=rand(N,N);
b=rand(N,1);
x = zeros(N,1);
sum = 0;
xold = x;
tic
for n_iter=1:1000
for i = 1:N
for j = 1:N
if (j ~= i)
sum = sum + (A(i,j)/A(i,i)) * xold(j);
else
continue;
end
end
x(i) = -sum + b(i)/A(i,i);
sum = 0;
end
if(abs(x(i)-xold(j))<0.001)
break;
end
xold = x;
end
gs_time=toc;
prompt1='Gauss-Seidel Method Time';
prompt2='x Matrix';
disp(prompt2);
disp(x);
disp(prompt1);
disp(gs_time);
%Gauss-seidal method for three equations
clc;
x1=0;
x2=0;
x3=0;
m=input('Enter number of iteration');
for i=1:1:m
x1(i+1)=(-0.01-0.52*x2(i)-x3(i))/0.3
x2(i+1)=0.67-1.9*x3(i)-0.5*x1(i+1)
x3(i+1)=(0.44-0.1*x1(i+1)-0.3*x2(i+1))/0.5
er1=abs((x1(i+1)-x1(i))/x1(i+1))*100
er2=abs((x2(i+1)-x2(i))/x2(i+1))*100
er3=abs((x3(i+1)-x3(i))/x3(i+1))*100
if er1<=0.01
er2<=0.01
er3<=0.01
break;
end
end
首先,一般性。 Gauß-Seidel 和 Jacobi 方法 仅适用于 对角占优矩阵 ,不适用于一般随机矩阵。因此,要获得正确的测试示例,您需要实际建设性地确保该条件,例如通过
A = rand(N,N)+N*eye(N)
或类似。
Else the method will diverge towards infinity in some or all components.
现在谈谈您的实施中的一些其他奇怪之处。什么
if(abs(x(i)-xold(j))<0.001)
是什么意思?请注意,此指令 在循环 之外,其中 i
和 j
是迭代变量,因此索引值可能是 undefined。由于惯性,它们会意外地同时具有值 N
,因此该标准至少有一点意义。
你要测试的是向量作为一个整体的差异范数,因此使用sum(abs(x-xold))/N
或max(abs(x-xold))
。在右侧,您可能希望乘以应用于 x
的相同范数构造,以便测试相对误差,同时考虑问题的规模。
根据给定代码中的说明,您正在实施 Jacobi 迭代,首先计算所有更新,然后推进迭代向量。对于 Gauß-Seidel 变体,您需要就地替换单个组件,以便立即使用新计算的值。
此外,您可以 shorten/simplify 内部循环
xold = x;
for i = 1:N
sum = b(i);
for j = 1:N
if (j ~= i)
sum = sum - A(i,j) * x(j);
end
end
x(i) = sum/A(i,i);
end
err = norm(x-xold)
使用 matlab 的语言特性甚至更短
xold = x
for i = 1:N
J = [1:(i-1) (i+1):N];
x(i) = ( b(i) - A(i,J)*x(J) )/A(i,i);
end
err = norm(x-xold)