0.999977856325231 1.999982030452898 0.999976030495353 1.999983560868664 0.999995121966767 1.999995923188666 最后迭代结果:
0.999962179037974 0.999986613620394 0.999990568926511
X =
0.999993879742566 0.999996851108214 1.999995832511947
由以上两种方法的迭代结果可知,SOR法有明显加快收敛的作用.SOR法是广泛使用的求解大型线性方程的迭代法,只要控制好松弛因子,常可获得较好的加速速度,而目前松弛因子的选取尚无有效的方法,只能通过经验或反复的试验才能获得,这也是SOR的一个不足之处.
例3.4 就例3.1中的系数阵A1和例3.3中的系数阵A2:
?10?2?1??4?2?1?????A1???210?1?,A2???24?2?
??1?15???1?23?????讨论Jacobi迭代法和Gauss-Seidel迭代法的收敛性.
解:
先讨论它对Jacobi法的收敛性,再求解它们对的收敛性,求解程序如下:
%计算例3.4
%example3_4.m clear all; clc;
A1=[10 -2 -1;-2 10 -1;-1 -2 5] A2=[4 -2 -1;-2 4 -2;-1 -2 3] %对A1求解
%构造A1的Jocobi迭代矩阵BJ n=max(size(A1)); D=zeros(n); for i=1:n
D(i,i)=A1(i,i); end
BJ_1=-1.0*inv(D)*(A1-D)
46
%利用定理3.2判断A1对Jacobi迭代法的收敛性, %为此求其无穷大范数 BJ1_inf=norm(BJ_1,inf) %对A2求解
%构造A2的Jocobi迭代矩阵BJ n=max(size(A2)); D=zeros(n); for i=1:n
D(i,i)=A2(i,i); end
BJ_2=-1.0*inv(D)*(A2-D)
%利用定理3.2判断A2对Jacobi迭代法的收敛性, %求无穷大范数
BJ2_inf=norm(BJ_2,inf)
%以下判断对Gauss_Seidel迭代法的收敛性 %对A1
D=zeros(n); L=D; U=D;
for i=1:n
for j=1:n if i==j
D(i,i)=A1(i,i); elseif i>j
L(i,j)=A1(i,j); elseif i U(i,j)=A1(i,j); end end end BS1=-1.0*inv(D+L)*U BS1_inf=norm(BS1,inf) %对A2 D=zeros(n); L=D; 47 U=D; for i=1:n for j=1:n if i==j D(i,i)=A2(i,i); elseif i>j L(i,j)=A2(i,j); elseif i U(i,j)=A2(i,j); end end end BS2=-1.0*inv(D+L)*U BS2_inf=norm(BS2,inf) 程序运行结果如下: A1 = 10 -2 -1 -2 10 -1 -1 -2 5 A2 = 4 -2 -1 -2 4 -2 -1 -2 3 BJ_1 = 0 0.2000 0.1000 0.2000 0 0.1000 0.2000 0.4000 0 BJ1_inf = 0.6000 BJ_2 = 0 0.5000 0.2500 0.5000 0 0.5000 0.3333 0.6667 0 BJ2_inf = 1 BS1 = 0 0.2000 0.1000 0 0.0400 0.1200 0 0.0560 0.0680 BS1_inf = 0.3000 BS2 = 0 0.5000 0.2500 48 0 0.2500 0.6250 0 0.3333 0.5000 BS2_inf = 0.8750 Jacobi从运行结果分析可知: BJ(A1)??0.6?1,Bs(A1对1)??0.3?1,因此A迭代法和Gauss?Seidel迭代法都收剑. Bs(A2)??0.875?1,A2对Gauss?Seidel迭代法收敛,而BJ(A2)??1,故A2无法用定理3.2判断其收敛性. 为判断A2对Jacobi的收敛性,改为求其谱半径的方法. >>R_A2=max(abs(eig(BJ_2))) R_A2 = 0.9207 其谱半径?(BJ(A2))?0.9207?1,可见A2对Jacobi迭代法还是收敛的. 例 3.5 就例3.1中系数矩阵A用从矩阵中判断收敛性的方法判断它对迭代法的收敛性,A矩阵为: ?10?2?1???A???210?1? ??1?25???解: 对该题我们可编制一个判断矩阵是否对角占优的函数IsDiagominant(A),程序编制如下： function x=IsDiagominant(A) %判断输入的矩阵A是否具有严格对角优势或A不可约且具有对角优势 %返回值为1,表示A是具有严格对角优势或A不可约且具有对角优势 %若返回值为0则没有上述性质 [m,n]=size(A); if m~=n error('A不是方阵'); return; end %先判断是否为严格对角占优 49 B=abs(A); for i=1:n b=B(i,i); temp=B(i,:); temp(i)=[]; if b>sum(temp) T(i)=1; else T(i)=0; end end if sum(T)==sum(ones(n,1)) Fg1=1;%输入方阵为严格对角占优矩阵 else Fg1=0;%输入方阵产不是严格对角占优矩阵 end %以下判断它是否为不可约且具有对角优势 %先判断是否为不可约矩阵 D=A; D=rref(D); d=D(n-1,:); if sum(d)~=0 flag=1;%A为不可约矩阵 else flag=0;%A为可约矩阵 end for i=1:n b=B(i,i); temp=B(i,:); temp(i)=[]; if b>=sum(temp) T1(i)=1; else T1(i)=0; end end Fg2=0;%先假定A为可约矩阵且不具有对角优势 50