数值分析上机试验——解线性方程组讲解

1、实验报告课程名称数值分析实验项目名称解线性方程组实验类型上机实验学时4班级20111131学号2011113130姓名张振指导教师沈艳实验室名称理学楼 407实验时间2013.12.9实验成绩预习部分实验过程 表现实验报告 部分总成绩教师签字日期哈尔滨工程大学教务处 制实验四 解线性方程组解线性方程组的基本思想1直接三角分解法：将系数矩阵 A 转变成等价两个矩阵 L 和 U 的乘积 ，其中 L 和 U 分别是下 三角和上三角矩阵。 当 A 的所有顺序主子式都不为 0时，矩阵 A 可以分解为 A=LU ，且分解唯一。其中 L 是单位下三角矩阵， U 是上三角矩阵。2平方根法：如果矩阵 A 为 n

2、 阶对称正定矩阵，则存在一个对角元素为正数的下三角实 矩阵 L ，使得： A=LLT 。当限定 L 的对角元素为正时，这种分解是唯一的，称为 平 方根法( Cholesky )分解。3追赶法：设系数矩阵为三对角矩阵000000000an 1bn 1cn 10anbnb1 c1 0 a2 b2 c2 0 a3 b3A可先依次求出000 000b100001100022000012000330000100L,U000n100000n1000nn0001则方程组 Ax=f 称为三对角方程组。设矩阵 A 非奇异， A 有 Crout 分解 A=LU ，其中 L 为下三角矩阵， U 为单位上三角矩阵，记

3、L ，U 中的元素后，令Ux=y ，先求解下三角方程组Ly=f 得出 y，再求解上三角方程组 Ux=y 。4雅克比迭代法：首先将方程组中的系数矩阵 A 分解成三部分， 即：A = L+D+U ，如图 1所示， 其中 D 为对角阵， L 为下三角矩阵， U 为上三角矩阵。之后确定迭代格式， X（k1） = BX （k） +f ，如图2所示，其中 B称为迭代矩阵， 雅克比迭代法中一般记为 J。（k = 0,1 ， ）再选取初始迭代向量 X （0），开始逐次迭代。5超松弛迭代法（ SOR） 它是在 GS 法基础上为提高收敛速度，采用加权平均而得到的新算法。 选取分裂矩阵 M 为带参数的下三角矩阵1M

4、 （ D L ），其中 0 为可选择的松弛因子，一般当 1 2时称为超松弛。. 实验题目及实验目的1（第五章习题 8）用直接三角分解（杜利特尔（ Doolittle ）分解）求线性方程 组111x1+x2+ x3= 9，452631118，x1+x2+ x3=342531x1+x2+ 2x3 =82的解。2（第五章习题 9）用追赶法解三对角方程组 Ax=b，其中210001121000A=01210，b=00012100001203（第五章习题 10）用改进的平方根法解线性方程组11x123x231x32114（第六章习题 7）用 SOR方法解线性方程组（分别取松弛因子 =1.03 ，=1，

5、=1.1 ）4x1 -x2- x1 +4 x2- x3= 4，- x2 +4 x3= -3.精确解 x*=( 1，1，- 1)T.要求当 x* x(k) 510 6时迭代终止，并且对每 22一个 值确定迭代次数 .5. (第六章习题 8)用 SOR方法解线性方程组(取 =0.9 )5x1 -2 x2 + x3= -12 ，- x1 +4 x2 - 2 x3= 20 ，2x1 -3 x2 +10x3= 3.要求当 x(k 1) x(k) 10 4时迭代终止 .6(第六章习题 9)设有线性方程组 Ax=b，其中 A 为对称正定阵，迭代公式x(k 1) x(k)+(b- A x(k)，k=0，1，2

6、，试证明当 0 2 时上述迭代法收敛(其中 0 A=1/4,1/5,1/6;1/3,1/4,1/5;1/2,1,2; b=9;8;8; x=ZJsanjiao(A,b)2.追赶法(文件 ZG_SDJ.m)function x=ZG_SDJ(a,b,c,f)%a ?a?%b?a?ya ? %c?a?ya ?%f?3 y?bN=length(a);b=b,0;c=0,c;a1=zeros(N,1);b1=zeros(N,1);y=zeros(N,1);x=zeros(N,1);a1(1)=a(1);b1(1)=b(1)/a1(1);y(1)=f(1)/a1(1);for j1=2:Na1(j1)=

7、a(j1)-c(j1)*b1(j1-1);b1(j1)=b(j1)/a1(j1);temp1=f(j1)-c(j1)*y(j1-1);y(j1)=temp1/a1(j1);endj1=N;x(j1)=y(j1);for j1=N-1:-1:1x(j1)=y(j1)-b1(j1)*x(j1+1);end控制台输入代码： a=2 2 2 2 2; b=-1 -1 -1 -1; c=-1 -1 -1 -1; f=1;0;0;0;0; x=ZG_SDJ(a,b,c,f)3.改进的平方根法(文件 GJPFG.m) function GJPFG(A,b) n=length(b);% n?a D?% LDL

8、 ?a? d(1)=A(1,1);for i=2:nfor j=1:i-1 sum1=0;for k=1:j-1 sum1=sum1+t(i,k)*l(j,k);endt(i,j)=A(i,j)-sum1; l(i,j)=t(i,j)/d(j); endsum2=0;for k=1:i-1 sum2=sum2+t(i,k)*l(i,k); endd(i)=A(i,i)-sum2;endfor i=1:nl(i,i)=1;enddisp( ￥ ?y?L?ao); %?a3?￥ ?y?L? ldisp( ? D?ao ); %?a3? D?d% LDLx=b? ?ax?% Ly=b ?y?% Lx

9、=inv D?y ?ax? y(1)=b(1);for i=2:nsum3=0;for k=1:i-1sum3=sum3+l(i,k)*y(k);endy(i)=b(i)-sum3;endx(n)=y(n)/d(n);for i=n-1:-1:1sum4=0;for k=i+1:nsum4=sum4+l(k,i)*x(k);endx(i)=(y(i)/d(i)-sum4;enddisp( Ly=b? ?ay ? o); ydisp( Ax=b ?ax?a o );x 控制台输入代码： A=2 -1 1;-1 -2 3;1 3 1; b=4;5;6; GJPFG(A,b) 4.SOR方法(文件

10、SOR_1.m) function SOR_1(A,b,x0,x_a,w) %x_a?a? ?aif (w=2)error(2? y ? ?); return ;endeps=5.0e-6;D=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);%? A ? %?A?y? %? A ? ? y?B=inv(D-L*w)*(1-w)*D+w*U);f=w*inv(D-L*w)*b;x=B*x0+f;n=1;% ?ywhile norm(x_a-x)=eps x0=x;x =B*x0+f;n=n+1;if (n=200) ?y? ?2?);2?disp( Warning

11、:returnendenddisp( Ax=b ?a?a o );xdisp( ? y?a); o n控制台输入代码： A=4 -1 0;-1 4 -1;0 -1 4; b=1;4;-3; x0=0;0;0; x_a=0.5;1;-0.5; w=1.03; SOR_1(A,b,x0,x_a,w) w=1; SOR_1(A,b,x0,x_a,w) w=1.1; SOR_1(A,b,x0,x_a,w) 5 SOR方法(文件 SOR_2.m) function SOR_2(A,b,x0,w,eps) if (w=2)error( 2? y ? ?); return ;endD=diag(diag(A

12、);L=-tril(A,-1);U=-triu(A,1);%? A ? %?A?y? %? A ? ? y?B=inv(D-L*w)*(1-w)*D+w*U); f=w*inv(D-L*w)*b; x=B*x0+f;n=1;% ?ywhile norm(x-x0)=eps x0=x;x =B*x0+f;n=n+1;if (n=200)disp( Warning: ?y? ?2?);2?returnendend disp( Ax=b ?a?a o );xdisp( ? y?a); on控制台输入代码： A=5 2 1;-1 4 2;2 -3 10; b=-12;20;3; x0=0;0;0; w

13、=0.9; eps=10e-4; SOR_2(A,b,x0,w,eps)6此题为证明题，无程序代码。7. 雅克比迭代法(文件 Jocabi.m )function x=Jocabi(n)A=hilb(n); x_a=ones(n,1); b=A*x_a;eps=1e-4;n=length(b);N=50;x=zeros(n,1);y=zeros(n,1);for k=1:Nsum=0;for i=1:ny(i)=(b(i)-A(i,1:n)*x(1:n)+A(i,i)*x(i)/A(i,i); endfor i=1:n sum=sum+(y(i)-x(i)2; endif sqrt(sum)e

14、psbreak ;elsefor i=1:nx(i)=y(i);endendend SOR方法(文件 SOR_3.m) function SOR_3(n,w) %x_a?a? ?aif (w=2)error(2? y ? ?); return ;end x0=zeros(n,1);A=hilb(n); x_a=ones(n,1); b=A*x_a;eps=1e-4;D=diag(diag(A);%? A ? L=-tril(A,-1);%?A?y?U=-triu(A,1);%? A ? ? y?B=inv(D-L*w)*(1-w)*D+w*U);f=w*inv(D-L*w)*b;x=B*x0+

15、f;n=1;% ?y while norm(x-x0)=epsx0=x;x =B*x0+f;n=n+1;if (n=2000)disp( Warning: ?y? ?2?);2? return ;endenddisp( Hx=b ?a?a o );xdisp( ? y?a); on控制台输入代码： x=Jocabi(6) x=Jocabi(8) x=Jocabi(10) SOR_3(6,1) SOR_3(6,1.25) SOR_3(6,1.5) SOR_3(8,1) SOR_3(8,1.25) SOR_3(8,1.5) SOR_3(10,1) SOR_3(10,1.25) SOR_3(10,1

16、.5)五实验结果比较与分析1.2.3.4.5.9. 证：x(k 1) (IA)x(k) + b，(k=0，1，2)故迭代矩阵 B=I- A，其特征值 =1- (A).由| |1 ，|1-(A)|1 得20 (A)故当 0 2 时，更有 0 2 ，从而有 | |1， ( B) 1，迭代格式收敛 (A)7.雅克比迭代法：这是因为希尔可以看到用雅克比迭代法求希尔伯特阵方程组的解是病态的， 伯特阵的谱半径大于 1，并不收敛SOR迭代法：n=6 , 10 4 时，松弛因子迭代次数近似解1620（ 1.0005 ，1.0045 ，0.9626，1.0441， 1.0285，0.9583）1.25588（

17、0.9997，1.0134 ，0.9362，1.0706， 1.0230，0.9555）1.5539（ 0.9991，1.0211，0.9082，1.1146，0.9899， 0.9656）n=8 , 10 4 时松弛因子迭代次数近似解1426（0.9970，1.0417，0.8967，1.0155，1.0654，1.0505 ，0.9991， 0.9309）1.251157（0.9971，1.0353，0.9174，1.0167，1.0408，1.0378，1.0022，0.9508）1.51701（0.9980，1.0232，0.9484，1.0045，1.0260，1.0324，1.0037，0.9623）n=10 , 10 4 时松弛因子迭代次数近似解11216（0.9977，1.0203，0.9797，0.9654，1.0010，1.0300，1.0367，1.0223，0.9924，0.9525）1.2