设计题目:超松弛迭代法解线性方程组
摘要
本文是在matlab环境下熟悉的运用计算机编程语言并结合超松弛变量超松弛迭代法的理论基础对方程组求解。
首先,本文以微分方程边值问题为例,导出了离散化后线性方程组即稀疏线性方程组,转化对稀疏线性方程组求解问题。其次,用超松弛( SOR) 迭代法编写matlab程序,对产生的稀疏线性方程组进行迭代法求解。然后,分别改变松弛因子ω和分段数n的值,分析其收敛性和收敛速度,做出各个方面的分析和比较得到相关结论。最后,将超松弛迭代算法在计算机上运用matlab语言实现, 得出了一组与精确解较接近的数值解,并画图比较,验证逐次超松弛( SOR) 迭代法的精确性。
关键词:稀疏线性方程组逐次超松弛迭代法松弛因子
matlab编程
一、问题提出
考虑两点边值问题
()()?????==<<=+.
11,00,
10,22y y a a dx
dy dx y d ε 容易知道它的精确解为
.1111ax e e a y x +???
? ??
---=
--εε
为了把微分方程离散,把[]1,0区间n 等分,令n
h 1
=,ih x i =,,1,,2,1-=n i 得到差分方程
,212
11a h y y h
y y y i i i i i =-++-++-ε 简化为
()(),2211ah y y h y h i i i =++-+-+εεε
从而离散后得到的线性方程组的系数矩阵为
()()()()?????
??
?
??
??????+-++-++-++-=h h h h h h h A εεεεεεεεεε2222
对1=ε,4.0=a ,200=n ,分别用1=ω、5.0=ω和5.1=ω的超松弛迭代法求解线性方程组,要求有4位有效数字,然后比较与精确解的误差,探讨使超松弛迭代法收敛较快的ω取值,对结果进行分析。改变n ,讨论同样问题。
二、超松弛迭代法产生的背景
对从实际问题中得到维数相当大的线性代数方程组的求解仍然十分困难, 以
至使人们不能在允许的时间用直接方法得到解, 因此, 客观上要求用新的方法来
解决大维数方程组的求解问题。
现有大多数迭代法不是对各类线性方程组都有收敛性, 在解题时, 要对原方程组矩阵作一根本的变换, 从而可能使条件数变坏, 也可能破坏了变换前后方程组的等价性, 以及丧失使原方程组的对称性等。探求新的有效的解题方法依然是迫切的任务。逐次超松弛(Successive Over Relaxation)迭代法是在高斯-赛德尔
(GS)迭代法基础上为提高收敛速度,采用加权平均而得到的新算法。
在求解过程中由于线性方程组的系数矩阵维数较大, 采用计算机编写算法来求解, 从而实现了对解析模型的计算机数值逼近的计算方法#本论文以逐次超松弛迭代法为主要的求解方法。
三、超松弛迭代法的理论基础
(一)逐次超松弛迭代法
逐次超松弛(Successive Over Relaxation)迭代法,简称SOR迭代法,它是在GS法基础上为提高收敛速度,采用加权平均而得到的新算法,设解方程
(7.1.3)的GS法记为
(1)
再由与加权平均得
这里ω>0称为松弛参数,将(1)代入则得
(2)
该法称为SOR迭代法,[WTBX]ω>0称为松弛因子,当ω=1时(2)式即为高斯-赛德尔迭代法,简记GS法,将(2)写成矩阵形式,则得
即
于是得SOR迭代的矩阵表示
(3)
其中
(4)
分解后,有.
(二)逐次超松弛迭代法的收敛性
根据迭代法收敛性定理,SOR法收敛的充分必要条件为,收敛的充分条件为,但要计算比较复杂,通常都不用此结论,而直接根据方程组的系数矩阵A 判断SOR迭代收敛性,下面先给出收敛必要条件.
定理1设,则解方程的SOR迭代法收敛的必要条件是0<ω<2.
该定理为SOR迭代法收敛的必要条件。
定理 2若对称正定,且0<ω<2,则解Ax=b的SOR迭代法对迭代收敛.
对于SOR迭代法,松弛因子的选择对收敛速度影响较大,关于最优松弛因子研究较为复杂,且已有不少理论结果.下面只给出一种简单且便于使用的结论。
定理3设为对称正定的三对角矩阵,是解方程的J法迭代矩阵,若,记,
则SOR法的最优松弛因子为
(5)
且
(6)
根据定理,,如图1所示.由(6)可知,当ω=1,时,收敛速度为
.
说明GS法比J法快一倍.
图1
定理4设,如果:
(1)A为严格对角占优矩阵;(2)0<ω<=1.
则解的SOR迭代法收敛。
四、实验容
1.自定义函数 sor(A, b, nm, e, w),以实现SOR方法求解线性方程组AX=B,其中
A——系数矩阵;
b ——常数列向量; w ——松弛因子; nm ——迭代的最大次数
e
——(1)
()k k X
X +∞
-达到的精度上限
由离散后的差分方程:
()(),2211ah y y h y h i i i =++-+-+εεε,1,,2,1-=n i
得到的线性方程组的系数矩阵为
()()()()?????
??
?
??
??????+-++-++-++-=h h h h h h h A εεεεεεεεεε2222
常数列向量b=???
??
??
?
????????????????+--y ah
ah ah y ah h 2002
2
202
)(..
.εε 其中
1
=ε,
4
.0=a ,
200
=n ,
n
h 1
=
,则有
00001.0,005.1,005.2)2(2==+-=+-ah h h εε。A 为(a ij )200*200型矩阵,b 为(bij )200*1
型矩阵。
在本次试验中,由于所提供数据较小,当最大迭代次数nm 较小时,在nm 迭代次数围,不能判断该超松弛迭代法是否收敛,此次取nm=30000。迭代精度e 也应取较小值才能使误差更小,此次取e=0.00001。由定理1可知,本次试验中,ω的取值围为:0<ω<2才能保证迭代法收敛。
取T x )1,1,,1,1,1()0( =,为1200?的矩阵。用SOR 迭代公式得
()
()()(
)
()
(
)
?????
?
???????+---=-+---=-+---=-+---=-+---=-+--=+++++++++++.005.2/005.200001.0;
005.2/005.1005.200001.0;005.2/005.1005.200001.0;
005.2/005.1005.200001.0;005.2/005.1005.200001.0;005.2/005.1005.200001.0)
(222)1(199)
(200)1(200)(200
)(199)1(198)(199)1(199)
(5)(4)1(3)(4)1(4
)(4
)(3)1(2)(3)1(3)(3)(2)1(1)(2)1(2
)
(2)(1)(1)1(1k k k k k k k k k k k k k k k k k k k k k k k k k k k k x x x x x x x x x x x x x x x x x x x x x x x x x x x x ωωωωωω ω取不同值时,对应的迭代次数、与精确解的误差如下表1。
表1 ω取不同值时对应的迭代次数与误差
图1 计算值与精确值图形比较
从本组的实验中,可以看出w 值的取定十分重要,它对求解的迭代次数影响十分明显。一个不好的w 值甚至会导致迭代超过10000次仍未能求得需要精度的值。
由表1可得,当ω=1.9时SOR 迭代法收敛速度最快,误差最小。
取ω=1.9,nm=30000,T x )6.0,6.0,,6.0,6.0,6.0()0( =等各个因子相同时,当分段点n 取不同值时,对应的迭代次数、与精确解的误差如下表2。
区间等分数n 松弛因子 ω=1.9
迭代次数
误差
120 416 0.2477 150 416 0.2477 200 416 0.2477 250 416 0.2477
满足误差
()42
*10-<-x x k 的迭代次数
从本组的实验中,当其他各个因子取适当值时,改变分段数时,对结果没有影响。
图3 精确图形
由图3可得,当各个参数取值适当时,用SOR迭代法所得线性方程组的解与精确解误差极小,从而验证了SOR迭代法的准确性。
五、结论
1.通过本次的课程设计,可知逐次超松弛迭代法与Jacobi 迭代法, Seidel 迭代
法相比, 收敛速度较快。由逐次超松弛迭代法求出的方程组的数值解与该方程组的精确解十分接近, 离散化后线性方程组的逐次超松弛迭代法的精确性较高。逐次超松弛迭代法可以广泛地应用于实际。该算法不仅可以用来求解高阶稀疏线性方程组, 还可以用来求解热传导问题这样可以大大减少计算量和计算机的存储量, 从而提高计算效率。本次的课程设计,我们运用了matlab语言来实现相关的计算,这样不仅对逐次超松弛迭代法有了更深层的了解掌握,还提高了对matlab的操作技术,深刻体会到了MATLAB功能的强大之处。通过本次试验,我掌握了用Jacobi、Gauss-Seidel、SOR迭代法求解线性方程组的方法;
六、参考文献
[1]庆扬,王能超,易大义.数值分析[M], 清华大学,2008.
[2]卫国. MATLAB程序设计与应用[M],高等教育,2008.
[3]王诗然. 稀疏线性方程组求解的逐次超松弛迭代法[J],师大学学
报,4,407-409,2006.
[4]建宇,黎燕. 牛顿一SOR迭代方法中最佳松弛因子的算法[J],大学学报,
4,381-382,1995.
[5]蔡大用.数值分析与实验学习指导[M],清华大学,2001.
附录
1.超松弛迭代法
function [n,x]=cscdd(A,b,X,nm,e,w)
n=1;
D=diag(diag(A)); %令A=D-L-U,计算矩阵D
L=tril(-A)+D; %令A=D-L-U,计算矩阵L
U=triu(-A)+D; %令A=D-L-U,计算矩阵U
M=inv(D-w*L)*((1-w)*D+w*U); %计算迭代矩阵
g=w*inv(D-w*L)*b; %计算迭代格式中的常数项
%下面是迭代过程
while n<=nm
x=M*X+g; %用迭代格式进行迭代
r=norm(x-X,'inf');
if r return; end X=x; n=n+1; end disp('在最大迭代次数不收敛!') 2.输入初始值并调用SOR迭代法 n0=200;m=1;a=0.4;h=1/n0;A=zeros(n0,n0); for i=1:n0 A(i,i)=-(2*m+h); end for i=2:n0-1 A(i,i-1)=m; A(i,i+1)=m+h; end A(1,2)=m+h;A(n0,n0-1)=m; for i=1:n0-1 b(i,1)=a*h^2; end b(n0,1)=a*h^2-(m+h); for i=1:n0 xi=i/n0; y0(i,1)=((1-0.4)/(1-exp(-1)))*(1-exp(-xi))+0.4*xi; x0(i,1)=1; end [n1,x1]=cscdd(A,b,x0,30000,0.00001,0.5); n1 u1=norm(x1-y0,2) [n2,x2]=cscdd(A,b,x0,30000,0.00001,1.0); n2 u2=norm(x2-y0,2) [n3,x3]=cscdd(A,b,x0,30000,0.00001,1.5); n3 u3=norm(x3-y0,2) t=1/200:1/200:1; plot(t,y0) hold on plot(t,x1,'g') hold on plot(t,x2,'r') hold on plot(t,x3,'k') legend('精确解','w=0.5','w=1','w=1.5'); title('计算值与精确值图形比较') hold off