作业_流程图形式

第一题 列主元Gauss 消去法（1）算法叙述：Gauss 消去法，实际上就是我们俗称的加减消元法，给出矩阵n n ij a ⨯=）（A 以及Tnb b b b )...,(21=，进而得到增广矩阵[A ，b]，通过加减消元法能够得到方程的解)...,(21n x x x x =和相应的消元后主对角线以下元素全部是0的[A ，b]。

（2）算法的具体描述（流程图形式）：求出矩阵行数N找出第i 列中最大值行（i 行到n 行）第i 行是否为最大值行第i 行与最大值行交换否消元处理是求出Xn 的值求出Xn-1...X1的值输出增广矩阵Ab 和X（3）具体算例：输入输出如下>> A=[12 -3 3 4;-18 3 -1 5;3 4 8 13;1 1 1 1]; >> b=[15;-15;6;8]; >> [Ab,x]=gauss[A,b]Ab = x = -18.0000 3.0000 -1.0000 5.0000 -15.0000 11.3117 0 4.5000 7.8333 13.8333 3.5000 21.2532 0 0 4.0741 10.4074 5.7778 -41.2792 0 -0.0000 0 0.4667 7.8000 16.7143（4）总结：使用gauss消去法是，要求主对角元素不能为零，且由于舍入误差较大，使方程的解不“好看”。

附程序：function [Ab,x]=Gauss(A,b)N_num = length(A);x=zeros(N_num,1);for i=1:N_num-1r_max=A(i,i);r_num = 0;for j=i+1:N_num %求出第i列中最大值的行数if(abs(A(j,i))>r_max)r_max = abs(A(j,i));r_num = j;endend;if(r_num >i)%如果A(i,i)非最大，则交换for k=i:N_numA_t = A(i,k);A(i,k) = A(r_num,k);A(r_num,k) = A_t;endB_t = b(i);b(i) = b(r_num);b(r_num) = B_t;end;for m=i+1:1:N_num %进行消元处理temp = A(m,i)/A(i,i);for n = i:1:N_numA(m,n) = A(m,n) - A(i,n)*temp;endb(m) = b(m) - b(i)*temp;endend;Ab = [A,b];x(N_num) = b(N_num)/A(N_num,N_num);%求得Xn的值for k = N_num-1:-1:1%求得其他x 的值for n=k+1:N_numb(k) = b(k) - A(k,n)*x(n);endx(k)=b(k)/A(k,k);end;第二题 一般三对角线性方程组的追赶法（1）对算法的叙述 线性方程组为Ax=f1112222111,0:,02,31n n n n i i ii i n n A f b A i n f b c f a b c c a f b c b a c a c b a -⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦⎧>>⎪⎪≥+≠=-⎨⎪>>⎪⎩ 并且满足此时对矩阵A 可LU 分解⎢⎢⎢⎣⎡=====2,3,...n)(i ,l c -b u 2,3,...n)(i 1,-/u a l b u i 1-i i i i i i 11解Ly=d 得到公式 ⎢⎣⎡=--==),...,3,2(,1l 11n i y d y d y i i i i解Ux=y 得到公式⎢⎣⎡--=+-==)1,2,...,2,1(,/)1(u /n n i u x c y x y x i i i i i nn n（2）对算法的具体描述输入A ，求出行数N输入A ，求出行数N LU 分解Ly=d ；Ux=yLU 分解Ly=d ；Ux=y 根据公式解Ly=d 根据公式解Ly=d 根据公式解Ux=y根据公式解Ux=y 输出解x输出解x图2-1 追赶法的算法流程图（3）具体算例，数值结果>> a=[2 2 1];>> b=[3 3 3 3];>> c=[1 1 1];>> d=[1 0 1 0];>> x=zhuigan(a,b,c,d)x = 0.5526 -0.6579 0.8684 -0.2895（4）总结三对角矩阵是一种特殊的矩阵，其求解类可以使用求解普通的矩阵的LU分解的方法附程序function x=zhuigan(A_down,A_mid,A_up,f)N = length(A_mid);%LU分解u(1) = A_mid(1);for i=2:Nl(i)=A_down(i-1)/u(i-1);u(i)=A_mid(i)-l(i)*A_up(i-1);if(u(i) == 0)s='u(i)=0无法分解';break;end;end%解ly=fy(1) = f(1);for i = 2:Ny(i) = f(i) - l(i)*y(i-1);end%解ux=yx(N)=y(N)/u(N);for i=N-1:-1:1x(i)=(y(i)-A_up(i)*x(i+1))/u(i);End第三题 Jacobi 、Gauss-Seidel 和SOR 迭代算法（1）对算法的叙述 （a ）Jacobi 迭代算法：如果A = D - L- U ，其中 ),...,,(2211nn a a a diag D =可得Jacobi 迭代公式为： 可得：b D x U L D x k k 11)1()()(--+++= k=0,1,2...其中)(1U L D B J +=-称为Jacobi 迭代矩阵。

（b ）GS 迭代算法：在计算时，通过对Jacobi 迭代格式进行修改，可以得到GS 迭代公式为：可得：b L D Ux L D x k k 11)1()()()(--+-+-= k=0,1,2... 其中U L D B G 1)(--= 称为高斯-塞德尔（Gauss-Seidel ）迭代矩阵。

这样能够得到更好的收敛效果。

（c ）SOR 迭代算法在GS 迭代算法基础上，为了得到更好的收敛效果，在修正项前乘以一个参数w ，得到SOR 迭代格式，其中w 为松弛因子。

bwL D w xwU D w wL D Dx Ux Lx b wD x x k k k k k k 11)1(1)1()(])1[()()()()()(--+-+-++--=-+++=])1[()(1U D w wL D B S +--=-称为SOR 迭代矩阵。

（2）对算法的具体描述输入A 、b参数x0、eps （Sor 多一个w ）输入A 、b 参数x0、eps （Sor 多一个w ）求得A 的对角矩阵D下三角L 、上三角U迭代矩阵B 求得A 的对角矩阵D 下三角L 、上三角U 迭代矩阵B 第一次迭代值x第一次迭代值x Norm 函数判定X-xx0<epsNorm 函数判定X-xx0<eps 输出x 和n输出x 和n 是Jacobi 迭代公式GS 迭代公式SOR 迭代公式Jacobi 迭代公式GS 迭代公式SOR 迭代公式否迭代次数n>M迭代次数n>M 否Stop 迭代Stop 迭代是图3-1 Jacobi 、GS 、SOR 算法流程图（3）具体算例，数值结果 （a ）Jacobi 迭代算法： >> A=[4 -1 1;-1 3 -1;0 -1 3]; b=[1;5;-3]; eps = 1e-5; x0 = [0;0;0];[x,n]=jacobi(A,b,x0,eps)输出x = n = 0.8000 19 1.8000 -0.4000 （b ）GS 迭代算法： >> A=[4 -1 1;-1 3 -1;0 -1 3]; b=[1;5;-3]; eps = 1e-5; x0 = [0;0;0]; [x,n]=gs(A,b,x0,eps)输出：x = n = 0.8000 7 1.8000-0.4000（c）SOR迭代算法：>> A=[4 -1 1;-1 3 -1;0 -1 3];b=[1;5;-3];eps = 1e-5;x0 = [0;0;0];w=1.1;[x,n]=SOR(A,b,x0,eps,w)输出：x = n =0.8000 81.8000-0.4000（4）总结由上例知，Jacobi算法、算法和SOR算法最后迭代得到的结果是相同的，但所需要的迭代次数不同，这说明这三种算法的说明收敛速度不同。

在要求相同的误差（5=eeps）条件下，所需迭代次数依次减小，1-由此得出SOR算法的收敛速度最快、Gauss-Seidel算法居中、Jacobi最慢。

其次，在使用SOR方法时，首先要确定权重w附程序程序一：Jacobi算法程序：function [x,n]=jacobi(A,b,x0,eps)M=200;%自定义迭代次数最大迭代次数D=diag(diag(A)); %求A的对角矩阵inv_D = inv(D); %求得D的逆L = -tril(A,-1); %求A的下三角阵U = -triu(A,1); %求A的上三角阵Bj=inv_D*(L+U); %GS迭代矩阵f=inv_D*b;x=Bj*x0+f; %GS迭代n=1; %迭代次数xx0 = x0;while norm(x-xx0)>=epsxx0=x;x=Bj*xx0+f;n = n+1;if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endEnd程序二：GS迭代算法function [x,n]=gs(A,b,x0,eps)M=200;%自定义迭代次数最大迭代次数D=diag(diag(A)); %求A的对角矩阵L = -tril(A,-1); %求A的下三角阵U = -triu(A,1); %求A的上三角阵inv_D_L = inv(D-L); %求得D-L的逆Bgs=inv_D_L*U; %GS迭代矩阵f=inv_D_L*b;x=Bgs*x0+f; %Jacobi迭代n=1; %迭代次数xx0 = x0;while norm(x-xx0)>=epsxx0=x;x=Bgs*xx0+f;n = n+1;if(n>=M)disp('迭代次数太多，可能不收敛！');return;endend程序三：SOR迭代算法程序function [x,n]=SOR(A,b,x0,eps,w)M=200;%自定义迭代次数最大迭代次数D=diag(diag(A)); %求A的对角矩阵L = -tril(A,-1); %求A的下三角阵U = -triu(A,1); %求A的上三角阵inv_D_wL = inv(D-w*L); %求得D-wL的逆Bsor=inv_D_wL*((1-w)*D+w*U); %SOR迭代矩阵f=w*inv_D_wL*b;x=Bsor*x0+f; %SOR迭代n=1; %迭代次数xx0 = x0;while norm(x-xx0)>=epsxx0=x;x=Bsor*xx0+f;n = n+1;if(n>=M)disp('迭代次数太多，可能不收敛！');return;endEnd第四题 幂法和反幂法算法一：幂法（1）幂法：求幂法是一种计算矩阵主特征值及对应特征向量的迭代方法, 特别是用于大型稀疏矩阵。