问题描述：
matlab 关于break的问题
写了一个关于Gauss-Seidel iteration的东西,题目说要用到break,看说明只知道是跳出当前for loop,于是这样写了.A是4*4的矩阵,b是4*1的矩阵.A*x =b:
error = 10.^(-8);
k = 0;
xk = zeros(4,1);
r = ones(4,1);
sum1 = 0;
sum2 = 0;
xadd = zeros(4,1);
while (k < 1001) && norm(r,2) >= (error * norm(b,2))
for i = 1:4
for j1 = 1:(i-1)
if i==1,break,end
sum1 = sum1 + A(i,j1)*xadd(j1,1);
end
for j2 = (i+1):4
if i==4,break,end
sum2 = sum2 + A(i,j2)*xk(j2,1);
end
xadd(i,1) = 1\A(i,i)*(b(i,1)-sum1-sum2);
end
xk = xadd;
r = b - A*xadd;
k = k + 1;
end
在workspace里j1是3（为什么不是4?）,j2是空白（【】）,其他数据看要么是NAN要么就很奇怪.是不是break用错了?还是有其他问题?
写了一个关于Gauss-Seidel iteration的东西,题目说要用到break,看说明只知道是跳出当前for loop,于是这样写了.A是4*4的矩阵,b是4*1的矩阵.A*x =b:
error = 10.^(-8);
k = 0;
xk = zeros(4,1);
r = ones(4,1);
sum1 = 0;
sum2 = 0;
xadd = zeros(4,1);
while (k < 1001) && norm(r,2) >= (error * norm(b,2))
for i = 1:4
for j1 = 1:(i-1)
if i==1,break,end
sum1 = sum1 + A(i,j1)*xadd(j1,1);
end
for j2 = (i+1):4
if i==4,break,end
sum2 = sum2 + A(i,j2)*xk(j2,1);
end
xadd(i,1) = 1\A(i,i)*(b(i,1)-sum1-sum2);
end
xk = xadd;
r = b - A*xadd;
k = k + 1;
end
在workspace里j1是3（为什么不是4?）,j2是空白（【】）,其他数据看要么是NAN要么就很奇怪.是不是break用错了?还是有其他问题?
问题解答：
我来补答
1、你这完全不是用MATLAB写Gauss-Seidel迭代程序的套路啊,给你两个别人写的程序供参考：
gaussSeidel.m
function [v,sN,vChain]=gaussSeidel(A,b,x0,errorBound,maxSp)
%Gauss-Seidel迭代法求解线性方程组
%A-系数矩阵 b-右端向量 x0-初始迭代点 errorBound-近似精度 maxSp-最大迭代次数
%v-近似解 sN-迭代次数 vChain-迭代过程的所有值
step=0;
error=inf;
s=size(A);
D=zeros(s(1));
vChain=zeros(15,3);%最多能记录15次迭代次数
k=1;
fx0=x0;
for i=1:s(1)
D(i,i)=A(i,i);
end;
L=-tril(A,-1);
U=-triu(A,1);
while error>=errorBound & step<maxSp
x0=inv(D)*(L+U)*x0+inv(D)*b;
vChain(k,:)=x0';
k=k+1;
error=norm(x0-fx0);
fx0=x0;
step=step+1;
end
v=x0;
sN=step;
gauseidel.m
function [y,n]=gauseidel(A,b,x0,eps)
% Gauss-Serdel迭代法的MATLAB函数文件gauseidel.m如下：
if nargin==3
eps=1.0e-6;
elseif nargin<3
error
return
end
D=diag(diag(A)); %求A的对角矩阵
L=-tril(A,-1); %求A的下三角阵
U=-triu(A,1); %求A的上三角阵
G=(D-L)\U;
f=(D-L)\b;
y=G*x0+f;
n=1; %迭代次数
while norm(y-x0)>=eps
x0=y;
y=G*x0+f;
n=n+1;
end
调用实例：
>> A=[10,-1,0;-1,10,-2;0,-2,10];
>> b=[9,7,6]';
>> x=gauseidel(A,b,[0,0,0]',1.0e-6)
x =
0.99579
0.95789
0.79158
>> x=gaussSeidel(A,b,x0,0.00001,11)
x =
0.99579
0.95789
0.79158
>> x=A\b
x =
0.99579
0.95789
0.79158
其中最后一个x=A\b是求解线性方程组的最正宗的解法.
2、关于你说的为什么最后变量是那些值的问题,楼上的回答应该比较清楚了.
补充一下你追问的内容：那个break是在for j2 = (i+1):4循环体里面,而那个循环如楼上所说,根本就不会执行,也就谈不上什么退出循环了.
3、我再考虑一下,按照你的思路写程序的话,应该怎么改（也就是尽量保留你的程序框架）.
再问： 是这样的，老师给出了写Gauss-Seidel iteration（其实我根本不知道这是啥）的基本步骤，要按照那些步骤写（当然中间可能我有理解错误的地方）。
再答： 这段算法的程序根本没必要搞那么复杂，可参见下面的代码：A=[10,-1,0 1;-1,10,-2 1;0,-2,10 1; 1 0 0 10];
b=[9,7,6,5]';
error = 10.^(-8);
n = length(b);
xk = zeros(n,1);
for k = 0 : 1000
for i = 1 : n
xk(i) = A(i,i) \ ( b(i) - sum(A(i,1:i-1)*xk(1:i-1)) ...
- sum(A(i,i+1:end)*xk(i+1:end)) );
end
r = b - A*xk;
if norm(r,2) <= (error * norm(b,2)), break, end
end
xk 和你原来的代码相比，主要改动有以下几处：1、按题意，把最外层的while循环改成for循环，把结束条件放在循环体里面用break；2、把递推公式中的两处求和用向量运算完成，大大减少了代码量，也更清晰，也就根本不涉及你原来纠结的两个break到底怎么回事了；3、删去很多没有必要的变量如sum1、sum2、xadd之类。4、把原来的4改成n，可适应不同阶次的矩阵A、b。 另外，你原来的代码是错的，比如xadd(i,1) = 1\A(i,i)*(b(i,1)-sum1-sum2);应该是除以A(i,i)，你给写成乘了，可能你没搞清楚左除（\）是什么意思吧。 错误还不止这一处，具体你自己再分析一下吧（sum1、sum2赋初值的位置不对）。
gaussSeidel.m
function [v,sN,vChain]=gaussSeidel(A,b,x0,errorBound,maxSp)
%Gauss-Seidel迭代法求解线性方程组
%A-系数矩阵 b-右端向量 x0-初始迭代点 errorBound-近似精度 maxSp-最大迭代次数
%v-近似解 sN-迭代次数 vChain-迭代过程的所有值
step=0;
error=inf;
s=size(A);
D=zeros(s(1));
vChain=zeros(15,3);%最多能记录15次迭代次数
k=1;
fx0=x0;
for i=1:s(1)
D(i,i)=A(i,i);
end;
L=-tril(A,-1);
U=-triu(A,1);
while error>=errorBound & step<maxSp
x0=inv(D)*(L+U)*x0+inv(D)*b;
vChain(k,:)=x0';
k=k+1;
error=norm(x0-fx0);
fx0=x0;
step=step+1;
end
v=x0;
sN=step;
gauseidel.m
function [y,n]=gauseidel(A,b,x0,eps)
% Gauss-Serdel迭代法的MATLAB函数文件gauseidel.m如下：
if nargin==3
eps=1.0e-6;
elseif nargin<3
error
return
end
D=diag(diag(A)); %求A的对角矩阵
L=-tril(A,-1); %求A的下三角阵
U=-triu(A,1); %求A的上三角阵
G=(D-L)\U;
f=(D-L)\b;
y=G*x0+f;
n=1; %迭代次数
while norm(y-x0)>=eps
x0=y;
y=G*x0+f;
n=n+1;
end
调用实例：
>> A=[10,-1,0;-1,10,-2;0,-2,10];
>> b=[9,7,6]';
>> x=gauseidel(A,b,[0,0,0]',1.0e-6)
x =
0.99579
0.95789
0.79158
>> x=gaussSeidel(A,b,x0,0.00001,11)
x =
0.99579
0.95789
0.79158
>> x=A\b
x =
0.99579
0.95789
0.79158
其中最后一个x=A\b是求解线性方程组的最正宗的解法.
2、关于你说的为什么最后变量是那些值的问题,楼上的回答应该比较清楚了.
补充一下你追问的内容：那个break是在for j2 = (i+1):4循环体里面,而那个循环如楼上所说,根本就不会执行,也就谈不上什么退出循环了.
3、我再考虑一下,按照你的思路写程序的话,应该怎么改（也就是尽量保留你的程序框架）.
再问： 是这样的，老师给出了写Gauss-Seidel iteration（其实我根本不知道这是啥）的基本步骤，要按照那些步骤写（当然中间可能我有理解错误的地方）。
再答： 这段算法的程序根本没必要搞那么复杂，可参见下面的代码：A=[10,-1,0 1;-1,10,-2 1;0,-2,10 1; 1 0 0 10];
b=[9,7,6,5]';
error = 10.^(-8);
n = length(b);
xk = zeros(n,1);
for k = 0 : 1000
for i = 1 : n
xk(i) = A(i,i) \ ( b(i) - sum(A(i,1:i-1)*xk(1:i-1)) ...
- sum(A(i,i+1:end)*xk(i+1:end)) );
end
r = b - A*xk;
if norm(r,2) <= (error * norm(b,2)), break, end
end
xk 和你原来的代码相比，主要改动有以下几处：1、按题意，把最外层的while循环改成for循环，把结束条件放在循环体里面用break；2、把递推公式中的两处求和用向量运算完成，大大减少了代码量，也更清晰，也就根本不涉及你原来纠结的两个break到底怎么回事了；3、删去很多没有必要的变量如sum1、sum2、xadd之类。4、把原来的4改成n，可适应不同阶次的矩阵A、b。 另外，你原来的代码是错的，比如xadd(i,1) = 1\A(i,i)*(b(i,1)-sum1-sum2);应该是除以A(i,i)，你给写成乘了，可能你没搞清楚左除（\）是什么意思吧。 错误还不止这一处，具体你自己再分析一下吧（sum1、sum2赋初值的位置不对）。
展开全文阅读