数值分析算法 MATLAB 实践 线性方程组迭代法
Gauss-Seidel迭代法
%% 求线性方程组的Gauss-Seidel迭代法，调用格式为
% [x, k] = guaseidel(A,b,x0,eps,it_max)
% 其中, A 为线性方程组的系数矩阵，b 为常数项，eps 为精度要求，默认为1e-5,
% it_max 为最大迭代次数，默认为100
% x 为线性方程组的解，k迭代次数
x0=[0,0,0]';%[x1;x2;x3]列向量
it_max = 1000;eps=1e-6;
[x3,k3_cnt] = GuaSeidelFunc(A,b,x0,eps,it_max);
disp('迭代次数:k3_cnt=');
disp(k3_cnt)
disp(['方程组的解：x3 = ']);
disp(x3)
%% 求线性方程组的Gauss-Seidel迭代法，调用格式为
% [x,k]=GuaSeideFunmethod(A,b,x0,ep,N)
% 其中, A 为线性方程组的系数矩阵，b 为常数项，ep为精度要求，默认为1e-6,
% N 为最大迭代次数，默认为500
% x 为线性方程组的解，k迭代次数
x0=[0,0,0]';%[x1;x2;x3]列向量
it_max = 1000;eps=1e-6;
[x4,k4_cnt]=GuaSeideFunmethod(A,b,x0,eps,it_max);
disp('迭代次数:k4_cnt=');
disp(k4_cnt)
disp(['方程组的解：x4 = ']);
disp(x4)
%% 求线性方程组的GuessSeidel迭代法，调用格式为
% function [x,k] = GuessSeidel(A,b,eps,it_max)
% 其中, A 为线性方程组的系数矩阵，b 为常数项，eps为精度要求，默认为1e-6,
% it_max 为最大迭代次数，默认为500
% x 为线性方程组的解，k迭代次数
it_max = 1000;eps=1e-6;
[x5,k5_cnt] = GuessSeidel(A,b,eps,it_max);
disp('迭代次数:k5_cnt=');
disp(k5_cnt)
disp(['方程组的解：x5 = ']);
disp(x5)
function [x,k] = GuessSeidel(A,b,eps,it_max)
%%GuessSeidel:高斯-赛德尔方法求解线性方程组
%高斯-赛德尔方法是一种迭代法，首先猜测各个xi的初始值（一个简单的方法是设各个xi为0）
%将这些初始值带入到第一个方程解出x1,然后更新x1，将xi带入第二个方程x2，更新x2
%依次迭代，直至数值解非常接近真实值为止
%判断条件：对任意的i，有ea(i)=abs((x(i)-xold(i))/x(i))<=es
%即：max(ea)<=es
%%输入
%A=系数矩阵
%b=右侧矩阵
%es=终止准则（default = 0.00001%）
%maxit=最大迭代次数（default = 500）
%输出：
%x=解向量
%%代码实现
%思路：解向量可以简单的表示为x=d-C*x
%其中di=b_i/a_ii,C的对角线元素为0。
if nargin<2
error('至少输入系数矩阵和右侧矩阵')
end
if nargin<4||isempty(it_max)
it_max=500;
end
if nargin<3||isempty(eps)
eps=1e-6;
end
[m,n]=size(A);
if m~=n
error('系数矩阵必须为方阵')
end
%求解C
k = 0;
C = A;
for i = 1:n
C(i,i)=0;
x(i) = 0;%顺便求初始x
end
x=x';
for i = 1:n
C(i,1:n) = C(i,1:n)/A(i,i);
end
%求解d
for i=1:n
d(i)=b(i)/A(i,i);
end
%开始迭代
iter=0;
while(1)
xprev=x;%记录上次的x
for i=1:n
x(i)=d(i) - C(i,:)*x;%求解并更新xi
if x(i)~=0
ea(i)=abs((x(i)-xprev(i))/x(i));
end
end
iter = iter+1;
if max(ea)<=eps || iter>=it_max
break
end
k=k+1;
end
end
function [x,k] = GuessSeidel(A,b,eps,it_max)
%%GuessSeidel:高斯-赛德尔方法求解线性方程组
%高斯-赛德尔方法是一种迭代法，首先猜测各个xi的初始值（一个简单的方法是设各个xi为0）
%将这些初始值带入到第一个方程解出x1,然后更新x1，将xi带入第二个方程x2，更新x2
%依次迭代，直至数值解非常接近真实值为止
%判断条件：对任意的i，有ea(i)=abs((x(i)-xold(i))/x(i))<=es
%即：max(ea)<=es
%%输入
%A=系数矩阵
%b=右侧矩阵
%es=终止准则（default = 0.00001%）
%maxit=最大迭代次数（default = 500）
%输出：
%x=解向量
%%代码实现
%思路：解向量可以简单的表示为x=d-C*x
%其中di=b_i/a_ii,C的对角线元素为0。
if nargin<2
error('至少输入系数矩阵和右侧矩阵')
end
if nargin<4||isempty(it_max)
it_max=500;
end
if nargin<3||isempty(eps)
eps=1e-6;
end
[m,n]=size(A);
if m~=n
error('系数矩阵必须为方阵')
end
%求解C
k = 0;
C = A;
for i = 1:n
C(i,i)=0;
x(i) = 0;%顺便求初始x
end
x=x';
for i = 1:n
C(i,1:n) = C(i,1:n)/A(i,i);
end
%求解d
for i=1:n
d(i)=b(i)/A(i,i);
end
%开始迭代
iter=0;
while(1)
xprev=x;%记录上次的x
for i=1:n
x(i)=d(i) - C(i,:)*x;%求解并更新xi
if x(i)~=0
ea(i)=abs((x(i)-xprev(i))/x(i));
end
end
iter = iter+1;
if max(ea)<=eps || iter>=it_max
break
end
k=k+1;
end
end
function [x,k]=GuaSeideFunmethod(A,b,x0,ep,N)
% 求线性方程组的Gauss-Seidel迭代法，调用格式为
% [x,k]=GuaSeideFunmethod(A,b,x0,ep,N)
% 其中, A 为线性方程组的系数矩阵，b 为常数项，ep为精度要求，默认为1e-6,
% N 为最大迭代次数，默认为500
% x 为线性方程组的解，k迭代次数
n=length(b);
if nargin<5
N=500;
end
if nargin<4
ep=1e-6;
end
if nargin<3
x0=zeros(n,1);
k=0;
end
x=zeros(n,1);
k=0;
while k<N
for i=1:n
if i==1
x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1); %开始迭代变量
elseif i==n
x(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n); %最后迭代变量
else %其它迭代变量
x(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);
end
end
if norm(x-x0,inf)<ep
break;
end
x0=x;
%disp('x=');
%disp(x); % 此两行代码可输出中间结果
k=k+1;
end
if k==N
warning('已到达迭代次数上限！');
end
end
