《matlab解线性方程组 线性方程组及MATLAB应用.doc》由会员分享，可在线阅读，更多相关《matlab解线性方程组 线性方程组及MATLAB应用.doc（35页珍藏版）》请在金锄头文库上搜索。
1、matlab解线性方程组 线性方程组及MATLAB应用 数值实验 线性方程组与MATLAB应用王1实验目的：理解矩阵的范数与条件数。 实验内容：已知矩阵111-11-1A=-1-111-1-111求A1，A2，A和cond2(A)。 11解：编写了一个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实验目的：研究高斯消去法的数值稳定性（出现小主元）。实验内容：设方程组Ax=b，其中
2、两个矩阵如下，分别对以上两个方程组0.310-155.291（1）A1=11.21159.1746.78-6.130-12 ，b1=1952221159.143-710-32.09999999999999（2）A2=5-1100185.9000000000000162 ，b2=55-1102（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的？解： 本题编写了一个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
3、-1 5 -1;0 1 0 2; cond(A1) cond(A2)求得两个矩阵的条件数分别为68.4296和8.9939，易知这矩阵A1的条件数远远大于1，而矩阵A2的条件数刚大于1，故这,矩阵A1为病态矩阵，矩阵A2为良态矩阵。 （2）用列主元消去法求得L和U及解向量x1,x2R；解：本题利用Matlab的列主元三角分解函数lu()；具体求解如下： 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;4 b1=
4、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.0000 U1 = 11.2000 9.0000 5.0000 2.0000 0 59.1400 3.0000 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.000
5、0 0 0 0 0.4000 -0.3333 1.0000 U2 =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=U1y1 x1 =3.8457 1.6095 -15.4761 10.4113 y2=L2b2; x2=U2y2 x2 =0.1337 -0.8218 0.8842 0.9109%1,x%2R4x用不选主元的高斯消去法求得L和U及解向量；解：编写一个LU_Fact的M文件储存不选主元的LU分解法然后调用求解。具体实现如下 LU_Fact.m如下所
6、示：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 = merror(矩阵A的行列必须相等！);return; endif m = nberror(矩阵A的行数必须和b的行数相同！); return; endL=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-1z=z
7、+L(k,q)*U(q,j); endU(k,j)=A(k,j)-z; endif abs(U(k,k)for i=k+1:n z=0;for q=1:k-1z=z+L(i,q)*U(q,k); endL(i,k)=(A(i,k)-z)/U(k,k); end endy(1)=b(1); for k=2:n z=0;for j=1:k-1z=z+L(k,j)*y(j); endy(k)=b(k)-z; endx(n)=y(n)/U(n,n); for k=n-1:-1:1 z=0;for j=k+1:nz=z+U(k,j)*x(j); endx(k)=(y(k)-z)/U(k,k); end求
8、解过程及结果如下： 对于方程组1的求解过程： 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 0 index =0对于方程组2的求解过程： x,L,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.0000 U = 10.0000 -7.0000 0
9、1.0000 0 -0.0000 6.0000 2.3000 0 0 0 0 0 0 0 0 index = 0（4）观察小主元并分析对计算结果的影响.总结：通过以上的求解过程容易发现，对于小主元的方程组求解问题，小主元导致计算的中间结果数量级增长，使得方程组无法求解或求得的解严重偏离了方程组的精确解。所以，对于小主元的方程组求解问题，应用选主元的LU分解法更为可靠。3实验目的：研究线性方程组直接法的时间复杂性.实验内容：分别用高斯消去法、LU分解法、追赶法、平方根法解方程组41x1114Ox12=，取n=500，比较所用时间. OO1MM14xn1解：利用tic和toc计算时间，然后调用不同
10、方法解方程组。具体如下：(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分解法解方程组，并计算时间如下所示： ticA=spdiags(ones(500,1),4*ones(500,1),ones(500,1),-1,0,1,500,500); b=ones(500
11、,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=eye(500); U=zeros(500); U(1,1)=A(1,1); for i=2:500U(i-1,i)=A(i-1,i);L(i,i-1)=A(i,i-1)/
12、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:500y(i)=b(i)-L(i,i-1)*y(i-1); end%赶的过程如下x(500)=y(500)/U(500,500); for i=499:-1:1x(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=i
