码字不易,觉得好,点个赞或收藏下吧!😊
数值分析编程作业2
对于迭代方法求解线性方程组的解:
- 首先系数矩阵A应当是非奇异方阵,这样能够保证AX=b不是超定方程组,且有唯一的非零解;
- 常用的方法有雅可比迭代法、高斯-赛德尔迭代法,超松弛迭代法和共轭梯度迭代法;
- 对于迭代方法要判定其收敛性,对于Jacobi和G-S迭代
- 系数矩阵A是严格对角占优矩阵的话,两者迭代必收敛;
- 若系数矩阵A对称正定,G-S迭代收敛;
- 然后是X(k+1)=BX(k)+f 中迭代矩阵B的范数(通常是1范数和∞范数)小于1,两迭代方法均收敛(充分条件);
- 求解迭代方程 X(k+1)=BX(k)+f 中迭代矩阵B的谱半径ρ,判断其是否小于1,若是则两迭代法均收敛(充要条件),否则不收敛。
- 对于超松弛迭代(SOR),当权值参数ω为1是,其退化为G-S迭代。
编程实现
对于该问题的MATLAB实现如下:
function xc = linearEqs_iterSolve( x0, tol, A, b, option ,omega)
%实现了Jacobbi迭代和超松弛迭代方法,当参数omega为1时,即为高斯-赛德尔迭代
%x0为线性方程组解的初始值,tol为误差限,A为线性方程组系数矩阵,b为常数向量,
%option为1表示使用Jacobbi迭代,2表示使用超松弛迭代,超松弛因子默认为1,可缺省,也可指定。
if nargin == 5
omega = 1; %如果缺省omega参数,则默认为高斯-赛德尔迭代
end
[D,L,U]=factorizeMatrix(A); %求得系数矩阵的分解矩阵,以便进行迭代矩阵的求解
x_previous = x0;
error = 0; %使用相邻两次迭代解的欧式距离作为误差值
itr_cnt = 1; %记录迭代次数
%disp('NO.ITER X1 X2 X3 ITER_ERROR X5 X6 ITER_ERROR')
disp('NO.ITER X1 X2 X3 X4 X5 X6 ITER_ERROR')
switch(option)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {1} %1代表使用Jacobbi迭代
Bj = D\(L+U);%Jacobbi迭代的迭代矩阵
Fj = D\b;
while 1
x_next = Jacobbi(Bj,Fj,x_previous);
if(itr_cnt ~= 1) %第一次迭代不用进行收敛判断
%error = sum(dist(x_next,x_previous)); %记录每次迭代的误差
error = norm((x_next-x_previous),inf);
%if(error < tol || itr_cnt >100)
if(error < tol ) %当误差小于设定的阈值时,判定为收敛,退出循环
display(itr_cnt, x_next, error);%打印出最后一次迭代得到的不动点信息
break
end
end
x_previous = x_next;%将该次迭代得到的迭代值保存到x_previous中去
display(itr_cnt, x_next, error);%打印迭代的数值信息
itr_cnt = itr_cnt+1;%该次迭代完成,迭代次数记录值加1
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {2}%2代表使用超松弛SOR迭代
L_SOR = inv((D-omega*L))*((1-omega)*D+omega*U);
while 1
x_next = SOR(L_SOR,D,L,b,omega,x_previous);
if(itr_cnt ~= 1) %第一次迭代不用进行收敛判断
%error = sum(dist(x_next,x_previous)); %记录每次迭代的误差
error = norm((x_next-x_previous),inf);
%if(error < tol || itr_cnt >100)
if(error < tol) %当误差小于设定的阈值时,判定为收敛,退出循环
display(itr_cnt, x_next, error);%打印出最后一次迭代得到的不动点信息
break
end
end
x_previous = x_next;%将该次迭代得到的迭代值保存到x_previous中去
display(itr_cnt, x_next, error);%打印迭代的数值信息
itr_cnt = itr_cnt+1;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xc = x_previous;
end
function [dia, low, upper] = factorizeMatrix(A)
dia = diag(diag(A)); %主对角元组成的矩阵
upper = -(triu(A)-dia); %上三角元素的相反数组成的矩阵
low = -(tril(A)-dia); %下三角元素的相反数组成的矩阵
end
function xk = Jacobbi(Bj,Fj,x) %Jacobbi迭代法的求解函数
xk = Bj*x + Fj;
end
function xk = SOR(L_SOR, D, L, b, omega, x) %超松弛迭代法的向量化计算函数
xk = L_SOR*x + omega*inv(D-omega*L)*b;
end
function display(itr, x_next, error) %迭代信息的格式控制输出函数
formatSpec = '%0.13f %0.13f %0.13f %0.13f %0.13f %0.13f';
if(itr < 10)
disp([' ',num2str(itr),' ',sprintf(formatSpec,x_next(1),x_next(2),x_next(3),x_next(4),x_next(5),x_next(6))...
,' ',sprintf('%0.9f',error)]);
else
disp([num2str(itr),' ',sprintf(formatSpec,x_next(1),x_next(2),x_next(3),x_next(4),x_next(5),x_next(6))...
,' ',sprintf('%0.9f',error)]);
end
% formatSpec = '%0.13f %0.13f %0.13f';
% if(itr < 10)
% disp([' ',num2str(itr),' ',sprintf(formatSpec,x_next(1),x_next(2),x_next(3))...
% ,' ',sprintf('%0.9f',error)]);
% else
% disp([num2str(itr),' ',sprintf(formatSpec,x_next(1),x_next(2),x_next(3))...
% ,' ',sprintf('%0.9f',error)]);
% end
end
测试代码:
%Homework 2
A=[3 -1 0 0 0 1/2;
-1 3 -1 0 1/2 0;
0 -1 3 -1 0 0;
0 0 -1 3 -1 0;
0 1/2 0 -1 3 -1;
1/2 0 0 0 -1 3]; %系数矩阵A
b = [5/2 3/2 1 1 3/2 5/2]'; %常数向量b
%x0 = zeros(6,1);
x0_2 = [2 4 3.2 4.1 0.6 1.1]';
tol = 10^(-5);
disp('Jacobi');
xj = linearEqs_iterSolve( x0_2, tol, A, b, 1 );
disp('SOR and omega = 1');
xg_s_1 = linearEqs_iterSolve( x0_2, tol, A, b, 2);
disp('SOR and omega = 1.1');
xg_s_2 = linearEqs_iterSolve( x0_2, tol, A, b, 2, 1.1);
disp('SOR and omega = 1.2');
xg_s_3 = linearEqs_iterSolve( x0_2, tol, A, b, 2, 1.2);
disp('SOR and omega = 1.3');
xg_s_4 = linearEqs_iterSolve( x0_2, tol, A, b, 2, 1.3);
disp('SOR and omega = 1.4');
xg_s_5 = linearEqs_iterSolve( x0_2, tol, A, b, 2, 1.4);
disp('SOR and omega = 1.5');
xg_s_6 = linearEqs_iterSolve( x0_2, tol, A, b, 2, 1.5);
实验结果
分别用Jacobi迭代法,SOR迭代法,其中(omega=1,1.1,1.2,1.3,1.4,1.5)进行迭代求解,可得到迭代结果如下:
可以看出,在选取适当的ω时,SOR迭代的收敛速度要比Jacobi迭代法要快的多,但是当ω选取的值不好的时候,SOR迭代的性能甚至不及Jacobi迭代。同时当ω取1时,SOR迭代就退化为Gauss-Seidel迭代了,G-S迭代可认为是对Jacobi迭代的一种改进,其思想是在迭代时尽量利用最新的信息。此外,对于SOR迭代而言,ω的选取可以通过简单代入迭代尝试的方法来确定。