数值分析实验报告1

1、实验一误差分析实验（病态问题）实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言， 所谓坏问题就是问题本身对扰动敏感者， 反之属于好问题。通过本实验可获得一个初步体会。数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。 病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价 （如耗用更多的机器时间、 占用更多的存储空间等）。问题提出：考虑一个高次的代数多项式20p(x) ( x 1)( x 2) (x 20)(x k)(1.1)k1显然该多项式的全部根为 1,2, ,20 共计 20 个，且每个根都是

2、单重的。现考虑该多项式的一个扰动p( x)x190(1.2)其中是一个非常小的数。这相当于是对（）中x19 的系数作一个小的扰动。我们希望比较（）和（）根的差别，从而分析方程（）的解对扰动的敏感性。实验内容： 为了实现方便，我们先介绍两个Matlab 函数：“roots”和“ poly”。uroots(a)其中若变量 a 存储 n+1 维的向量，则该函数的输出 u 为一个 n 维的向量。设 a 的元素依次为 a1, a2 , , an 1 ，则输出 u 的各分量是多项式方程a1 xna2 xn 1an xan 10的全部根；而函数bpoly(v)的输出 b 是一个 n+1 维变量，它是以n 维

3、变量 v 的各分量为根的多项式的系数。可见“ roots”和“ poly”是两个互逆的运算函数。ess0.000000001;vezeros(1,21);ve(2)ess;roots ( poly (1: 20)ve)上述简单的 Matlab 程序便得到（）的全部根，程序中的“ess”即是（）中的。实验要求：（1）选择充分小的 ess，反复进行上述实验，记录结果的变化并分析它们。如果扰动项的系数 很小，我们自然感觉（）和（）的解应当相差很小。计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？（2）将方程（）中的扰动项改成x18 或其它形式，实验中又有怎样的现象出现？（3）（选作

4、部分）请从理论上分析产生这一问题的根源。注意我们可以将方程（）写成展开的形式，p( x,)x20x190(1.3)同时将方程的解 x 看成是系数 的函数，考察方程的某个解关于 的扰动是否敏感，与研究它关于 的导数的大小有何关系？为什么？你发现了什么现象，哪些根关于 的变化更敏感？思考题一：（上述实验的改进）在上述实验中我们会发现用roots 函数求解多项式方程的精度不高，为此你可以考虑用符号函数 solve 来提高解的精确度， 这需要用到将多项式转换为符号多项式的函数 poly2sym, 函数的具体使用方法可参考 Matlab 的帮助。实验过程：程序：a=poly(1:20);rr=roots

5、(a);for n=2:21nfor m=1:9ess=10(-6-m);ve=zeros(1,21);ve(n)=ess;r=roots(a+ve);-6-ms=max(abs(r-rr)endend利用符号函数：（思考题一）a=poly(1:20);y=poly2sym(a);rr=solve(y)for n=2:21nfor m=1:8ess=10(-6-m);ve=zeros(1,21);ve(n)=ess;a=poly(1:20)+ve;y=poly2sym(a);r=solve(y);-6-ms=max(abs(r-rr)endend数值实验结果及分析：format long-6-

6、m-7-8-9-10n234050060000700008000090000100000110000120000130000140000150000160000170000180000190000200000210000-6-m-11-12-13-14n203000400005000060000700008000090000100000110000120000130000140000150000160000170000180000190000200000210000讨论：利用这种方法进行这类实验， 可以很精确的扰动敏感性的一般规律。 即当对扰动项的系数越来越小时， 对其多项式扰动的结果也就越来

7、越小， 即扰动敏感性与扰动项的系数成正比， 扰动项的系数越大， 对其根的扰动敏感性就越明显， 当扰动的系数一定时， 扰动敏感性与扰动的项的幂数成正比， 扰动的项的幂数越高， 对其根的扰动敏感性就越明显。实验总结：利用 MATLAB 来进行病态问题的实验，虽然其得出的结果是有误差的，但是可以很容易的得出对一个多次的代数多项式的其中某一项进行很小的扰动，对其多项式的根会有一定的扰动的，所以对于这类病态问题可以借助于 MATLAB 来进行问题的分析。学号： 06450210姓名：万轩实验二插值法实验（多项式插值的振荡现象）问题提出 ：考虑一个固定的区间上用插值逼近一个函数。显然拉格朗日插值中使用的节

8、点越多，插值多项式的次数就越高。 我们自然关心插值多项式的次数增加时，L(x) 是否也更加靠近被逼近的函数。 龙格给出了一个极着名例子。 设区间-1 ， 1 上函数f(x)=1(1+25x2)实验内容： 考虑区间 -1 ，1 的一个等距划分，分点为：x(i)=-1+2i/n,i=0,1,2,n泽拉格朗日插值多项式为：L(x)=l(i)(x)/(1+25x(j)2 ) i=0,1,n其中 l(i)(x), i=0,1,n,n 是 n 次拉格朗日插值基函数。实验要求： 选择不断增大的分点数目 n=2,3 ，画出 f(x) 及插值多项式函数 L(x) 在 -1 ，1 上的图象，比较分析实验结果。（2

9、）选择其它的函数，例如定义在区间-5，5 上的函数h(x)=x/(1+x4) , g(x)=arctanx重复上述的实验看其结果如何。（3）区间 a,b 上切比雪夫点的定义为：xk=（b+a）/2+(b-a)/2)cos(2k-1)/(2(n+1),k=1,2,n+1以 x1,x2x(n+1) 为插值节点构造上述各函数的拉格朗日插值多项式，比较其结果。实验过程：程序：多项式插值的震荡现象（实验）for m=1:6subplot(2,3,m)largrang(6*m)%把窗口分割成 2*3 大小的窗口%对 largrang 函数进行运行if m=1title('longn=6')

10、elseif m=2title('longn=12')elseif m=3title('longn=18')elseif m=4title('longn=24')elseif m=5title('longn=30')elseif m=6title('longn=36')end%对每个窗口分别写上标题为插值点的个数end保存为：function largrang(longn)mm=input('please input mm(运行第几个函数就输入mm 为几 ):mm=')if mm=1%d 表示定义域

11、的边界值d=1;elseif mm=2|mm=3d=5;endx0=linspace(-d,d,longn); %x 的节点 if mm=1y0=1./(1.+25.*x0.2);elseif mm=2y0=x0./(1.+x0.4);elseif mm=3y0=atan(x0);endx=sym('x');n=length(x0); s=;for k=1:np=;for j=1:nif j=kp=p*(x-x0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;endy=s;if mm=1ezplot('1/(1+25*x2)')elseif

12、 mm=2ezplot('x/(1+x4)')elseif mm=3ezplot('atan(x)')endhold onezplot(y,-d,d)hold off保存为：数值实验结果及分析：对于第一个函数f(x)=1/(1+25x 2)0.60.40.20-11.510.50-0.5-10.50-0.5-5420-2-4-5longn=6longn=12longn=18110.50.500-0.501-101-101xxxlongn=24longn=30longn=3621.51100.5-10-2-0.501-101-101xxx对于第二个函数h(x)=x

13、/(1+x 4)longn=6longn=12longn=181100-1-105-505-505xxxlongn=24longn=3020longn=36200-2-2005-505-505xxx对于第三个函数g(x)=arctan(x)longn=6longn=12longn=181221000-1-1-2-2-505-505-505xxxlongn=24longn=30longn=36222000-2-2-2-505-505-505xxx讨论：通过对三个函数得出的 largrang 插值多项式并在数学软件中的运行， 得出函数图象，说明了对函数的支点不是越多越好，而是在函数的两端而言支点越

14、多，而 largrang 插值多项式不是更加靠近被逼近的函数，反而更加远离函数，在函数两端的跳动性更加明显， argrang 插值多项式对函数不收敛。实验总结：利用 MATLAB 来进行函数的 largrang 插值多项式问题的实验，虽然其得出的结果是有误差的，但是增加支点的个数进行多次实验，可以找出函数的 largrang 插值多项式的一般规律，当支点增加时， largrang 插值多项式对函数两端不收敛，不是更加逼近，而是更加远离，跳动性更强。所以对于函数的 largrang 插值多项式问题可以借助于 MATLAB 来进行问题的分析，得到比较准确的实验结规律。学号： 06450210姓名：

15、万轩实验五解线性方程组的直接方法实验 （主元的选取与算法的稳定性）问题提出： Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的， 如何才能确保 Gauss消去法作为数值算法的稳定性呢？ Gauss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。实验内容：考虑线性方程组Axb,AR n n , bR n编制一个能自动选取主元， 又能手动选取主元的求解线性方程组的 Gauss消去过程。实验要求：61786115（ 1）取矩阵 A,b，则方程有解 x*(1,1, ,1)T 。取 n=

16、10861158614计算矩阵的条件数。让程序自动选取主元，结果如何？（ 2）现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元， 观察并记录计算结果。 若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。（ 3）取矩阵阶数 n=20 或者更大，重复上述实验过程， 观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异， 说明主元素的选取在消去过程中的作用。（ 4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。实验过程：程序：建立 M 文件：function x=gauss(n,r

17、)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)p=input('条件数对应的范数是p-范数： p=')pp=cond(A,p)pausem,n=size(A);nb=n+1;Ab=A br=input(' 请输入是否为手动，手动输入1，自动输入 0： r=')for i=1:n-1if r=0pivot,p=max(abs(Ab(i:n,i);ip=p+i-1;if ip=iAb(i ip,:

18、)=Ab(ip i,:);disp(Ab); pauseendendif r=1i=iip=input(' 输入 i 列所选元素所处的行数：ip=');Ab(i ip,:)=Ab(ip i,:);disp(Ab); pauseendpivot=Ab(i,i);for k=i+1:nAb(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb);enddisp(Ab); pauseendx=zeros(n,1);x(n)=Ab(n,nb)/Ab(n,n);for i=n-1:-1:1x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1

19、:n)/Ab(i,i);end数值实验结果及分析：取矩阵 A 的阶数 :n=10，自动选取主元：>> format long>> gauss请输入矩阵 A 的阶数 :n=10n =10条件数对应的范数是p-范数： p=1p =1pp =+003请输入是否为手动，手动输入1，自动输入 0：r=0r =0取矩阵 A 的阶数 :n=10，手动选取主元：选取绝对值最大的元素为主元：>> gauss请输入矩阵 A 的阶数 :n=10n =10条件数对应的范数是p- 范数： p=2p =2pp=+003请输入是否为手动，手动输入1，自动输入 0：r=1r =1ans=1

20、111111111选取绝对值最小的元素为主元：>> gauss请输入矩阵 A 的阶数 :n=10n =10条件数对应的范数是p- 范数： p=2p =2pp =+003请输入是否为手动，手动输入1，自动输入 0：r=1r =1ans =取矩阵 A 的阶数 :n=20，手动选取主元： 选取绝对值最大的元素为主元：>> gauss请输入矩阵 A 的阶数 :n=20条件数对应的范数是p- 范数： p=1p =1pp =+006ans =11111111111111111111 选取绝对值最小的元素为主元：>> gauss请输入矩阵 A 的阶数 :n=20.n =2

21、0条件数对应的范数是p- 范数： p=2p =2pp =+006请输入是否为手动，手动输入1，自动输入 0：r=1r =1ans =将 M文件中的第三行：A=diag(6*ones(1,n)+diag(ones(1,n-1),1)+diag(8*ones(1,n-1),-1)改为：A=hilb(n) >> gauss请输入矩阵 A 的阶数 :n=7n =7条件数对应的范数是p- 范数： p=1p =1pp =+008请输入是否为手动，手动输入1，自动输入 0：r=1r =1ans = >> gauss请输入矩阵 A 的阶数 :n=7n =7条件数对应的范数是p- 范数：

22、 p=2p =2pp =+008请输入是否为手动，手动输入1，自动输入 0：r=1r =1ans =该问题在主元选取与算出结果有着很大的关系， 取绝对值大的元素作为主元比取绝对值小的元素作为主元时产生的结果比较准确， 即选取绝对值小的主元时结果产生了较大的误差，条件数越大产生的误差就越大。讨论：在 gauss消去法解线性方程组时， 主元的选择与算法的稳定性有密切的联系，选取绝对值大的元素作为主元比绝对值小的元素作为主元时对结果产生的误差较小。条件数越大对用 gauss消去法解线性方程组时， 对结果产生的误差就越大。实验总结：对用 gauss消去法解线性方程组时，主元的选取与算法的稳定性有密切的

23、联系，选取适当的主元有利于得出稳定的算法， 在算法的过程中， 选取绝对值较大的主元比选取绝对值较小的主元更有利于算法的稳定， 选取绝对值最大的元素作为主元时，得出的结果相对较准确较稳定。 条件数越小， 对用这种方法得出的结果更准确。在算除法的过程中要尽量避免使用较小的数做为除数， 以免发生结果数量级加大，使大数吃掉小数，产生舍入误差。学号： 06450210姓名：万轩实验（线性代数方程组的性态与条件数的估计）问题提出：理论上，线性代数方程组Axb 的摄动满足xcond ( A)Abx1 A 1AAb矩阵的条件数确实是对矩阵病态性的刻画， 但在实际应用中直接计算它显然不现实，因为计算 A 1 通

24、常要比求解方程 Ax b 还困难。实验内容： Matlab 中提供有函数“ condest”可以用来估计矩阵的条件数，它给出的是按 1-范数的条件数。首先构造非奇异矩阵 A 和右端，使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动A和 b ，使得A 和b 充分小。实验要求：（ 1）假设方程Ax=b 的解为 x，求解方程 ( AA)x?bb ，以 1-范数，给出x x? x 的计算结果。x x（ 2）选择一系列维数递增的矩阵（可以是随机生成的） ，比较函数“ condest”所需机器时间的差别 .考虑若干逆是已知的矩阵， 借助函数“eig”很容易给出 cond2(A) 的数值。将它与函

25、数“ cond(A,2)”所得到的结果进行比较。（ 3）利用“condest”给出矩阵A 条件数的估计， 针对（1）中的结果给出x的x理论估计，并将它与（ 1）给出的计算结果进行比较，分析所得结果。注意，如果给出了 cond(A)和 A 的估计，马上就可以给出 A 1 的估计。（ 4）估计着名的 Hilbert 矩阵的条件数。H(hi , j ) n n , hi, j1, i , j 1,2, , nij 1实验过程：程序：n=input('please input n:n=')%输入矩阵的阶数a=fix(100*rand(n)+1x=ones(n,1)b=a*x%随机生成一

26、个矩阵a%假设知道方程组的解全为%用矩阵 a 和以知解得出矩阵1bdata=rand(n)*%随即生成扰动矩阵datadatb=rand(n,1)*%随即生成扰动矩阵datbA=a+dataB=b+datbxx=geshow(A,B)%解扰动后的解x0=norm(xx-x,1)/norm(x,1)%得出?x 的理论结果xxxx保存为：function x=geshow(A,B)%用高斯消去法解方程组m,n=size(A);nb=n+1;AB=A B;for i=1:n-1pivot=AB(i,i);for k=i+1:nAB(k,i:nb)=AB(k,i:nb)-(AB(k,i)/pivot)

27、*AB(i,i:nb);endendx=zeros(n,1);x(n)=AB(n,nb)/AB(n,n);for i=n-1:-1:1x(i)=(AB(i,nb)-AB(i,i+1:n)*x(i+1:n)/AB(i,i);end保存为：function cond2(A)%自定义求二阶条件数B=A'*A;V1,D1=eig(B);V2,D2=eig(B(-1);cond2A=sqrt(max(max(D1)*sqrt(max(max(D2)end保存为：format longfor n=10:10:100n=n%n 为矩阵的阶A=fix(100*randn(n);condestA=con

28、dest(A)cond2(A)condA2=cond(A,2)pause%随机生成矩阵 A%用 condest求条件数%用自定义的求条件数%用 cond 求条件数%运行一次暂停end保存为：n=input('please input n:n=')a=fix(100*rand(n)+1;%输入矩阵的阶数%随机生成一个矩阵ax=ones(n,1);%假设知道方程组的解全为 1b=a*x;%用矩阵 a 和以知解得出矩阵 bdata=rand(n)*;%随即生成扰动矩阵 datadatb=rand(n,1)*;%随即生成扰动矩阵 datbA=a+data;B=b+datb;xx=ges

29、how(A,B);%利用第一小问的求出解阵x0=norm(xx-x,1)/norm(x,1)%得出x?x x 的理论结果xxx00=cond(A)/(1-norm(inv(A)*norm(xx-x)*(norm(xx-x)/(norm(A)+norm(datb)/norm(B)%得出x?x x 的估计值xxdatx=abs(x0-x00)%求两者之间的误差保存为：format longfor n=4:11n=n%n 为矩阵的阶数Hi=hilb(n);%生成 Hilbert 矩阵cond1Hi=cond(Hi,1)%求 Hilbert 矩阵得三种条件数cond2Hi=cond(Hi,2)cond

30、infHi=cond(Hi,inf)pauseend数值实验结果及分析： >> fanshuplease input n:n=6n =6a =142516881989329385489260144088501316235219292324010100737241437227701x =111111b =251410221157218187data =*datb =*A =+002 *B =+002 *xx =x0 =x x? x 的计算结果为：x x（ 2）NcondestAcond2AcondA210e+00220e+00230e+002e+002e+00240e+00250e+

31、00260e+004e+003e+00370e+003e+002e+00280e+00290e+003e+002e+002100e+003e+002e+002>> sy5_2please input n:n=8n = 8x0 =x00 =datx =给出对x?的估计是：x xx xx x? x 的理论结果是：x x结果相差：(4)ncond1Hicond2HicondinfHi4e+004e+004e+0045e+005e+005e+0056e+007e+007e+0077e+008e+008e+0088e+010e+010e+0109e+012e+011e+01210e+013e

32、+013e+01311e+015e+014e+015讨论：线性代数方程组的性态与条件数有着很重要的关系， 既矩阵的条件数是刻画矩阵性质的一个重要的依据，条件数越大，矩阵“病态”性越严重，在解线性代数方程组的过程中较容易产生比较大的误差， 则在实际问题的操作过程中， 我们必须要减少对条件数来求解， 把条件数较大的矩阵化成条件数较小的矩阵来进行求解。实验总结：在本次实验中， 使我们知道了矩阵条件数对线性代数方程组求解的影响， 条件数越大，对最后解的影响的越大， hilbert 矩阵是一个很 ”病态 ”的矩阵，他的条件数随着阶数的增加而增大， 每增加一阶， 条件数就增大一个数量级， 在求解的过程中要

33、尽量避免 hilbert 矩阵学号： 06450210姓名：万轩实验七非线性方程求根实验（迭代法、初始值与收敛性）实验目的：初步认识非线性问题的迭代法与线性问题迭代法的差别， 探讨迭代法及初始值与迭代收敛性的关系。问题提出：迭代法是求解非线性方程的基本思想方法，与线性方程的情况一样，其构造方法可以有多种多样，但关键是怎样才能使迭代收敛且有较快的收敛速度。实验内容：考虑一个简单的代数方程x 2x10针对上述方程，可以构造多种迭代法，如xn 1xn21(7.1)xn 11(7.2)1xnxn 1xn 1(7.3)在实轴上取初始值x0，请分别用迭代（） -（）作实验，记录各算法的迭代过程。实验要求：

34、（ 1）取定某个初始值，分别计算（） -（）迭代结果，它们的收敛性如何？重复选取不同的初始值，反复实验。请自选设计一种比较形象的记录方式（如利用Matlab 的图形功能），分析三种迭代法的收敛性与初值选取的关系。（ 2）对三个迭代法中的某个，取不同的初始值进行迭代，结果如何？试分析迭代法对不同的初值是否有差异？（ 3）线性方程组迭代法的收敛性是不依赖初始值选取的。比较线性与非线性问题迭代的差异，有何结论和问题。实验过程：程序：clearclcs=input('请输入要运行的方程，运行第几个输入几s=');clfif s=1% 决定坐标轴的范围和初始值a=;b=; y00=0;

35、x00=input('请输入第一个函数的初值：x00=');elseif s=2a=;b=; y00=0; x00=input('请输入第二个函数的初值：x00=');elseif s=3a=0;b=2; y00=0; x00=input('endx=linspace(a,b,80);y0=x;y1=zxy7f(x,s);clear y;y=y0;y1;if s=1plot(x,y,'linewidth',1)legend('y=x','y=f1')title('x(n+1)=x(n)2-1'

36、;)elseif s=2plot(x,y,'linewidth',2)legend('y=x','y=f2')title('x(n+1)=1+1/x(n)')elseif s=3plot(x,y,'linewidth',3)legend('y=x','y=f3')请输入第三个函数的初值：x00=');%计算直线y=x% 计算迭代函数 y=f(x)%画图% 输出标题title('x(n+1)=sqrtx(n)+1')endhold onplot(a b,0,0,

37、'k-',0 0,a b,'k-')axis(a,b,a,b)%画坐标轴z=;for i=1:15% 画蛛网图，迭代过程为n=15 次xt(1)=x00;yt(1)=y00;%决定始点坐标xt(2)=zxy7f(xt(1),s);%决定终点坐标yt(2)=zxy7f(xt(1),s);zxyplot7(xt,yt,%画蛛网图if i<=5pause%按任意键逐次观察前5 次迭代的蛛网图endx00=xt(2);y00=yt(2);% 将本次迭代的终点作为下次的始点z=z,xt(1);%保存迭代点end保存为：function y=zxy7f(x,s)if

38、s=1y=(x.*x-1);elseif s=2y=(1+1./x);elseif s=3y=sqrt(x+1);end保存为：functionout=zxyplot7(x,y,p)%画一次迭代的蛛网图，改变p 调节箭头的大小u(1)=0;v(1)=(y(2)-y(1);%画出始点 (x(1),y(1)终点 (x(2),y(2)的有向折线段u(2)=eps;v(2)=eps;h=quiver(x(1) x(1),y(1) y(2),u,v,p);set(h,'color','red');hold onu(1)=(x(2)-x(1);v(1)=0;u(2)=eps;v(2)=eps;h=quiver(x(1) x(2),y(2) y(2),u,v,p);set(h,'color','red');plot(x(1) x(1) x(2),y(1) y(2) y(2),'')保存为：数