迭代法实验
实验五 线性方程组的迭代法实验
一. 实验目的
（1）深入理解线性方程组的迭代法的设计思想，学会利用系数矩阵的性质以保证迭 代过程的收敛性，以及解决某些实际的线性方程组求解问题。
（2）熟悉Matlab编程环境，利用Matlab解决具体的方程求根问题。 二. 实验要求
建立Jacobi迭代公式、Gauss-Seidel迭代公式和超松弛迭代公式，用Matlab软件实现线性方程组求解的Jacobi迭代法、Gauss-Seidel迭代法和超松弛迭代法，并用实例在计算机上计算。 三. 实验内容
1. 实验题目
（1）分别利用Jacobi迭代和Gauss-Seidel迭代求解下列线性方程组，取x0={0 ,0,0,0,0-,o}t （2）分别取w=1、1.05、1.1、1.25和1.8，用超松弛法求解上面的方程组，要求精度 为510 。
2. 设计思想
1.Jacobi迭代: Jacobi迭代的设计思想是将所给线性方程组逐步对角化，将一般形式的线性方程组的求解归结为对角方程组求解过程的重复。
2.Gauss-Seidel迭代: Gauss-Seidel迭代的设计思想是将一般形式的线性方程组的求解过程归结为下三角方程组求解过程的重复。
3.超松弛迭代： 基于Gauss-Seidel迭代，对i=1，2，…反复执行计算迭代公式，即为超松弛迭代。
3. 对应程序 1.Jacobi迭代：
function [x,k]=Jacobimethod(A,b,x0,N,emg)
%A是线性方程组的左端矩阵，b是右端向量，x0是迭代初始值
% N表示迭代次数上限，emg表示控制精度，k表示迭代次数，x是解 n=length(A); x1=zeros(n,1); x2=zeros(n,1); x1=x0;k=0;
r=max(abs(b-A*x1)); while r>emg for i=1:n sum=0; for j=1:n
if i~=j
sum=sum+A(i,j)*x1(j);
end end
x2(i)=(b(i)-sum)/A(i,i); end
r=max(abs(x2-x1)); x1=x2; k=k+1; if k>N
disp('迭代失败，返回'); return; end end
x=x1;
2.Gauss-Seidel迭代：
function [x,k]=Gaussmethod(A,b,x0,N,emg)
%A是线性方程组的左端矩阵，b是右端向量，x0是迭代初始值
% N表示迭代次数上限，emg表示控制精度，k表示迭代次数，x是解 n=length(A); x1=zeros(n,1); x2=zeros(n,1); x1=x0;
r=max(abs(b-A*x1)); k=0;
while r>emg
for i=1:n sum=0;
for j=1:n if j>i
sum=sum+A(i,j)*x1(j); elseif j
sum=sum+A(i,j)*x2(j); end end
x2(i)=(b(i)-sum)/A(i,i); end
r=max(abs(x2-x1)); x1=x2; k=k+1; if k>N
disp('迭代失败，返回'); return; end end x=x1;
3.超松弛(SOR)迭代：
function [x,k]=SORmethod(A,b,x0,N,emg,w)
%A是线性方程组的左端矩阵，b是右端向量，x0是迭代初始值
% N表示迭代次数上限，emg表示控制精度，k表示迭代次数，x是解 %w表示松弛因子 n=length(A);
x1=zeros(n,1);x2=zeros(n,1); x1=x0;
r=max(abs(b-A*x1)); k=0;
while r>emg for i=1:n
sum=0; for j=1:n if j>=i
sum=sum+A(i,j)*x1(j); else
if j
end
x2(i)=x1(i)+w*(b(i)-sum)/A(i,i); end
r=max(abs(x2-x1)); x1=x2; k=k+1; if k>N
disp('迭代失败，返回'); return; end end x=x1;
四. 实验体会 在同等精度下，Gauss-Seidel迭代法比Jacobi迭代法收敛速度快。一
般来说，Gauss-Seidel迭代法比Jacobi迭代法收敛要快，但有时反而比Jacobi迭代法要慢，而且Jacobi迭代法更易于优化。因此，两种方法各有优缺点，使用时要根据所需适当选取。 当松弛因子为1时，超松弛迭代方法等同于Gauss-Seidel迭代法，这和理论推导完全相同。另外，超松弛迭代法的收敛速度完全取决于松弛因子的选取，一个适当的因子能大大提高收敛速度。
实验四 线方程组的直接解法
一、问题提出
给出下列几个不同类型的线性方程组，请用适当算法计算其解。 1、 设线性方程组
⎡4⎢8⎢⎢4⎢⎢0⎢-4⎢⎢8⎢0⎢⎢16⎢4⎢⎢⎣0
26
2-22621060
-3-5-216-8-1-112-1
-1-3-15-153-9-78
263-167-41713-3
[1**********]-24
00-1-1-32522-8
01013
00392
0⎤⎡x1⎤⎡5⎤
⎢x⎥⎢⎥0⎥⎥⎢2⎥⎢12⎥1⎥⎢x3⎥⎢3⎥⎥⎢⎥⎢⎥4⎥⎢x4⎥⎢2⎥ 3⎥⎢x5⎥⎢3⎥⎥⎢⎥=⎢⎥5⎥⎢x6⎥⎢46⎥1⎥⎢x7⎥⎢13⎥⎥⎢⎥⎢⎥2⎥⎢x8⎥⎢38⎥4⎥⎢x9⎥⎢19⎥⎥⎢⎥⎢⎥⎢⎥x⎢⎥-1⎥-21⎦⎣10⎦⎣⎦
6-3
30-120126
3
x*=(1,-1,0,1,2,0,3,1,-1,2)T
2、 设对称正定阵系数阵线方程组
⎡4⎢2⎢⎢-4⎢⎢0⎢2⎢⎢4⎢0⎢⎢⎣0
22-1-21320
-4-1141-8-356
0-216-1-4-33
21-8-1224-10-3
43-3-44111-4
025-3-101142
0⎤⎡x1⎤⎡0⎤
⎢x⎥⎢⎥0⎥⎥⎢2⎥⎢-6⎥6⎥⎢x3⎥⎢20⎥ ⎥⎢⎥⎢⎥3⎥⎢x4⎥⎢23⎥
=
-3⎥⎢x5⎥⎢9⎥
⎥⎢⎥⎢⎥-4⎥⎢x6⎥⎢-22⎥2⎥⎢x7⎥⎢-15⎥⎥⎢⎥⎢⎥
x19⎥45⎢⎥⎢⎥⎦⎣8⎦⎣⎦
x*=(1,-1,0,2,1,-1,0,2)T
三对角形线性方程组
0000000⎤⎡x1⎤⎡7⎤⎡4-10
⎢x⎥⎢⎢-14-10⎥000000⎥⎢⎥⎢2⎥⎢5⎥
⎢0-14-1000000⎥⎢x3⎥⎢-13⎥⎢⎥⎢⎥⎢⎥
x00-14-10000024⎥⎢⎥⎢⎥ ⎢
⎢⎥⎢0⎥⎢x500-14-100006⎥
⎢⎥⎢⎥=⎢⎥
000-14-1000⎥⎢x6⎥⎢-12⎥⎢0⎢00000-14-100⎥⎢x7⎥⎢14⎥⎢⎥⎢⎥⎢⎥
x000000-14-10-4⎢⎥⎢8⎥⎢⎥⎢⎥⎢0⎥⎢000000-14-1x95⎥
⎢⎥⎢⎥⎢⎥
⎢⎥x⎢⎥⎢⎥00000000-14-5⎣⎦⎣10⎦⎣⎦
x*=(2,1, ,1,1)-3,0-,1,2,3-,T0
二、要求
1、 对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法；平方根法与改进平方根法；追赶法求解（选择其一）； 2、 应用结构程序设计编出通用程序；
3、 比较计算结果，分析数值解误差的原因；
4、 尽可能利用相应模块输出系数矩阵的三角分解式。
三、目的和意义
1、通过该课题的实验，体会模块化结构程序设计方法的优点； 2、运用所学的计算方法，解决各类线性方程组的直接算法； 3、提高分析和解决问题的能力，做到学以致用；
3、 通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点。
四、实验学时：2学时 五、实验步骤：
1．进入C或matlab开发环境； 2．根据实验内容和要求编写程序； 3．调试程序； 4．运行程序；
5．撰写报告,讨论分析实验结果．
实验五 解线性方程组的迭代法
一、问题提出
对实验四所列目的和意义的线性方程组，试分别选用Jacobi 迭代法，Gauss-Seidel迭代法和SOR方法计算其解。
二、要求
1、体会迭代法求解线性方程组，并能与消去法做以比较；
2、分别对不同精度要求，如ε=10-3,10-4,10-5由迭代次数体会该迭代法的收敛快慢； 3、对方程组2，3使用SOR方法时，选取松弛因子ω=0.8，0.9，1，1.1，1.2等，试看对算法收敛性的影响，并能找出你所选用的松弛因子的最佳者； 4、给出各种算法的设计程序和计算结果。
三、目的和意义
1、通过上机计算体会迭代法求解线性方程组的特点，并能和消去法比较； 2、运用所学的迭代法算法，解决各类线性方程组，编出算法程序； 3、体会上机计算时，终止步骤x
敛散性的意义；
4、 体会初始解x，松弛因子的选取，对计算结果的影响。
(k+1)
-xk
∞
（给予的迭代次数），对迭代法
四、实验学时：2学时 五、实验步骤：
1．进入C或matlab开发环境；
2．根据实验内容和要求编写程序； 3．调试程序； 4．运行程序；
5．撰写报告,讨论分析实验结果．
例3 例3 用平方根法分解对称正定矩阵
-11⎤⎡4
A=⎢-14.252.75⎥
⎢⎥⎢⎣12.753.5⎥⎦
l11=a11=4=2l21=
解
a21-1==-0.5l112 a1
l31=31==0.5
l112
2
l22=a22-l21=4.25-0.25=2
a-ll2.75-0.5(-0.5)l32=323121==1.5l222 22
l33=a33-l31-l32=.5-0.25-2.25=1
于是A=LL，其中
T
1
n(n+1)
由于A为对称矩阵,因此,在电算时只要存储A的下三角部分,其需要存储2个元素,
可用一维数组存放,即
00⎤⎡2
L=⎢-0.500⎥
⎢⎥⎢⎣0.51.51⎥⎦
⎡1⎤
A⎢n(n+1)⎥={a11,a21,...,an1,an2,...ann}⎣2⎦ ⎡1⎤1A⎢n(n+1)⎥i(i-1)+ja2ij⎦的第2矩阵元素存放在⎣个位置,L的元素存放在A的相应位
置上.另外,平方根法的运算量是
开平方 n次；
当n比较大时,平方根法的运算量和存贮量约为高斯消元法的二分之一，因此它是求解对称正定
矩阵比较好的方法.
为了避免开方运算,我们可以采用下面的分解式
13321n+n+n
23次； 乘除法 6137n+n2-n
6次. 加减法 6
A=LDLT
其中L是单位下三角阵,D是对角阵,由矩阵乘法,可得L与D的计算公式.
对于i
(24)
(25)
=1,2,...,n,有
k-1j=1
lik=(aik-∑lijdjlkj)
k
i-1j=1
k=1,2,...,i-1
di=aii-∑lij2dj,
为了避免重复计算,我们引入
(26)
tij=lijdj
于是上述公式可改写成
对于i
(27)
k=1,2,...,i-1
=1,2,...,n,有
tik=aik-∑tijlkj,
j=1
k-1
(28)
(29)
tik
lik=,
dk
i-1j=1
k=1,2,...,n
di=aii-∑tijlij,
计算出T
(30)
=LD的第i行元素tik,k=1,2,...i-1后,存放在A的第i行相应位置，然后再
ltalt计算L的第i行元素ik仍然存放在A的第i行，即用ik冲掉ik,再用ik冲掉ik,D的对角线
元素存放在A的相应位置上.
对称正定矩阵A按LDL的分解和按LL分解其计算量差不多,但LDL分解不需要开方计算，它称为改进的平方根法. 四 追赶法
在计算样条函数,解常微分方程边值问题,解热传导方程等都会要求解系数矩阵呈三对角线T
T
T
形的线性方程组,这时
⎡a11a⎢12
⎤a⎥22a23A=⎢a21
⎢
⎥⎢
⎥⎢an-1n-2
an-1n-1a⎥n-1n⎢⎥⎣
ann-1
ann⎥⎦
的LU分解中,矩阵L和U分别取下二对角线和上二对角线形式,设
⎡⎢l11⎤⎡1u12⎤
L=⎢l⎥⎢21l ⎥22
⎢⎢ ⎥⎥U=⎢
⎢1u⎥n-1n⎥⎣l⎥⎢⎥nn-1lnn⎦, ⎣1⎦由A=LU得计算公式
a11=l11 aii-1=lii-1,i=2,3,...,n aii=lii-1ui-1i+lii,,i=2,3,...,n aii+1=li1uii+1,i=1,2,...,n-1即
l11=a11
ua12
12=l11
lii-1=aii-1
lii=aii-lii-1ui-1i
uaii+1
ii+1=
liii=2,3,...,
nb
此时,求解Ax=等价于解两个二对角线方程组
⎧⎨
Ly=b⎩Ux=y
自上而下解方程组
Ly=b形象地称为“追”.
y11=
bl11
(31)
yi=
自下而上解方程组
bi-lii-1yi-1
,lii
i=2,3,...n
(32)
Ux=y称为“赶”.
i=n-1,...,2,1
xn=yn
xi=yi-uii+1xi+1,
习惯,上述求解方法称为“追赶法”. 例4 用追赶法解三对角线方程组
(33)
解 由三对角分解公式有
=1⎧2x1-x2
⎪-x+2x-x=0⎪123⎨
-x2+2x3-x4=0⎪⎪-x3+2x4=1 ⎩
而由“追”公式有
l11=a11=2
u12=a12l11=-2 l21=a21=-1
l22=a22-l21u12=2-2=2 u23=a23l22=-3 l32=a32=-1
l33=a33-l32u23=3 u34=a34l33=-4 l43=a43=-1
l44=a44-l43u34=4
y1=
b11=l112
y2=
b2-l21y11
=l223
y3=
b3-l32y21
=l334
最后，由“赶”公式得原方程组的解
b4-l43y3
y4==1
l44
x4=y4=1x3=y3-u34x4=1x2=y2-u23x3=1 x1=y1-u12x2=1
追赶法公式实际上就是把高斯消元法用到求解三对角线方程组上去的结果，这时由于A特别
简单，因此使得求解的计算公式非常简单，而且计算量仅有5n-4次乘除法,3n-3次加减法,
仅占5n-2个存贮单元,所以可以在小机器上解高阶三对角线形的线性代数方程组.
求解线性方程组的直接解法
二 实验部分
本章实验内容：
实验题目：Gauss消元法,追赶法,范数。
实验内容：①编制用Gauss消元法求解线性方程组Ax=f的程序。
②编制用追赶法求解线性方程组Ax=f的程序。 ③编
制向量和矩阵的范数程序。
实验目的：①了解Gauss消元法原理及实现条件,熟练掌握Gauss消元法解方程
组的算法,并能计算行列式的值。
②掌握追赶法，能利用追赶法求解线性方程组。
③理解向量和矩阵范数定义,性质并掌握其计算方法.
编程要求：利用Gauss消元法,追赶法解线性方程组。分析误差。 计算算法:①Gauss消元法:
1. 消元过程
(k)设akk≠0,对k=1,2, ,n-1计算 (k)(k)⎧mik=aik/akk⎪(k+1)
(k)(k)
⎨aij=aij-mikakj
⎪(k+1)(k)(k)
=bi-mikbk⎩bi
i,j=k+1,k+2, ,n
⒉回代过程
(n)(n)
⎧xn=bn/ann⎪n⎨(i)(i)(i)x=(b-ax)/a∑iijjii ⎪i
j=i+1⎩
i=n-1, ,2,1
②追赶法：
1.分解Ax=f: β1=c1/b1,βi=ci/(bi-aiβi-1)
( i=2,3, ,n-1)
2.解Lg=f,求g:
g1=f1/b1,gi=(fi-aigi-1)/(bi-aiβi-1)
(i=2,3, ,n)
3．解Ux=g,求x:xn=yn,xi=yi-βixi+1 （i=n-1,n-2, ,2,1） ③范数：
常用向量范数有：（令x=( x1,x2,…,xn)T） 1-范数： ║x║1=│x1│+│x2│+…+│xn│
2-范数： ║x║2=(│x1│2+│x2│2+…+│xn│2)1/2 ∞-范数： ║x║∞=max(│x1│,│x2│,…,│xn│) 常用的三种向量范数导出的矩阵范数是：
1-范数:║A║1= max{║Ax║1/║x║1=1}=max∑ajj 1≤j≤ni=1n
2-范数:║A║2=max{║Ax║2/║x║2=1}=1,λ1是ATA的最
大特征值.
∞-范数:║A║∞=max{║Ax║∞/║x║∞=1}=max∑aij 1≤i≤nn
实验例题⑴：用Gauss消元法解方程组
⎛2-11⎫⎛x1⎫⎛4⎫ ⎪ ⎪ ⎪ -1-23⎪ x2⎪= 5⎪
1 ⎪ ⎪31⎪⎝⎭⎝x3⎭⎝6⎭j=1
实验例题⑵: 用追赶法解三对角方程组Ax=b,其中
00⎤⎡2-10⎡1⎤⎢-12-10⎥⎢0⎥0⎢⎥⎢⎥ A=⎢0-12-10⎥,b=⎢0⎥ ⎢⎥⎢⎥00-12-1⎢⎥⎢0⎥⎢⎢00-12⎥⎣0⎦⎣0⎥⎦
实验例题⑶：设
⎛0.60.5⎫ A= 0.10.3⎪⎪， ⎝⎭
计算A的行列范数.
程序①:Gauss消元法
function x=Gauss(A,b)
%A是线性方程组的系数矩阵,b为自由项.
n=length(A)
for k=1:n-1
m(k+1:n,k)=A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m(k+1:n,k)*A(k,k+1:n);
b(k+1:n)=b(k+1:n)-m(k+1:n,k)*b(k);
end
x=zeros(n,1);
x(n)=b(n)/A(n,n);
for k=n-1:-1:1
x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);
end
x=x';
disp(sprintf('k x(k)'));
for i=0:n
disp(sprintf('%d %f ',i,x(i+1)));
end
数值结果：x=Gauss(A,b)
n =3
k x(k)
0 1.111111
1 0.777778
2 2.555556
程序②:追赶法
function [x,y,beta]=zhuiganfa(a,b,c,f)
%a,b,c是三对角阵的对角线上的元素,f是自由项.
n=length(b);
beta(1)=c(1)/b(1);
for i=2:n
beta(i)=c(i)/(b(i)-a(i)*beta(i-1));
end
y(1)=f(1)/b(1);
for i=2:n
y(i)=(f(i)-a(i)*y(i-1))/(b(i)-a(i)*beta(i-1));
end
x(n)=y(n);
for i=n-1:-1:1
x(i)=y(i)-beta(i)*x(i+1);
end
disp(sprintf('k x(k) y(k) beta(k)')); for i=0:n
disp(sprintf('%d %f ',i,x(i+1),y(i+1),beta(i+1))); end
数值结果:
a=[0 -1 -1 -1 -1]';
b=[2 2 2 2 2]';
c=[-1 -1 -1 -1 0]';
f=[1 0 0 0 0]';
[x,y,beta]=zhuiganfa(a,b,c,f)
k x(k) y(k) beta(k)
0 0.833333 5.000000e-001 -0.500000
1 0.666667 3.333333e-001 -0.666667
2 0.500000 2.500000e-001 -0.750000
3 0.333333 2.000000e-001 -0.800000
4 0.166667 1.666667e-001 0.000000