迭代求解线性方程组的解

发布时间：2024-12-10 21:35

了解并应用问题解决的迭代方法 #生活技巧# #领导力技巧# #问题解决能力#

码字不易，觉得好，点个赞或收藏下吧！

数值分析编程作业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 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990

测试代码：

%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); 1234567891011121314151617181920212223242526

实验结果

分别用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迭代而言，ω的选取可以通过简单代入迭代尝试的方法来确定。