数值分析实验报告.doc
《数值分析实验报告.doc》由会员分享,可在线阅读,更多相关《数值分析实验报告.doc(40页珍藏版)》请在三一文库上搜索。
1、实验五 解线性方程组的直接方法实验5.1 (主元的选取与算法的稳定性)问题提出:Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的,如何才能确保Gauss消去法作为数值算法的稳定性呢?Gauss消去法从理论算法到数值算法,其关键是主元的选择。主元的选择从数学理论上看起来平凡,它却是数值分析中十分典型的问题。实验内容:考虑线性方程组 编制一个能自动选取主元,又能手动选取主元的求解线性方程组的Gauss消去过程。实验要求:(1)取矩阵,则方程有解。取n=10计算矩阵的条件数。让程序自动选取主元,结果如何?(2)现选择程序中手动选取主元的功能。每步消
2、去过程总选取按模最小或按模尽可能小的元素作为主元,观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元,结果又如何?分析实验的结果。(3)取矩阵阶数n=20或者更大,重复上述实验过程,观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异,说明主元素的选取在消去过程中的作用。(4)选取其他你感兴趣的问题或者随机生成矩阵,计算其条件数。重复上述实验,观察记录并分析实验结果。思考题一:(Vadermonde矩阵)设 ,其中,(1)对n=2,5,8,计算A的条件数;随n增大,矩阵性态如何变化?(2)对n=5,解方程组Ax=b;设A的最后一个元素有扰动10-4,再求解Ax=b(3
3、计算(2)扰动相对误差与解的相对偏差,分析它们与条件数的关系。(4)你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗?相关MATLAB函数提示:zeros(m,n) 生成m行,n列的零矩阵ones(m,n) 生成m行,n列的元素全为1的矩阵eye(n) 生成n阶单位矩阵rand(m,n) 生成m行,n列(0,1)上均匀分布的随机矩阵diag(x) 返回由向量x的元素构成的对角矩阵tril(A) 提取矩阵A的下三角部分生成下三角矩阵triu(A) 提取矩阵A的上三角部分生成上三角矩阵rank(A) 返回矩阵A的秩det(A) 返回方阵A的行列式inv(A)
4、 返回可逆方阵A的逆矩阵V,D=eig(A) 返回方阵A的特征值和特征向量norm(A,p) 矩阵或向量的p范数cond(A,p) 矩阵的条件数L,U,P=lu(A) 选列主元LU分解R=chol(X) 平方根分解Hi=hilb(n) 生成n阶Hilbert矩阵5.1实验过程:5.1.1程序:function x=gauss(n,r)n=input(请输入矩阵A的阶数:n=)A=diag(6*ones(1,n)+diag(ones(1,n-1),1)+diag(8*ones(1,n-1),-1)b=A*ones(n,1)for i=1:4p=input(条件数对应的范数是p-范数:p=)pp=
5、cond(A,p)endpausem,n=size(A);nb=n+1;Ab=A br=input(请输入是否为手动,手动输入1,自动输入0:r=)for i=1:n-1 if r=0 pivot,p=max(abs(Ab(i:n,i); ip=p+i-1; if ip=i Ab(i ip,:)=Ab(ip i,:);disp(Ab); pause end end if r=1 i=i ip=input(输入i列所选元素所处的行数:ip=); Ab(i ip,:)=Ab(ip i,:);disp(Ab); pause end pivot=Ab(i,i); for k=i+1:n Ab(k,i:
6、nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb); end disp(Ab); pauseendx=zeros(n,1);x(n)=Ab(n,nb)/Ab(n,n);for i=n-1:-1:1 x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n)/Ab(i,i);end5.1.2实验结果如下: 1.按照实验要求一:取矩阵A的阶数:n=10且自动选取主元,程序结果运行如下:(2) 现选择程序中手动选取主元的功能,观察并记录计算结果。 选取绝对值最大的元素为主元:程序运行开始如第一问的截图也是求范数故这里不在给出。The answer is :
7、1 1 1 1 1 1 1 1 1 1 选取绝对值最小的元素为主元:The answer is: 1.0e+003*(INF 0.007 0.0057 -0.0410 -0.0303 0.3430 0.2577 -2.7290 -2.0463 2.7308)取矩阵A的阶数:n=20,手动选取主元: 选取绝对值最大的元素为主元:The answer is :1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 选取绝对值最小的元素为主元:The answer is: 1.0e+007*(-Inf 0.0000 0.0000 -0.0000 -0.0000 0.0000
8、0.0000 -0.0003 -0.0002 0.0022 0.0016 -0.0175 -0.0131 0.1398 0.1049 -1.1185 -0.8389 8.9478 6.7109 -8.9478)修改程序如下:function x=gaussong(n,r)n=input(请输入矩阵A的阶数:n=)A=hilb(n)b=A*ones(n,1)for i=1:4p=input(条件数对应的范数是p-范数:p=)pp=cond(A,p)endpausem,n=size(A);nb=n+1;Ab=A br=input(请输入是否为手动,手动输入1,自动输入0:r=)for i=1:n-
9、1 if r=0 pivot,p=max(abs(Ab(i:n,i); ip=p+i-1; if ip=i Ab(i ip,:)=Ab(ip i,:);disp(Ab); pause end end if r=1 i=i ip=input(输入i列所选元素所处的行数:ip=); Ab(i ip,:)=Ab(ip i,:);disp(Ab); pause end pivot=Ab(i,i); for k=i+1:n Ab(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb); end disp(Ab); pauseendx=zeros(n,1);x(n)=A
10、b(n,nb)/Ab(n,n);for i=n-1:-1:1 x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n)/Ab(i,i);end所求范数为:自动输入结果为:ans = 1.0000 1.00001 .0000 1.0000 1.0002 0.9996 1.0007 0.9993 1.0004 0.9999选取绝对值最大的元素为主元结果为:The answer is :NaN NaN NaN NaN NaN Inf -Inf -Inf 281.3945 -283.3708选取绝对值最小的元素为主元结果为:The answer is :NaN NaN NaN -Inf
11、5.8976 -1.9243 -2.0291 -4.9972 23.4548 -11.10125.1.3 对实验结果进行分析:5.1.3.1 对实验要求一的结果进行分析:对于Gauss消去法就是用行的初等变换将原线性方程组系数矩阵转化为简单形式,从而进行求解,缺点是迭代次数可能较多,效率不高,且在消去过程中不可以将主元素很小的做除数,否则将导致其他元素数量级的严重增长和舍入误差的扩散,使得计算解不可靠。5.1.3.2 对实验要求二的结果进行分析:通过每次选取最大或最小的主元可以发现取绝对值大的元素作为主元比取绝对值小的元素作为主元时产生的结果比较准确,即选取绝对值小的主元时结果产生了较大的误
12、差,条件数越大产生的误差就越大。所以应尽量避免很小的数作为除数。5.1.3.3 对实验要求三的结果进行分析:此要求是对要求一和要求二的一个延续,通过实验结果可以看出若采用很小的数作为主元迭代次数越多导致的结果越不可靠,甚至出现错误。5.1.3.4 对实验要求四的结果进行分析:对新矩阵进行实验发现依然符合上述规律,可以知道,在进行迭代时主元的选择与算法的稳定性有密切的联系选取绝对值大的元素作为主元比绝对值小的元素作为主元时对结果产生的误差较小。条件数越大对用gauss消去法解线性方程组时,对结果产生的误差就越大。5.1.4实验总结:1. 在用gauss消去法解线性方程组时,主元的选取与算法的稳定
13、性有密切的联系,选取适当的主元有利于得出稳定的算法,2. 在算法的过程中,选取绝对值较大的主元比选取绝对值较小的主元更有利于算法的稳定,选取绝对值最大的元素作为主元时,得出的结果相对较准确较稳定。3. 条件数越小,对用这种方法得出的结果更准确。4. 在算除法的过程中要尽量避免使用较小的数做为除数,以免发生结果数量级加大,使大数吃掉小数,产生舍入误差。实验5.2(线性代数方程组的性态与条件数的估计)问题提出:理论上,线性代数方程组的摄动满足 矩阵的条件数确实是对矩阵病态性的刻画,但在实际应用中直接计算它显然不现实,因为计算通常要比求解方程还困难。实验内容:Matlab中提供有函数“condest
14、可以用来估计矩阵的条件数,它给出的是按1-范数的条件数。首先构造非奇异矩阵A和右端,使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动,使得充分小。实验要求:(1)假设方程Ax=b的解为x,求解方程,以1-范数,给出的计算结果。(2)选择一系列维数递增的矩阵(可以是随机生成的),比较函数“condest”所需机器时间的差别.考虑若干逆是已知的矩阵,借助函数“eig”很容易给出cond2(A)的数值。将它与函数“cond(A,2)”所得到的结果进行比较。(3)利用“condest”给出矩阵A条件数的估计,针对(1)中的结果给出的理论估计,并将它与(1)给出的计算结果进行比较,分析所得结
15、果。注意,如果给出了cond(A)和的估计,马上就可以给出的估计。(4)估计著名的Hilbert矩阵的条件数。5.2 实验过程如下: 5.2.1.1 实验要求一的程序如下:function n=jisuan(n)a=fix(100*rand(n)+1 x=ones(n,1) b=a*x data=rand(n)*0.00001 datb=rand(n,1)*0.00001 A=a+dataB=b+datbx0=get(A,B) x1=norm(x0-x,1)/norm(x,1) function x=get(A,B)m,n=size(A);nb=n+1;AB=A B;for i=1:n-1 p
16、ivot=AB(i,i); for k=i+1:n AB(k,i:nb)=AB(k,i:nb)-(AB(k,i)/pivot)*AB(i,i:nb); endendx=zeros(n,1);x(n)=AB(n,nb)/AB(n,n);for i=n-1:-1:1 x(i)=(AB(i,nb)-AB(i,i+1:n)*x(i+1:n)/AB(i,i);End5.2.1.2 实验要求一程序运行结果如下:系数矩阵a为:7018142952536398247346580289074421962620992338538830595879897467437769加扰动后的系数矩阵A为:70.0000046
17、 18.0000042 14.0000099 29.0000032 52.0000044 53.0000013 63.0000057 98.0000030 2.0000079 47.0000096 34.0000093 65.0000021 80.0000079 28.0000087 90.0000044 7.0000073 44.0000068 21.0000061 96.0000006 26.0000002 20.0000050 99.0000041 23.0000021 38.0000063 53.0000060 88.0000077 30.0000021 59.0000074 58.0
18、000084 79.0000037 89.0000005 74.0000097 67.0000064 43.0000027 77.0000063 69.0000058 b值为:236309270302367419加扰动后的b值为:236.00000451309.00000044270.00000027302.00000313367.00000013419.00000384 data的值为:4.610951267E-064.153748604E-069.900825926E-063.200355775E-064.399243096E-061.337727485E-065.678287124E-0
19、63.04998677E-067.888616922E-069.600986004E-069.333801082E-062.071327296E-067.942106514E-068.743671716E-064.386585338E-067.266317666E-066.833323243E-066.071989445E-065.91825935E-071.50094987E-074.983113034E-064.119532082E-062.125598643E-066.298878488E-066.028690857E-067.6795039E-062.139633316E-067.44
20、5657831E-068.392382403E-063.704768261E-065.02688037E-079.708449393E-066.434922879E-062.679472507E-066.287846E-065.75147779E-06 datb的值为:4.514248268E-064.38953253E-072.7185123E-073.126850481E-061.28625747E-073.839672885E-06 xx的值为:0.999999684491.00000038780.999999143481.00000065171.00000221560.99999754
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 实验 报告
