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

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

2、量和求；比较两种方法所得结果。2、用列主高斯消元法解线性方程组。（1）、（2）、分别输出，解向量，（1）中的条件数。分析比较（1）、（2）的计算结果3、线性方程组的和分别为,则解. 用MATLAB内部函数求和的所有特征值和. 若令,求解，输出向量和，从理论结果和实际计算两方面分析线性方程组解的相对误差以及的相对误差的关系。4、 希尔伯特矩阵，其中，（1）分别对计算，分析条件数作为的函数如何变化。（2）令，计算，然后用高斯消去法解线性方程组求出，计算剩余向量以及。分析当增加时解分量的有效位数如何随变化，它与条件数有何关系？当多大时连一位有效数字也没有了？将每种情形的两个结果进行表格对比，如：n=

3、6时：GAUSS列主消去法求得的的有效数字四、实验结果五、讨论分析（对上述算例的计算结果进行比较分析，主要说清matlab的算符与消去法的适用范围不同，自己补充）六、改进实验建议（自己补充）1.列主元的高斯消去法 利用列主元的高斯消去法matlab程序源代码：首先建立一个gaussMethod.m的文件，用来实现列主元的消去方法。function x=gaussMethod(A,b)%高斯列主元消去法，要求系数矩阵非奇异的， %n = size(A,1);if abs(det(A)<= 1e-8 error('系数矩阵是奇异的'); return;end%for k=1:

4、n ak = max(abs(A(k:n,k); index = find(A(:,k)=ak); if length(index) = 0 index = 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; %消元过程 for i=k+1:n m=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

5、(k); endend %回代过程 x(n)=b(n)/A(n,n); for k=n-1:-1:1; x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)')/A(k,k); end x=x'end 然后调用gaussMethod函数，来实现列主元的高斯消去法。在命令框中输入下列命令：输出结果如下：利用LU分解法及matlab程序源代码：function L,U=myLU(A) %实现对矩阵A的LU分解，L为下三角矩阵 An,n=size(A); L=zeros(n,n); U=zeros(n,n); for i=1:n L(i,i)=1; end for k=1:n

6、 for j=k:n U(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j)'); end for i=k+1:n L(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)')/U(k,k); end end 在命令框中输入下列命令：>> a=10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2a = 10.0000 -7.0000 0 1.0000 -3.0000 2.1000 6.0000 2.0000 5.0000 -1.0000 5.0000 -1.0000 2.0000

7、1.0000 0 2.0000>> l,u=lu(a)l = 1.0000 0 0 0 -0.3000 -0.0000 1.0000 0 0.5000 1.0000 0 0 0.2000 0.9600 -0.8000 1.0000u = 10.0000 -7.0000 0 1.0000 0 2.5000 5.0000 -1.5000 0 0 6.0000 2.3000 0 0 0 5.0800>> b=8 5.900001 5 1'b = 8.0000 5.9000 5.00001.0000>> y=lby =8.00001.00008.30005.

8、0800 >> x1=Uxx1 =0.0000-1.00001.00001.0000>>det1= det(a)det1 =-762.00012、（1）在MATLAB窗口：>>A=3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34A = 3.0100 6.0300 1.9900 1.2700 4.1600 -1.2300 0.9870 -4.8100 9.3400>> b=1 1 1'b = 1 1 1>>x1,det1,index=Gauss (A,b)x1 = 1.0e+03 *

9、 1592.599624841381 -631.9113762025488 -493.6177247593899det1 = -0.0305index = 1 (2) 在MATLAB窗口：>> A=3.00 6.03 1.99;1.27 4.16 -1.23;0.990 -4.81 9.34A = 3.0000 6.0300 1.9900 1.2700 4.1600 -1.2300 0.9900 -4.8100 9.3400>> b=1 1 1'b = 1 1 1>>x2,det2,index=Gauss5555(A,b)x2 = 119.5273

10、-47.1426 -36.8403det2 = -0.4070index = 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=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=norm(A-A1)/norm(A)delta_x1=norm(x-x1)/norm(x)delta_x2=norm

11、(x-x2)/norm(x)cond_A=cond(A)x = 1.0000 1.0000 1.0000 1.0000x1 = 9.2000 -12.6000 4.5000 -1.1000x2 = -9.5863 18.3741 -3.2258 3.5240delta_b = 0.0033delta_A = 0.0076delta_x1 = 8.1985delta_x2 = 10.4661cond_A = 2.9841e+033、在MATLAB窗口： A=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10; b=32 23 33 31' x=Abb1=32.1 22.

12、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=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.0000 1.0000 1.0000 1.0000x1 = 9.2000 -12.6000 4.5000 -1.1000x2 = -9.5863 18.3741 -3.2258

13、 3.5240delta_b = 0.0033delta_A = 0.0076delta_x1 = 8.1985delta_x2 = 10.4661cond_A = 2.9841e+034、k=13for n=2:6 a=hilb(n); co(n)=cond(a,inf);endx=1:6;plot(x,co);b=zeros(k);x11=b;x0=b;r=b;for i=2:k x=(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

14、-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时：xrn有效数字14.44E-160161-6.66E-16015n=3时：xrn有效数字1-1.33E-15-2.22E-161519.55E-150141-9.99E-15014n=4时：xrn有效数字1-2.35E-1401412.56E-130131-6.11E-131.11E-161213.96E-13013n=5时：xrn有效数字1-1

15、.21E-144.44E-161416.97E-1401311.18E-132.22E-16131-6.17E-1301214.59E-13-1.11E-1613n=6时：xrn有效数字1-9.28E-1301212.67E-110111-1.82E-1001014.75E-100101-5.26E-100912.08E-10010n=7时：xrn有效数字1-9.26E-1201113.71E-100101-3.59E-092.22E-1691.000000011.40E-08080.99999997-2.57E-08081.000000022.22E-081.11E-1680.9999999

16、9-7.30E-0908n=8时：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时：xrn有效数字1-2.75E-104.44E-16101.000000021.90E-08080.99999968-3.20E-07071.000002282.28E-064.4

17、4E-1660.99999164-8.36E-06051.000017061.71E-05050.99998042-1.96E-05-1.11E-1651.000011821.18E-05050.99999708-2.92E-0606n=10时：xrn有效数字1-9.05E-104.44E-1691.000000087.76E-08070.99999835-1.65E-06061.000014911.49E-054.44E-1650.99992911-7.09E-05041.000194430.000194426040.99968164-0.000318365-1.11E-1641.00030

18、7150.000307149040.99983898-0.000161019041.000035373.54E-0505n=11时：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.00965881-2.22E-1620.991981908-0.0080181031.003707

19、6540.00370765030.999267974-0.000732-2.22E-163n=12时：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.3039733630.30397336010.643862006-0.3561381.11E-1611.260

20、6840180.260684022.22E-1610.891666924-0.10833311.11E-1611.0195107530.0195107502n=13时：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.669261244-0.33073882.22E-1612.1490741311.1490741300-1.654

21、17119-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可逆时，能够将计算进行到底，列主元法就能确保算法的稳定，而且计算量不大。直接三角消去过程，实质上是将A分解为两个三角矩阵的乘积A

22、=LU，并求解Ly=b的过程。回带过程就是求解上三角方程组Ux=y。所以在实际的运算中，矩阵L和U可以直接计算出，而不需要任何中间步骤，从而在计算过程中将高斯消去法的步骤进行了进一步的简略，大大提高了运算速度，这就是三角分解法。通过以上的计算比较，2.题方程组具有严重的病态性。当系数矩阵有微小的变化时，wucha =-401.8918 159.5435 124.6330，所得的解与原方程组的解有很大的相对误差。1题方程组中当系数矩阵A和b有微小变化时，wucha =0 0 0 0，所得的解与方程组的解没有相对误差。所以1题方程组是良性的。用MATLAB内部函数inv通过求逆矩阵，然后通过x=i

23、nv（A）*b也可以求出方程组的解，但是没有列主元高斯消去法具有良好的稳定性。det函数求方程组系数矩阵的行列式时所得结果和高斯消去法和三角法所得结果相同，具有方便快捷的优点。题四可以看出，条件数越大，有效位数越少，当n=13时，出现有效位数为0的情况。附：高斯列主消去法源程序代码function x,det,index=Gauss(A,b)% ÇóÏßÐÔ·½³Ì×éµÄÁÐÖ÷ÔªGaussÏûÈ¥·¨£¬ÆäÖÐ% A - 方程组矩阵% b - 方程组右端% x -