最新Gauss列主元素消去法实验.docx
- 文档编号:28693495
- 上传时间:2023-07-19
- 格式:DOCX
- 页数:32
- 大小:64.10KB
最新Gauss列主元素消去法实验.docx
《最新Gauss列主元素消去法实验.docx》由会员分享,可在线阅读,更多相关《最新Gauss列主元素消去法实验.docx(32页珍藏版)》请在冰豆网上搜索。
最新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.给定方程组
(1)
(2)
先用编写的程序计算,再将
(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.Matlab6.0及以上版本。
实验讲评:
实验成绩:
评阅教师:
200年月日
Lab06.Gauss列主元素消去法实验
第一题:
1、算法描述:
Ⅰ、Gauss消去法由书上定理5可知设Ax=b,其中A∈R
^(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、det
1
2、对于k=1,2,…,n-1
(1)按列主元
(2)如果
,则停止(det(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,1
5、
第二题:
编写用不选主元的直接三角分解法解线性方程组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)/U(1,1);
fork=2:
n
fori=k:
n
U(k,i)=A(k,i)-L(k,1:
(k-1))*U(1:
(k-1),i)
end
forj=(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);
functiony=SolveDownTriangle(L,b)
N=size(L);
n=N
(1);
fori=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
functionx=SolveUpTriangle(U,y)
N=size(U);
n=N
(1);
fori=n:
-1:
1
if(i 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 第三题: 编写用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); ifn~=m error('TherowsandcolumnsofnatrixAmustbeequal! '); return; end ifm~=nb error('ThecolumnsofAmustbeequalthelengthb! '); return; end index=1;det=1;x=zeros(n,1); fork=1: n-1 maxa=0; fori=k: n ifabs(A(i,k))>maxa maxa=abs(A(i,k));r=i; end end ifmaxa<1e-10 index=0;return; end ifr>k forj=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 fori=k+1: n m=A(i,k)/A(k,k); forj=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); ifabs(A(n,n))<1e-10 index=0;return; end fork=n: -1: 1 forj=k+1: n b(k)=b(k)-A(k,j)*x(j); end x(k)=b(k)/A(k,k); end 第四题: Ⅰ、运用Gauss列主元素消去法解上述程序计算矩阵: ⑴、计算第一个矩阵 1、计算源程序: 计算程序如下: 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); ifn~=m error('TherowsandcolumnsofnatrixAmustbeequal! '); return; end ifm~=nb error('ThecolumnsofAmustbeequalthelengthb! '); return; end index=1;det=1;x=zeros(n,1); fork=1: n-1 maxa=0; fori=k: n ifabs(A(i,k))>maxa maxa=abs(A(i,k));r=i; end end ifmaxa<1e-10 index=0;return; end ifr>k forj=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 fori=k+1: n m=A(i,k)/A(k,k); forj=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); ifabs(A(n,n))<1e-10 index=0;return; end fork=n: -1: 1 forj=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); ifn~=m error('TherowsandcolumnsofnatrixAmustbeequal! '); return; end ifm~=nb error('ThecolumnsofAmustbeequalthelengthb! '); return; end index=1;det=1;x=zeros(n,1); fork=1: n-1 maxa=0; fori=k: n ifabs(A(i,k))>maxa maxa=abs(A(i,k));r=i; end end ifmaxa<1e-10 index=0;return; end ifr>k forj=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 fori=k+1: n m=A(i,k)/A(k,k); forj=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); ifabs(A(n,n))<1e-10 index=0;return; end fork=n: -1: 1 forj=k+1: n b(k)=b(k)-A(k,j)*x(j); end x(k)=b(k)/A(k,k); end 计算结果如下: 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); ifn~=m error('TherowsandcolumnsofnatrixAmustbeequal! '); return; end ifm~=nb error('ThecolumnsofAmustbeequalthelengthb! '); return; end index=1;det=1;x=zeros(n,1); fork=1: n-1 maxa=0; fori=k: n ifabs(A(i,k))>maxa maxa=abs(A(i,k));r=i; end end ifmaxa<1e-10 index=0;return; end ifr>k forj=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 fori=k+1: n m=A(i,k)/A(k,k); forj=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); ifabs(A(n,n))<1e-10 index=0;return; end fork=n: -1: 1 forj=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 >> ②、计算修改后的程序: 计算程序如下: 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); ifn~=m error('TherowsandcolumnsofnatrixAmustbeequal! '); return; end ifm~=nb error('ThecolumnsofAmustbeequalthelengthb! '); return; end index=1;det=1;x=zeros(n,1); fork=1: n-1 maxa=0; fori=k: n ifabs(A(i,k))>maxa maxa=abs(A(i,k));r=i; end end ifmaxa<1e-10 index=0;return; end ifr>k forj=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 fori=k+1: n m=A(i,k)/A(k,k); forj=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); ifabs(A(n,n))<1e-10 index=0;return; end fork=n: -1: 1 forj=k+1: n b(k)=b(k)-A(k,j)*x(j); end x(k)=b(k)/A(k,k); end 计算结果如下: ans= -0.3543 -1.4252 1.3827 1.5669 >> Ⅱ、运用不选主元的直接三角分解法解上述程: 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); fork=2: n fori=k: n U(k,i)=A(k,i)-L(k,1: (k-1))*U(1: (k-1),i) end forj=(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); functiony=SolveDownTriangle(L,b) N=size(L); n=N (1); fori=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 functionx=SolveUpTriangle(U,y) N=size(U); n=N (1); fori=n: -1: 1 if(i 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= 3.01006.03001.9900 01.61580 000 U= 3.01006.03001.9900 01.6158-2.0696 000 L= 1.000000 0.42191.00000 0.3279-4.20061.0000 U= 3.01006.03001.9900 01.6158-2.0696 00-0.0063 ans= 1.0e+003* 1.5926 -0.6319 -0.4936 >> 2、计算修改后的程序: 计算程序如下: 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); fork=2: n fori=k: n U(k,i)=A(k,i)-L(k,1: (k-1))*U(1: (k-1),i) end forj=(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); functiony=SolveDownTriangle(L,b) N=size(L); n=N (1); fori=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 functionx=SolveUpTriangle(U,y) N=size(U); n=N (1); fori=n: -1: 1 if(i 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= 3.00006.03001.9900 01.60730 000 U= 3.00006.03001.9900 01.6073-2.0724 000 L= 1.000000 0.42331.00000 0.3300-4.23061.0000 U= 3.00006.03001.9900 01.6073-2.0724 00-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); fork=2: n fori=k: n U(k,i)=A(k,i)-L(k,1: (k-1))*U(1: (k-1),i) end forj=(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); functiony=SolveDownTriangle(L,b) N=size(L); n=N (1); fori=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 functionx=SolveUpTriangle(U,y) N=size(U); n=N (1); fori=n: -1: 1 if(i 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.000001.0000 0-0.000000 0000 0000 U= 10.0000-7.000001.0000 0-0.00006.00000 0000 0000 U= 10.0000-7.000001.0000 0-0.00006.00002.3000 0000 0000 L= 1.0e+005* 0.0000000 -0.00000.000000 0.0000-2.50000.00000 0.0000000.0000 L= 1.0e+005* 0.0000000 -0.00000.000000 0.0000-2.50000.00000 0.0000-2.400000.0000 U= 1.0e+006* 0.0000-0.000000.0000 0-0.00000.00000.0000 001.50000 0000 U= 1.0e+006* 0.0000-0.000000.0000 0
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 最新 Gauss 元素 消去 实验