Gauss列主元素消去法实验

Lab06．Gauss列主元素消去法实验

【实验目的和要求】

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及解向量x。 3．编写用Gauss列主元素消去法解线性方程组Ax=b的M文件。要求输出Ax=b中矩阵A及向量b、PA=LU分解的L与U、detA及解向量x，交换顺序。

4．给定方程组

3.011.27(1) 0.9876.034.164.811.99x111.23x219.34x31

103(2) 5272.099991106501x182x25.9000011x352x41 先用编写的程序计算，再将(1)中的系数3.01改为3.00，0.987改为0.990；将(2)中的系数2.099999改为2.1，5.900001改为9.5,再用Gauss列主元素消去法解，并将两次计算的结果进行比较。

【实验仪器与软件】

1．CPU主频在1GHz以上，内存在128Mb以上的PC；

2．Matlab 6.0及以上版本。 实验讲评：

实验成绩：

评阅教师：

09级数信（1）班 姓名：马勇 学号：2009051430

200 年 月 日

Lab06．Gauss列主元素消去法实验

第一题：

1、算法描述：

Ⅰ、Gauss消去法由书上定理5可知 设Ax=b，其中A∈R

^(n

(k)0(k1,2,....,n1)，则可通过高斯消去法将Ax=b约化为等价的 (1)如果akk角形线性方程组，且计算公式为：

① 消元计算（k=1,2，….，n-1）

mikaik/akk,ik1,...,n,aijbi(k1)(k)(k)aijmikakj,i,jk1,...,n,

(k)(k)(k1)bi(n)(k)mikbk,ik1,...,n.(n)n(k)② 回带公式

xnbnxi(b/ann,(i)iji1aiixj)/aii,in1,...,2,1.(i)(i)

（2）如果A为非奇异矩阵，则可通过高斯消去法将方程组Ax=b约化方程组为上三角矩阵

以上消元和回代过程总的乘除法次数为

n3n33n2n3n33，加减法次数为

3n225n6n33

以上过程就叫高斯消去法。 Ⅱ、Gauss列主元素消去法 设Ax=b.本算法用A的具有行交换的列主元素消去法，消元结果冲掉A，乘法mik冲掉aij，计算解x冲掉常数项b，行列式存放在det中。

1、 det1

2、对于k=1,2，…，n-1 （1） 按列主元

aik,kmaxaikkin

（2） 如果aik,k0，则停止（det（A）=0）

09级数信（1）班 姓名：马勇 学号：2009051430

（3） 如果ik=k则转（4）

mikakkaijaijmik*akj换行：

nbi(biji1aij*bj)/aiimik*bk

detann*detann0bnbn/ann（4） 消元计算对于 i=k+1，…，n ① aikmik=aik/akk

② 对于j=k+1，…，n

aijaijmik*akj

③ bibimik*bk

（5） detakk*d et3、如果ann0，则计算停止（det（A）=0） 4、回代求解 （1）bnbn/ann （2）对于i=n-1，…，2,1

nbi(biji1aij*bj)/aii

5、detann*det

第二题：

编写用不选主元的直接三角分解法解线性方程组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);

09级数信（1）班 姓名：马勇 学号：2009051430

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:(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) end end

y = SolveDownTriangle(L,b); x = SolveUpTriangle(U,y);

function y=SolveDownTriangle(L,b) N=size(L); n=N(1); for i=1:n if(i>1)

s=L(i,1:(i-1))*y(1:(i-1),1); else s=0; end

y(i,1)=(b(i)-s)/L(i,i); end

function x=SolveUpTriangle(U,y) N=size(U); n=N(1); for i=n:-1:1 if(is=U(i,(i+1):n)*x((1+i):n,1); else s=0; end

x(i,1)=(y(i)-s)/U(i,i); end

第三题：

编写用Gauss列主元素消去法解线性方程组Ax=b的M文件。要求输出Ax=b中矩阵A及向量b、PA=LU分解的L与U、detA及解向量x，交换顺序。 编写M文件如下：

function [x,det,index]=gauss(A,b) [n,m]=size(A); nb=length(b); if n~=m

09级数信（1）班 姓名：马勇 学号：2009051430

error('The rows and columns of natrix A must be equal !'); return; end if m~=nb

error('The columns of A must be equal the length b !'); return; end

index=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; end end

if maxa<1e-10 index =0;return; end if r>k 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; end

for i=k+1:n m=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); end

det=det*A(n,n); if abs(A(n,n))<1e-10 index=0;return; end

for k=n:-1:1 for j=k+1:n

b(k)=b(k)-A(k,j)*x(j); end

x(k)=b(k)/A(k,k); end

09级数信（1）班 姓名：马勇 学号：2009051430

第四题：

Ⅰ、运用Gauss列主元素消去法解上述程序计算矩阵：

⑴、计算第一个矩阵 ① 、计算源程序： 计算程序如下：

function [x,det,index]=gauss(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,m]=size(A); nb=length(b); if n~=m

error('The rows and columns of natrix A must be equal !'); return; end if m~=nb

error('The columns of A must be equal the length b !'); return; end

index=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; end end

if maxa<1e-10 index =0;return; end if r>k 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; end

for i=k+1:n m=A(i,k)/A(k,k); for j=k+1:n

09级数信（1）班 姓名：马勇 学号：2009051430

A(i,j)=A(i,j)-m*A(k,j); end

b(i)=b(i)-m*b(k); end

det=det*A(k,k); end

det=det*A(n,n); if abs(A(n,n))<1e-10 index=0;return; end

for k=n:-1:1 for j=k+1:n

b(k)=b(k)-A(k,j)*x(j); end

x(k)=b(k)/A(k,k); end

运行结果如下：

ans =

1.0e+003 *

1.5926 -0.6319 -0.4936 >>

②、计算修改后的程序： 计算程序如下：

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);

09级数信（1）班 姓名：马勇 学号：2009051430

if n~=m

error('The rows and columns of natrix A must be equal !'); return; end if m~=nb

error('The columns of A must be equal the length b !'); return; end

index=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; end end

if maxa<1e-10 index =0;return; end if r>k 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; end

for i=k+1:n m=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); end

det=det*A(n,n); if abs(A(n,n))<1e-10 index=0;return; end

for k=n:-1:1 for j=k+1:n

b(k)=b(k)-A(k,j)*x(j); end

x(k)=b(k)/A(k,k); end

09级数信（1）班 姓名：马勇 学号：2009051430

计算结果如下： ans =

119.5273 -47.1426 -36.8403 >>

⑵：计算第二个矩阵： ①：计算源程序： 计算程序如下：

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; end if m~=nb

error('The columns of A must be equal the length b !'); return; end

index=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; end end

09级数信（1）班 姓名：马勇 学号：2009051430

if maxa<1e-10 index =0;return; end if r>k 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; end

for i=k+1:n m=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); end

det=det*A(n,n); if abs(A(n,n))<1e-10 index=0;return; end

for k=n:-1:1 for j=k+1:n

b(k)=b(k)-A(k,j)*x(j); end

x(k)=b(k)/A(k,k); end

计算结果如下：

ans =

0.0000 -1.0000 1.0000 1.0000

09级数信（1）班 姓名：马勇学号：2009051430

>>

②、计算修改后的程序： 计算程序如下：

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; end if m~=nb

error('The columns of A must be equal the length b !'); return; end

index=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; end end

if maxa<1e-10 index =0;return; end if r>k 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; end

for i=k+1:n m=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

09级数信（1）班 姓名：马勇 学号：2009051430

det=det*A(k,k); end

det=det*A(n,n); if abs(A(n,n))<1e-10 index=0;return; end

for k=n:-1:1 for j=k+1:n

b(k)=b(k)-A(k,j)*x(j); end

x(k)=b(k)/A(k,k); end

计算结果如下： ans =

-0.33 -1.4252 1.3827 1.5669 >>

Ⅱ、运用不选主元的直接三角分解法解上述程：

⑴ 、计算第一个矩阵： ①、 计算源程序： 计算程序如下：

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);

09级数信（1）班 姓名：马勇 学号：2009051430

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:(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) end end

y = SolveDownTriangle(L,b); x = SolveUpTriangle(U,y);

function y=SolveDownTriangle(L,b) N=size(L); n=N(1); for i=1:n if(i>1)

s=L(i,1:(i-1))*y(1:(i-1),1); else s=0; end

y(i,1)=(b(i)-s)/L(i,i); end

function x=SolveUpTriangle(U,y) N=size(U); n=N(1); for i=n:-1:1 if(is=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 =

3.0100 6.0300 1.9900 0 1.6158 0 0 0 0

09级数信（1）班 姓名：马勇 学号：2009051430

U =

3.0100 6.0300 1.9900 0 1.6158 -2.0696 0 0 0 L =

1.0000 0 0 0.4219 1.0000 0 0.3279 -4.2006 1.0000 U =

3.0100 6.0300 1.9900 0 1.6158 -2.0696 0 0 -0.0063 ans =

1.0e+003 *

1.5926 -0.6319 -0.4936 >>

②、 计算修改后的程序： 计算程序如下：

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

09级数信（1）班 姓名：马勇 学号：2009051430

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) end end

y = SolveDownTriangle(L,b); x = SolveUpTriangle(U,y);

function y=SolveDownTriangle(L,b) N=size(L); n=N(1); for i=1:n if(i>1)

s=L(i,1:(i-1))*y(1:(i-1),1); else s=0; end

y(i,1)=(b(i)-s)/L(i,i); end

function x=SolveUpTriangle(U,y) N=size(U); n=N(1); for i=n:-1:1 if(is=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 =

3.0000 6.0300 1.9900 0 1.6073 0 0 0 0 U =

09级数信（1）班 姓名：马勇 学号：2009051430

3.0000 6.0300 1.9900 0 1.6073 -2.0724 0 0 0 L =

1.0000 0 0 0.4233 1.0000 0 0.3300 -4.2306 1.0000 U =

3.0000 6.0300 1.9900 0 1.6073 -2.0724 0 0 -0.0844 ans =

119.5273 -47.1426 -36.8403 >>

⑵、计算第二个矩阵的程序： ①、计算源程序： 计算程序如下：

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,1:(k-1))*U(1:(k-1),i) end

09级数信（1）班 姓名：马勇 学号：2009051430

for j=(k+1):n

L(j,k) = (A(j,k)-L(j,1:(k-1))*U(1:(k-1),k))/U(k,k) end end

y = SolveDownTriangle(L,b); x = SolveUpTriangle(U,y);

function y=SolveDownTriangle(L,b) N=size(L); n=N(1); for i=1:n if(i>1)

s=L(i,1:(i-1))*y(1:(i-1),1); else s=0; end

y(i,1)=(b(i)-s)/L(i,i); end

function x=SolveUpTriangle(U,y) N=size(U); n=N(1); for i=n:-1:1 if(is=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 0 U =

10.0000 -7.0000 0 1.0000 0 -0.0000 6.0000 0 0 0 0 0

09级数信（1）班 姓名：马勇学号：2009051430

0 0 0 0 U =

10.0000 -7.0000 0 1.0000 0 -0.0000 6.0000 2.3000 0 0 0 0 0 0 0 0 L =

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.0000 L =

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.0000 U =

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 0 U =

1.0e+006 *

09级数信（1）班 姓名：马勇学号：2009051430

0.0000 -0.0000 0 0.0000 0 -0.0000 0.0000 0.0000 0 0 1.5000 0.5750 0 0 0 0 L =

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.0000 0.0000 U =

1.0e+006 *

0.0000 -0.0000 0 0.0000 0 -0.0000 0.0000 0.0000 0 0 1.5000 0.5750 0 0 0 0.0000 ans =

0.0000 -1.0000 1.0000 1.0000 >>

②、计算修改后的程序：计算程序如下:

09级数信（1）班 姓名：马勇学号：2009051430

function [x,L,U]=Doolittle(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 = 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:(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) end end

y = SolveDownTriangle(L,b); x = SolveUpTriangle(U,y);

function y=SolveDownTriangle(L,b) N=size(L); n=N(1); for i=1:n if(i>1)

s=L(i,1:(i-1))*y(1:(i-1),1); else s=0; end

y(i,1)=(b(i)-s)/L(i,i); end

function x=SolveUpTriangle(U,y) N=size(U); n=N(1); for i=n:-1:1 if(is=U(i,(i+1):n)*x((1+i):n,1); else s=0; end

x(i,1)=(y(i)-s)/U(i,i); end

09级数信（1）班 姓名：马勇 学号：2009051430

计算结果如下：

U =

10 -7 0 1 0 0 0 0 0 0 0 0 0 0 0 U =

10 -7 0 0 0 6 0 0 0 0 0 0 U =

10.0000 -7.0000 0 0 0 0 0 0

Warning: Divide by zero. > In hghj465 at 14 L =

1.0000 0 -0.3000 1.0000 0.5000 Inf 0.2000 0

Warning: Divide by zero. > In hghj465 at 14 L =

1.0000 0 -0.3000 1.0000 09级数信（1）班 0 1 0 0 0 0 6.0000 0 0 0 0 1.0000 0 0 0 2.3000 0 0 0 0 0 1.0000 0 0

姓名：马勇学号：2009051430 1.0000

0.5000 Inf 1.0000 0 0.2000 Inf 0 1.0000 U =

10.0000 -7.0000 0 1.0000 0 0 6.0000 2.3000 0 0 -Inf 0 0 0 0 U =

10.0000 -7.0000 0 0 0 6.0000 0 0 -Inf 0 0 0 L =

1.0000 0 0 -0.3000 1.0000 0 0.5000 Inf 1.0000 0.2000 Inf NaN U =

10.0000 -7.0000 0 0 0 6.0000 0 0 -Inf 0 0 0

Warning: Divide by zero.

> In hghj465>SolveUpTriangle at 39 In hghj465 at 18 ans = NaN NaN NaN

09级数信（1）班 0 2.3000 -Inf 0 0 0 0 1.0000 2.3000 -Inf NaN 姓名：马勇学号：2009051430 1.0000 1.0000

NaN >>

09级数信（1）班

姓名：马勇学号：2009051430