实验一-解线性方程组的直接法(1)

实验报告课程名称数值分析实验项目解线性方程组的直接法专业班级姓名学号指导教师成绩日期月日一.实验目的1、掌握程序的录入和matlab的使用和操作；2、了解影响线性方程组解的精度的因素——方法与问题的性态。3、学会Matlab提供的“\”的求解线性方程组。二.实验要求1、按照题目要求完成实验内容；2、写出相应的Matlab程序；3、给出实验结果(可以用表格展示实验结果)；4、分析和讨论实验结果并提出可能的优化实验。5、写出实验报告。三.实验步骤1、用LU分解及列主元高斯消去法解线性方程组a)15900001.582012151526099999.23107104321xxxx，输出bAx中系数LUA分解的矩阵L和U，解向量x和)det(A;用列主元法的行交换次序解向量x和求)det(A；比较两种方法所得结果。2、用列主高斯消元法解线性方程组bAx。（1）、11134.981.4987.023.116.427.199.103.601.3321xxx（2）、11134.981.4990.023.116.427.199.103.600.3321xxx分别输出)det(,,AbA，解向量x，（1）中A的条件数。分析比较（1）、（2）的计算结果3、线性方程组bAx的A和b分别为1095791068565778710A,31332332b则解Tx),1,1,1,1(.用MATLAB内部函数求)det(A和A的所有特征值和2)(Acond.若令98.99599.6989.998.585604.508.72.71.8710AA,求解bxxAA))((，输出向量x和2x，从理论结果和实际计算两方面分析线性方程组bAx解的相对误差22/xx以及A的相对误差22/AA的关系。4、希尔伯特矩阵nnijnRhH)(，其中11jihij，（1）分别对6,,3,2n计算)(nHcond，分析条件数作为n的函数如何变化。（2）令nTRx),1,,1,1(，计算xHbnn，然后用高斯消去法解线性方程组nnbxH求出x，计算剩余向量xHbrnnn以及xxx。分析当n增加时解x分量的有效位数如何随n变化，它与条件数有何关系？当n多大时x连一位有效数字也没有了？将每种情形的两个结果进行表格对比，如：n=6时：GAUSS列主消去法求得的xx的有效数字四、实验结果五、讨论分析（对上述算例的计算结果进行比较分析，主要说清matlab的算符与消去法的适用范围不同，自己补充）六、改进实验建议（自己补充）1.列主元的高斯消去法利用列主元的高斯消去法matlab程序源代码：首先建立一个gaussMethod.m的文件，用来实现列主元的消去方法。functionx=gaussMethod(A,b)%高斯列主元消去法，要求系数矩阵非奇异的，%n=size(A,1);ifabs(det(A))=1e-8error('系数矩阵是奇异的');return;end%fork=1:nak=max(abs(A(k:n,k)));index=find(A(:,k)==ak);iflength(index)==0index=find(A(:,k)==-ak);end%交换列主元temp=A(index,:);A(index,:)=A(k,:);A(k,:)=temp;temp=b(index);b(index)=b(k);b(k)=temp;%消元过程fori=k+1:nm=A(i,k)/A(k,k);%消除列元素A(i,k+1:n)=A(i,k+1:n)-m*A(k,k+1:n);b(i)=b(i)-m*b(k);endend%回代过程x(n)=b(n)/A(n,n);fork=n-1:-1:1;x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)')/A(k,k);endx=x';end然后调用gaussMethod函数，来实现列主元的高斯消去法。在命令框中输入下列命令：输出结果如下：利用LU分解法及matlab程序源代码：function[L,U]=myLU(A)%实现对矩阵A的LU分解，L为下三角矩阵A[n,n]=size(A);L=zeros(n,n);U=zeros(n,n);fori=1:nL(i,i)=1;endfork=1:nforj=k:nU(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j)');endfori=k+1:nL(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)'))/U(k,k);endend在命令框中输入下列命令：a=[10-701;-32.09999962;5-15-1;2102]a=10.0000-7.000001.0000-3.00002.10006.00002.00005.0000-1.00005.0000-1.00002.00001.000002.0000[l,u]=lu(a)l=1.0000000-0.3000-0.00001.000000.50001.0000000.20000.9600-0.80001.0000u=10.0000-7.000001.000002.50005.0000-1.5000006.00002.30000005.0800b=[85.90000151]'b=8.00005.90005.00001.0000y=l\by=8.00001.00008.30005.0800x1=U\xx1=0.0000-1.00001.00001.0000det1=det(a)det1=-762.00012、（1）在MATLAB窗口：A=[3.016.031.99;1.274.16-1.23;0.987-4.819.34]A=3.01006.03001.99001.27004.1600-1.23000.9870-4.81009.3400b=[111]'b=111[x1,det1,index]=Gauss(A,b)x1=1.0e+03*1592.599624841381-631.9113762025488-493.6177247593899det1=-0.0305index=1(2)在MATLAB窗口：A=[3.006.031.99;1.274.16-1.23;0.990-4.819.34]A=3.00006.03001.99001.27004.1600-1.23000.9900-4.81009.3400b=[111]'b=111[x2,det2,index]=Gauss5555(A,b)x2=119.5273-47.1426-36.8403det2=-0.4070index=13、在MATLAB窗口：A=[10787;7565;86109;75910];b=[32233331]';x=A\bb1=[32.122.933.130.9]';x1=A\b1A1=[1078.17.2;7.085.0465;85.989.899;6.99599.98];x2=A1\bdelta_b=norm(b-b1)/norm(b)delta_A=norm(A-A1)/norm(A)delta_x1=norm(x-x1)/norm(x)delta_x2=norm(x-x2)/norm(x)cond_A=cond(A)x=1.00001.00001.00001.0000x1=9.2000-12.60004.5000-1.1000x2=-9.586318.3741-3.22583.5240delta_b=0.0033delta_A=0.0076delta_x1=8.1985delta_x2=10.4661cond_A=2.9841e+033、在MATLAB窗口：A=[10787;7565;86109;75910];b=[32233331]';x=A\bb1=[32.122.933.130.9]';x1=A\b1A1=[1078.17.2;7.085.0465;85.989.899;6.99599.98];x2=A1\bdelta_b=norm(b-b1)/norm(b)delta_A=norm(A-A1)/norm(A)delta_x1=norm(x-x1)/norm(x)delta_x2=norm(x-x2)/norm(x)cond_A=cond(A)x=1.00001.00001.00001.0000x1=9.2000-12.60004.5000-1.1000x2=-9.586318.3741-3.22583.5240delta_b=0.0033delta_A=0.0076delta_x1=8.1985delta_x2=10.4661cond_A=2.9841e+034、k=13forn=2:6a=hilb(n);co(n)=cond(a,inf);endx=1:6;plot(x,co);b=zeros(k);x11=b;x0=b;r=b;fori=2:kx=(linspace(1,1,i))';x0(1:i,(i-1))=x;H=hilb(i);b0=H*x;b(1:i,(i-1))=b0;x1=gauss2(H,b0);r(1:i,(i-1))=b0-H*x1;x11(1:i,(i-1))=x1;enddx=x11-x0;结果如下：co=[0,27,748,28375,943656,29070279]可见，条件数随着n的增大，急剧增加，如图1所示。将（2）求得的结果（dx，x11）整理得n=2时：xΔxrn有效数字14.44E-160161-6.66E-16015n=3时：xΔxrn有效数字1-1.33E-15-2.22E-161519.55E-150141-9.99E-15014n=4时：xΔxrn有效数字1-2.35E-1401412.56E-130131-6.11E-131.11E-161213.96E-13013n=5时：xΔxrn有效数字1-1.21E-144.44E-161416.97E-1401311.18E-132.22E-16131-6.17E-1301214.59E-13-1.11E-1613n=6时：xΔxrn有效数字1-9.28E-1301212.67E-110111-1.82E-1001014.75E-100101-5.26E-100912.08E-10010n=7时：xΔxrn有效数字1-9.26E-1201113.71E-100101-3.59E-092.22E-1691.000000011.40E-08080.99999997-2.57E-08081.000000022.22E-081.11E-1680.99999999-7.30E-0908n=8时：xΔxrn有效数字1-2.82E-114.44E-161111.53E-09090.99999998-2.01E-082.22E-1681.000000111.09E-07070.9999997-2.95E-07-2.22E-1671.000000424.19E-07-1.11E-1670.9999997-2.99E-071.11E-1671.000000088.44E-0807n=9时：xΔxrn有效数字1-2.75E-104.44E-16101.000000021.90E-08080.99999968-3.20E-07071.000002282.28E-064.44E-1660.99999164-8.36E-06051.000017061.71E-05050.99998042-1.