MATLAB实验一解线性方程组的直接法.docx

1、MATLAB实验一解线性方程组的直接法实验报告课程名称 数值分析实验项目 解线性方程组的直接法专业班级 姓名 学号指导教师 成绩 日期 月曰1.实验目的1掌握程序的录入和matlab的使用和操作；2、 了解影响线性方程组解的精度的因素方法与问题的性态。3、 学会Matlab提供的“”的求解线性方程组。2.实验要求1按照题目要求完成实验内容；2、 写出相应的 Matlab程序；3、 给出实验结果(可以用表格展示实验结果)；4、 分析和讨论实验结果并提出可能的优化实验。5、 写出实验报告。3.实验步骤1用LU分解及列主元高斯消去法解线性方程组,z10-701、1f 、X1，z 8 -32.0999

2、9962 1X25.900001a)5-15一1X352102丿k丿 1丿输出Ax二b中系数A二LU分解的矩阵L和U，解向量x和det(A);用列主元法 的行交换次序解向量x和求det( A);比较两种方法所得结果。2、用列主高斯消元法解线性方程组 Ax二bp.016.031.99、p1、I(1)、1.274.16-1.23x=10987-4.819.34 Z丿l1丿广 3.006.031.99、行(2)、1.274.16-1.23X2=10990-4.819.34 ”W丿J丿分别输出A, b,det(A),解向量x , (1)中A的条件数。分析比较(1)、(2)的计算结果con d(A)2.

3、若令求解(A 、：A)(x，x)二b，输出向量 x和、胡2，从理论结果和实际计算两方面分析线分量的有效位数如何随 n变化，它与条件数有何关系？当 n多大时x连一位有效数字也3、线性方程组Ax = b的A和b分别为n=6 时：没有了？GAUSS列主消去法求得的xx的有效数字将每种情形的两个结果进行表格对比，如:四、 实验结果五、 讨论分析(对上述算例的计算结果进行比较分析， 主要说清matlab的算符与消去法的适用范围不同，自己补充)六、 改进实验建议(自己补充)1列主元的高斯消去法利用列主元的高斯消去法matlab程序源代码：首先建立一个gaussMethod.m的文件，用来实现列主元的消去方

4、法。fun cti on x=gaussMethod(A,b)%高斯列主元消去法，要求系数矩阵非奇异的， %n = size(A,1);if abs(det(A) a=10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2 a =10.0000-7.000001.0000-3.00002.10006.00002.00005.0000-1.00005.0000-1.00002.00001.000002.0000l =l,u=lu(a)11.0000000-0.3000-0.00001.000000.50001.0000000.20000.9600-0.80001.

5、0000u =10.0000-7.000001.000002.50005.0000-1.5000006.00002.30000005.0800b=8 5.900001 5 1b =8.00005.90005.00001.0000 y=lb y =8.00001.00008.30005.0800 x1=Ux x1 =0.0000-1.00001.00001.0000det仁 det(a)detl =-762.00012、(1)在 MATLAB 窗口：A=3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34 A =3.01006.03001.99001.2

6、7004.1600-1.23000.9870-4.81009.3400 b=1 1 1b =111x1,det1,i ndex=Gauss (A,b) x1 =1.0e+03 *1592.599624841381-631.9113762025488-493.6177247593899det1 =-0.0305in dex =1在MATLAB窗口： A=3.00 6.03 1.99;1.27 4.16 -1.23;0.990 -4.81 9.34 A =3.00006.03001.99001.27004.1600-1.23000.9900-4.81009.3400 b=1 1 1b =111x2

7、,det2,i ndex=Gauss5555(A,b)x2 =119.5273-47.1426-36.8403det2 =-0.4070 in dex =13、在 MATLAB 窗口：A=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10;b=32 23 33 31;x=Abb仁32.1 22.9 33.1 30.9;x1=Ab1A1=10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98; x2=A1bdelta_b=n orm(b-b1)/no rm(b)delta_A=norm(A-A1)/norm(A)delta_

8、x1= no rm(x-x1)/no rm(x)delta x2=no rm(x-x2)/no rm(x)con d_A=c ond(A)x =1.00001.00001.00001.0000 x1 =9.2000 -12.60004.5000-1.1000 x2 =-9.586318.3741 -3.22583.5240 delta_b =0.0033 delta_A =0.0076 delta_x1 =8.1985 delta_x2 = 10.4661 cond_A =2.9841e+033、在 MATLAB 窗口：A=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10;

9、b=32 23 33 31;x=Abb1=32.1 22.9 33.1 30.9;x1=Ab1A1=10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98; x2=A1bdelta b=norm(b-b1)/norm(b)delta_A=n orm(A-A1)/norm(A) delta_x1= no rm(x-x1)/norm(x) delta_x2=no rm(x-x2)/norm(x) con d_A=c on d(A)x =1.00001.00001.00001.0000x1 =9.2000-12.60004.5000-1.1000

10、x2 =-9.586318.3741 -3.22583.5240delta_b =0.0033delta_A =0.0076delta_x1 =8.1985delta_x2 =10.4661 cond_A =2.9841e+034、k=13for n=2:6a=hilb( n);co(n)=con d(a,i nf); end x=1:6;plot(x,co);b=zeros(k);x1 仁b;x0=b;r=b;for i=2:kx=(li nspace(1,1,i);x0(1:i,(i-1)=x;H=hilb(i);b0=H*x;b(1:i,(i-1)=b0;x1=gauss2(H,bO);

11、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-16121

12、3.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

13、.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

14、.99999968-3.20E-07071.000002282.28E-064.44E-1660.99999164-8.36E-06051.000017061.71E-05050.99998042-1.96E-05-1.11E-1651.000011821.18E-05050.99999708-2.92E-0606n=10 时:x xrn有效数字1-9.05E-104.44E-1691.000000087.76E-08070.99999835-1.65E-06061.000014911.49E-054.44E-1650.99992911-7.09E-05041.000194430.000194

15、426040.99968164-0.000318365-1.11E-1641.000307150.000307149040.99983898-0.000161019041.000035373.54E-0505n=11 时:x xrn有效数字0.999999995-5.16E-09081.0000005395.39E-07-4.44E-1660.999986041-1.40E-05051.0001558750.00015588040.999072325-0.0009277031.0032587180.00325872030.992910161-0.0070898021.0096588070.00

16、965881-2.22E-1620.991981908-0.0080181031.0037076540.00370765030.999267974-0.000732-2.22E-163n=12 时:x xrn有效数字0.999999965-3.49E-088.88E-1681.0000044094.41E-064.44E-1660.999861744-0.0001383041.0018797110.00187971030.986241983-0.013758021.0603759740.060375972.22E-1620.831939173-0.16806082.22E-1611.30397

17、33630.30397336010.643862006-0.3561381.11E-1611.2606840180.260684022.22E-1610.891666924-0.10833311.11E-1611.0195107530.0195107502n=13 时:x xrn有效数字1.0000000696.89E-088.88E-1670.999989224-1.08E-05-4.44E-1651.0004135380.00041354-2.22E-1640.993147495-0.00685252.22E-1631.0612462080.06124621-2.22E-1620.6692

18、61244-0.33073882.22E-1612.1490741311.1490741300-1.65417119-2.6541712005.118548084.11854808-1.11E-160-3.243049939-4.24304991.11E-1603.7830048962.783004900-0.051792817-1.05179281.11E-1601.1743291170.174329121.11E-161五、分析讨论:实验的数学原理很容易理解，也容易上手。把运算的结果带入原 方程组，可以发现符合的还是比较好。这说明列主元消去法计算这 类方程的有效性。当A可逆时，能够将计算进

19、行到底，列 主元法就 能确保算法的稳定，而且计算量不大。直接三角消去过程，实质上 是将A分解为两个三角矩阵的乘积 A=LU并求解Ly=b的过程。回 带过程就是求解上三角方程组 Ux=yo所以在实际的运算中，矩阵L 和U可以直接计算出，而不需要任何中间步骤，从而在计算过程中 将高斯消去法的步骤进行了进一步的简略，大大提高了运算速度， 这就是三角分解法。通过以上的计算比较，2.题方程组具有严重的病态性。当系数矩阵 有微小的变化时，wucha =-401.8918 159.5435 124.6330 ，所得 的解与原方程组的解有很大的相对误差。1题方程组中当系数矩阵A 和b有微小变化时，wucha

20、=0 0 0 0，所得的解与方程组的解没有 相对误差。所以1题方程组是良性的。用MATLA内部函数inv通过 求逆矩阵，然后通过x=inv（A）*b也可以求出方程组的解，但是没 有列主元高斯消去法具有良好的稳定性。det函数求方程组系数矩 阵的行列式时所得结果和高斯消去法和三角法所得结果相同，具有 方便快捷的优点。题四可以看出，条件数越大，有效位数越少，当 n=13时，出现有效位数为0的情况。附：高斯列主消去法源程序代码function x,det,index=Gauss(A,b)% ?o ?D? ?3 iXe ? a D?- ?aGauss? OY ?D% A -方程组矩阵% b -方程组右

21、端% x -方程组的解% det -方程组行列式% index -index=0 表示求解失败， index=1表示求解成功n,m=size(A); nb=length(b); if n=merror( The rows and columns of matrix A must be equal! ); return ;endif m=nberror( The columns of A must be equal the dimension of b! );return ;endindex=1; det=1; x=zeros(n,1);for k=1:n-1% 选主元a_max=0;for i

22、=k:nif abs(A(i,k)a_maxa_max=abs(A(i,k); r=i;endendif a_maxkfor j=k:nz=A(k,j); A(k,j)=A(r,j); A(r,j)=z;endz=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:nA(i,j)=A(i,j)-m*A(k,j);end b(i)=b(i)-m*b(k); end det=det*A(k,k);enddet=det*A(n,n);% 回代过程if abs(A(n,n)1e-20index=0;return ;endfor k=n:-1:1for j=k+1:n b(k)=b(k)-A(k,j)*x(j); end x(k)=b(k)/A(k,k);end