matlab第三版第七章,数值分析matlab版第七章

数值分析matlab版第七章

实验报告七 题目：Gauss列主元消去法和LU三角分解法解线性方程组 摘要：掌握用Gauss列主元消去法和LU三角分解法设计程序，从而实现线性组方程求解 原理： 1. Gauss列主元消去法是通过Gauss消去，将线性方程组的系数矩阵化为简单的上三角阵，然后进行回代求出x； 2.LU三角分解法的实质是将方程组的系数矩阵化为一个下三角矩阵L和一个上三角阵U，利用A=LU将Ax=b化简为LUx=b,进一步化简为Ly=b和Ux=y。 例题7.9 Hilbert法： h=hilb(i); condA(i)=cond(h,2); end disp( n cond ); for i=3:10 s=sprintf( %d %f ,i,condA(i)); disp(s); End >> Untitled2 n cond 3 524.056778 4 15513.738739 5 476607.250242 6 14951058.641005 7 475367356.370392 8 15257575270.772364 9 493154209891.651430 10 16024675274036.547000 LC分解法： A=，b=[1/2,1,1/3,1/4,0,0,1/2,0,1,1] >> A=[1,1/2,1/3,1/4,1/5,1/6,1/7,1/8,1/9,1/10;1/2,1/3,1/4,1/5,1/6,1/7,1/8,1/9,1/10,1/11;. 1/3,1/4,1/5,1/6,1/7,1/8,1/9,1/10,1/11,1/12;1/4,1/5,1/6,1/7,1/8,1/9,1/10,1/11,1/12,1/13;. 1/5,1/6,1/7,1/8,1/9,1/10,1/11,1/12,1/13,1/14;1/6,1/7,1/8,1/9,1/10,1/11,1/12,1/13,1/14,1/15;. 1/7,1/8,1/9,1/10,1/11,1/12,1/13,1/14,1/15,1/16;1/8,1/9,1/10,1/11,1/12,1/13,1/14,1/15,1/16,1/17;. 1/9,1/10,1/11,1/12,1/13,1/14,1/15,1/16,1/17,1/18;1/10,1/11,1/12,1/13,1/14,1/15,1/16,1/17,1/18,1/19] A = Columns 1 through 9 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625 0.0588 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625 0.0588 0.0556 Column 10 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625 0.0588 0.0556 0.0526 >> b=[1/2,1,1/3,1/4,0,0,1/2,0,1,1] b = 0.5000 1.0000 0.3333 0.2500 0 0 0.5000 0 1.0000 1.0000 >> [L,D]=lu(A) L = Columns 1 through 9 1.0000 0 0 0 0 0 0 0 0 0.5000 1.0000 -0.8571 1.0000 0 0 0 0 0 0.3333 1.0000 0 0 0 0 0 0 0 0.2500 0.9000 0.5143 -0.2857 0.7500 0.2837 1.0000 0 0 0.2000 0.8000 0.7837 -0.2857 1.0000 0 0 0 0 0.1667 0.7143 0.9184 -0.1984 0.8333 -0.2579 -0.5051 1.0000 0 0.1429 0.6429 0.9796 -0.0952 0.4545 -0.2579 -0.3497 0.9643 -0.5204 0.1250 0.5833 1.0000 0 0 0 0 0 0 0.1111 0.5333 0.9974 0.0808 -0.4569 0.4444 0.1740 -0.7000 1.0000 0.1000 0.4909 0.9818 0.1469 -0.8811 1.0000 0 0 0 Column 10 0 0 0 0 0 0 1.0000 0 0 0 D = Columns 1 through 9 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0 0.0833 0.0889 0.0833 0.0762 0.0694 0.0635 0.0583 0.0539 0 0 0.0065 0.0110 0.0139 0.0156 0.0165 0.0170 0.0172 0 0 0 0.0011 0.0024 0.0034 0.0043 0.0049 0.0053 0 0 0 0 -0.0000 -0.0001 -0.0002 -0.0003 -0.0003 0 0 0 0 0 0.0000 0.0000 0.0000 0.0000 0 0 0 0 0 0 -0.0000 -0.0000 -0.0000 0 0 0 0 0 0 0 0.0000 0.0000 0 0 0 0 0 0 0 0 -0.0000 0 0 0 0 0 0 0 0 0 Column 10 0.1000 0.0500 0.0172 0.0056 -0.0004 0.0001 -0.0000 0.0000 -0.0000 -0.0000 1． (1) Gauss消去法 A=[0.4096,0.1234,0.3678,0.2943;0.2246,0.3872,0.4015,0.1129;0.3645,0.1920,. 0.3781,0.0643;0.1784,0.4001,0.2786,0.3927]; b=[0.4043,0.1550,0.4240,-0.2557] ; [m,n]=size(A); if m~=n error( 矩阵A的行数和列数必须相等 ); return; end if m~=size(b) error( b的大小必须和A的行数或A的列数相同 ); return; end if rank(A)~=rank([A,b]) error( A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解 ); return; end c=n+1; A(:,c)=b; for k=1:n-1 A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c); end x=zeros(length(b),1); x(n)=A(n,c)/A(n,n); for k=n-1:-1:1 x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k); end disp( x= ); disp(x); >>Untitled7 x= -0.1831 -1.6641 2.2189 -0.4467 (2)列主元素消去法： at long; A=[0.4096,0.1234,0.3678,0.2943;0.2246,0.3872,0.4015,0.1129;0.3645,0.1920,. 0.3781,0.0643;0.1784,0.4001,0.2786,0.3927]; b=[0.4043,0.1550,0.4240,-0.2557] ; [m,n]=size(A); if m~=n error( 矩阵A的行数和列数必须相等 ); return; end if m~=size(b) error( b的大小必须和A的行数或A的列数相同 ); return; end if rank(A)~=rank([A,b]) error( A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解 ); return; end c=n+1; A(:,c)=b; for k=1:n-1 [r,m]=max(abs(A(k:n,k))); m=m+k-1; if(A(m,k)~=0) if(m~=k) A([k m],:)=A([m k],:); end A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c); end end x=zeros(length(b),1); x(n)=A(n,c)/A(n,n); for k=n-1:-1:1 x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k); end disp( X= ); disp(x); at short; >> Untitled X= -0.183103037000047 -1.664079887224428 2.218904588362461 -0.446709634933741 结论： 高斯的计算步骤很多要进行一系列矩阵的计算，而且必须满足矩阵A的行、列相等才有唯一解，不然就会有不同的结果，而且计算时间很长，速度很慢，而追赶法计算比较快，不用进行矩阵的计算就可以得出结果