《线性方程组及MATLAB应用》由会员分享,可在线阅读,更多相关《线性方程组及MATLAB应用(16页珍藏版)》请在金锄头文库上搜索。
1、数值实验 线性方程组与MATLAB应用 王1实验目的:理解矩阵的范数与条件数。 实验内容:已知矩阵 求,和。 解:编写了一个M文件来求矩阵A的范数与条件数: test3_1.m如下: A=1 1 1 1;-1 1 -1 1;-1 -1 1 1;1 -1 -1 1; norm(A,1) norm(A,2) norm(A,inf) cond(A,2) 计算结果依次是: 4 2 4 1.00002实验目的:研究高斯消去法的数值稳定性(出现小主元)。 实验内容:设方程组,其中两个矩阵如下,分别对以上两个方程组(1),(2),(1) 计算矩阵的条件数,判断系数矩阵是良态的还是病态的?解: 本题编写了一个
2、test3_21的M文件如下:A1=0.3*1e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1;A2=10 -7 0 1;-3 2.099999999999999 6 2;5 -1 5 -1;0 1 0 2;cond(A1)cond(A2)求得两个矩阵的条件数分别为68.4296和8.9939,易知这矩阵A1的条件数远远大于1,而矩阵A2的条件数刚大于1,故这,矩阵A1为病态矩阵,矩阵A2为良态矩阵。(2)用列主元消去法求得L和U及解向量; 解:本题利用Matlab的列主元三角分解函数lu();具体求解如下: A1=0.3*1e-15 59
3、.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1; A2=10 -7 0 1;-3 2.099999999999999 6 2;5 -1 5 -1;0 1 0 2; b1=59.17;46.78;1;2; b2=8;5.0000000000001;5;1; L1,U1=lu(A1) L1 = 0.0000 1.0000 0 0 0.4724 -0.1755 1.0000 0 1.0000 0 0 0 0.0893 0.0202 -0.1738 1.0000U1 = 11.2000 9.0000 5.0000 2.0000 0 59.1400 3.0000
4、 1.0000 0 0 -2.8354 1.2307 0 0 0 1.0151 L2,U2=lu(A2)L2 =1.0000 0 0 0 -0.3000 -0.0000 1.0000 0 0.5000 1.0000 0 0 0 0.4000 -0.3333 1.0000U2 =10.0000 -7.0000 0 1.0000 0 2.5000 5.0000 -1.5000 0 0 6.0000 2.3000 0 0 0 3.3667 y1=L1b1; x1=U1y1x1 =3.8457 1.6095 -15.4761 10.4113 y2=L2b2; x2=U2y2x2 =0.1337 -0.
5、8218 0.8842 0.9109用不选主元的高斯消去法求得L和U及解向量;解:编写一个LU_Fact的M文件储存不选主元的LU分解法然后调用求解。具体实现如下LU_Fact.m如下所示:function x,L,U,index=LU_Fact(A,b)%普通LU分解。A为要分解的矩阵,b为方程组右端项;%x为方程组的解,L为单位下三角阵,U为上三角阵,index为指示变量n,m=size(A); nb=length(b);if n = m error(矩阵A的行列必须相等!);return;endif m = nb error(矩阵A的行数必须和b的行数相同!); return;endL=
6、eye(n);U=zeros(n);index=1;x=zeros(n,1);y=zeros(n,1);for k=1:n for j=k:n z=0; for q=1:k-1 z=z+L(k,q)*U(q,j); end U(k,j)=A(k,j)-z; end if abs(U(k,k) x,L,U,index=LU_Fact(A1,b1)x = 0 0 0 0L = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1U = 0.0000 59.1400 3.0000 1.0000 0 0 0 0 0 0 0 0 0 0 0 0index =0对于方程组2的求解过程: x,L,
7、U,index=LU_Fact(A2,b2)x = 0 0 0 0L = 1.0000 0 0 0 -0.3000 1.0000 0 0 0.5000 0 1.0000 0 0 0 0 1.0000U = 10.0000 -7.0000 0 1.0000 0 -0.0000 6.0000 2.3000 0 0 0 0 0 0 0 0index = 0(4) 观察小主元并分析对计算结果的影响.总结:通过以上的求解过程容易发现,对于小主元的方程组求解问题,小主元导致计算的中间结果数量级增长,使得方程组无法求解或求得的解严重偏离了方程组的精确解。所以,对于小主元的方程组求解问题,应用选主元的LU分解
8、法更为可靠。3实验目的:研究线性方程组直接法的时间复杂性.实验内容:分别用高斯消去法、LU分解法、追赶法、平方根法解方程组,取,比较所用时间.解:利用tic和toc计算时间,然后调用不同方法解方程组。具体如下:(1).高斯消去法解方程组,并计算时间如下所示:ticA=spdiags(ones(500,1),4*ones(500,1),ones(500,1),-1,0,1,500,500);b=ones(500,1);x,detA,index=Gauss(A,b)toc输出:Elapsed time is 45.369000 seconds. 这个运行过多次,时间较长。(2).LU分解法解方程组
9、,并计算时间如下所示:ticA=spdiags(ones(500,1),4*ones(500,1),ones(500,1),-1,0,1,500,500);b=ones(500,1);x,L,U,index=LU_Fact(A,b)toc输出:Elapsed time is 1.929000 seconds.(3).追赶法解方程组,并计算时间如下所示:对于追赶法解方程组的问题,编写了一个追赶法算法的M文件Persue.m如下:ticA=spdiags(ones(500,1),4*ones(500,1),ones(500,1),-1,0,1,500,500);b=ones(500,1);L=ey
10、e(500);U=zeros(500);U(1,1)=A(1,1);for i=2:500 U(i-1,i)=A(i-1,i); L(i,i-1)=A(i,i-1)/U(i-1,i-1); U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i);end%P=sparse(A); Q=sparse(L); R=sparse(U); 此行目的是输出,看上述计算是否有错y=zeros(500,1);x=zeros(500,1);%追的过程如下y(1)=b(1);for i=2:500 y(i)=b(i)-L(i,i-1)*y(i-1);end%赶的过程如下x(500)=y(500)/U(500,500);for i=499:-1:1 x(i)=(y(i)-U(i,i+1)*x(i+1)/U(i,i);end toc运行如下:Elapsed time is 0.013000 seconds.(4).平方根法解方程组,并计算时间如下所示:ticA=spdiags(ones(500,1),4*ones(500,1),ones(500,1),-1,0,1,500,500);b=ones(500,1);L=chol(A);y=inv(L)b;x=Ly;toc输