最新Gauss列主元素消去法实验.docx

1、最新Gauss列主元素消去法实验Lab06Gauss列主元素消去法实验【实验目的和要求】1使学生深入理解并掌握Gauss消去法和Gauss列主元素消去法步骤；2通过对Gauss消去法和Gauss列主元素消去法的程序设计，以提高学生程序设计的能力；3对具体问题，分别用Gauss消去法和Gauss列主元素消去法求解。通过对结果的分析比较，使学生感受Gauss列主元素消去法优点。【实验内容】1根据Matlab语言特点，描述Gauss消去法和Gauss列主元素消去法步骤。2编写用不选主元的直接三角分解法解线性方程组Ax=b的M文件。要求输出Ax=b中矩阵A及向量b，A=LU分解的L与U，detA及解向

2、量x。3编写用Gauss列主元素消去法解线性方程组Ax=b的M文件。要求输出Ax=b中矩阵A及向量b、PA=LU分解的L与U、detA及解向量x，交换顺序。4给定方程组(1) (2) 先用编写的程序计算，再将(1)中的系数3.01改为3.00，0.987改为0.990；将(2)中的系数2.099999改为2.1，5.900001改为9.5,再用Gauss列主元素消去法解，并将两次计算的结果进行比较。【实验仪器与软件】 1CPU主频在1GHz以上，内存在128Mb以上的PC；2Matlab 6.0及以上版本。实验讲评：实验成绩： 评阅教师：200 年 月 日 Lab06Gauss列主元素消去法实

3、验第一题：1、算法描述：、Gauss消去法由书上定理5可知 设Ax=b，其中AR(n*n)(1)如果，则可通过高斯消去法将Ax=b约化为等价的角形线性方程组，且计算公式为：1消元计算（k=1,2，.，n-1）2回带公式（2）如果A为非奇异矩阵，则可通过高斯消去法将方程组Ax=b约化方程组为上三角矩阵以上消元和回代过程总的乘除法次数为，加减法次数为以上过程就叫高斯消去法。、Gauss列主元素消去法 设Ax=b.本算法用A的具有行交换的列主元素消去法，消元结果冲掉A，乘法冲掉，计算解x冲掉常数项b，行列式存放在det中。1、det12、对于k=1,2，n-1（1）按列主元（2）如果，则停止（det

4、（A）=0）（3）如果=k则转（4）换行：（4）消元计算对于 i=k+1，n =/ 对于j=k+1，n （5） 3、如果，则计算停止（det（A）=0）4、回代求解（1）（2）对于i=n-1，2,15、第二题：编写用不选主元的直接三角分解法解线性方程组Ax=b的M文件。要求输出Ax=b中矩阵A及向量b，A=LU分解的L与U，detA及解向量x。编写M文件如下：function x,L,U=Doolittle(A,b)N = size(A);n = N(1);L = eye(n,n); U = zeros(n,n);U(1,1:n) = A(1,1:n); L(1:n,1) = A(1:n,1)

5、/U(1,1);for k=2:n for i=k:n U(k,i) = A(k,i)-L(k,1:(k-1)*U(1:(k-1),i) end for j=(k+1):n L(j,k) = (A(j,k)-L(j,1:(k-1)*U(1:(k-1),k)/U(k,k) endendy = SolveDownTriangle(L,b);x = SolveUpTriangle(U,y);function y=SolveDownTriangle(L,b)N=size(L);n=N(1);for i=1:n if(i1) s=L(i,1:(i-1)*y(1:(i-1),1); else s=0; e

6、nd y(i,1)=(b(i)-s)/L(i,i);endfunction x=SolveUpTriangle(U,y) N=size(U);n=N(1);for i=n:-1:1 if(imaxa maxa=abs(A(i,k);r=i; endendif maxak for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z;det=-det;endfor i=k+1:nm=A(i,k)/A(k,k);for j=k+1:n A(i,j)=A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(

7、k);end det=det*A(k,k); enddet=det*A(n,n);if abs(A(n,n)maxa maxa=abs(A(i,k);r=i; endendif maxak for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z;det=-det;endfor i=k+1:nm=A(i,k)/A(k,k);for j=k+1:n A(i,j)=A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k);end det=det*A(k,k); enddet=det*A(n,n);

8、if abs(A(n,n)、计算修改后的程序：计算程序如下：function x,det,index=gauss(A,b) A=3.00,6.03,1.99;1.27,4.16,-1.23;0.990,-4.81,9.34;b=1;1;1;n,m=size(A);nb=length(b);if n=m error(The rows and columns of natrix A must be equal !); return;endif m=nb error(The columns of A must be equal the length b !); return;endindex=1;d

9、et=1;x=zeros(n,1);for k=1:n-1 maxa=0; for i=k:n if abs(A(i,k)maxa maxa=abs(A(i,k);r=i; endendif maxak for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z;det=-det;endfor i=k+1:nm=A(i,k)/A(k,k);for j=k+1:n A(i,j)=A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k);end det=det*A(k,k); enddet=det*

10、A(n,n);if abs(A(n,n)：计算第二个矩阵：计算源程序：计算程序如下：function x,det,index=gauss(A,b) A=10,-7,0,1;-3,2.09999,6,2;5,-1,5,-1;2,1,0,2;b=8;5.900001;5;1;n,m=size(A);nb=length(b);if n=m error(The rows and columns of natrix A must be equal !); return;endif m=nb error(The columns of A must be equal the length b !); ret

11、urn;endindex=1;det=1;x=zeros(n,1);for k=1:n-1 maxa=0; for i=k:n if abs(A(i,k)maxa maxa=abs(A(i,k);r=i; endendif maxak for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z;det=-det;endfor i=k+1:nm=A(i,k)/A(k,k);for j=k+1:n A(i,j)=A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k);end det=det*A(k

12、,k); enddet=det*A(n,n);if abs(A(n,n)、计算修改后的程序：计算程序如下：function x,det,index=gauss(A,b) A=10,-7,0,1;-3,2.1,6,2;5,-1,5,-1;2,1,0,2;b=8;9.5;5;1;n,m=size(A);nb=length(b);if n=m error(The rows and columns of natrix A must be equal !); return;endif m=nb error(The columns of A must be equal the length b !); r

13、eturn;endindex=1;det=1;x=zeros(n,1);for k=1:n-1 maxa=0; for i=k:n if abs(A(i,k)maxa maxa=abs(A(i,k);r=i; endendif maxak for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z;det=-det;endfor i=k+1:nm=A(i,k)/A(k,k);for j=k+1:n A(i,j)=A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k);end det=det*A

14、(k,k); enddet=det*A(n,n);if abs(A(n,n)、运用不选主元的直接三角分解法解上述程：1、计算第一个矩阵：1、计算源程序：计算程序如下：function x,L,U=Doolittle(A,b)A=3.01,6.03,1.99;1.27,4.16,-1.23;0.987,-4.81,9.34;b=1,1,1;N = size(A);n = N(1);L = eye(n,n); U = zeros(n,n);U(1,1:n) = A(1,1:n); L(1:n,1) = A(1:n,1)/U(1,1);for k=2:n for i=k:n U(k,i) = A(k

15、,i)-L(k,1:(k-1)*U(1:(k-1),i) end for j=(k+1):n L(j,k) = (A(j,k)-L(j,1:(k-1)*U(1:(k-1),k)/U(k,k) endendy = SolveDownTriangle(L,b);x = SolveUpTriangle(U,y);function y=SolveDownTriangle(L,b)N=size(L);n=N(1);for i=1:n if(i1) s=L(i,1:(i-1)*y(1:(i-1),1); else s=0; end y(i,1)=(b(i)-s)/L(i,i);endfunction x=

16、SolveUpTriangle(U,y) N=size(U);n=N(1);for i=n:-1:1 if(i2、计算修改后的程序：计算程序如下：function x,L,U=Doolittle(A,b)A=3.00,6.03,1.99;1.27,4.16,-1.23;0.990,-4.81,9.34;b=1,1,1;N = size(A);n = N(1);L = eye(n,n); U = zeros(n,n);U(1,1:n) = A(1,1:n); L(1:n,1) = A(1:n,1)/U(1,1);for k=2:n for i=k:n U(k,i) = A(k,i)-L(k,1:

17、(k-1)*U(1:(k-1),i) end for j=(k+1):n L(j,k) = (A(j,k)-L(j,1:(k-1)*U(1:(k-1),k)/U(k,k) endendy = SolveDownTriangle(L,b);x = SolveUpTriangle(U,y);function y=SolveDownTriangle(L,b)N=size(L);n=N(1);for i=1:n if(i1) s=L(i,1:(i-1)*y(1:(i-1),1); else s=0; end y(i,1)=(b(i)-s)/L(i,i);endfunction x=SolveUpTri

18、angle(U,y) N=size(U);n=N(1);for i=n:-1:1 if(i、计算第二个矩阵的程序：、计算源程序：计算程序如下：function x,L,U=Doolittle(A,b)A=10,-7,0,1;-3,2.09999,6,2;5,-1,5,-1;2,1,0,2;b=8,5.900001,5,1;N = size(A);n = N(1);L = eye(n,n); U = zeros(n,n);U(1,1:n) = A(1,1:n); L(1:n,1) = A(1:n,1)/U(1,1);for k=2:n for i=k:n U(k,i) = A(k,i)-L(k,

19、1:(k-1)*U(1:(k-1),i) end for j=(k+1):n L(j,k) = (A(j,k)-L(j,1:(k-1)*U(1:(k-1),k)/U(k,k) endendy = SolveDownTriangle(L,b);x = SolveUpTriangle(U,y);function y=SolveDownTriangle(L,b)N=size(L);n=N(1);for i=1:n if(i1) s=L(i,1:(i-1)*y(1:(i-1),1); else s=0; end y(i,1)=(b(i)-s)/L(i,i);endfunction x=SolveUpT

20、riangle(U,y) N=size(U);n=N(1);for i=n:-1:1 if(in) s=U(i,(i+1):n)*x(1+i):n,1); else s=0; end x(i,1)=(y(i)-s)/U(i,i);end计算结果如下：U = 10.0000 -7.0000 0 1.0000 0 -0.0000 0 0 0 0 0 0 0 0 0 0U = 10.0000 -7.0000 0 1.0000 0 -0.0000 6.0000 0 0 0 0 0 0 0 0 0U = 10.0000 -7.0000 0 1.0000 0 -0.0000 6.0000 2.3000 0

21、 0 0 0 0 0 0 0L = 1.0e+005 * 0.0000 0 0 0 -0.0000 0.0000 0 0 0.0000 -2.5000 0.0000 0 0.0000 0 0 0.0000L = 1.0e+005 * 0.0000 0 0 0 -0.0000 0.0000 0 0 0.0000 -2.5000 0.0000 0 0.0000 -2.4000 0 0.0000U = 1.0e+006 * 0.0000 -0.0000 0 0.0000 0 -0.0000 0.0000 0.0000 0 0 1.5000 0 0 0 0 0U = 1.0e+006 * 0.0000 -0.0000 0 0.0000 0