数值分析实验报告1

1、 实验1.1（病态问题）实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。题。通过本实验可获得一个初步体会。数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非问题提出：考虑一个高次的代数多项式显然该多项式的全部根为 1,2,20 共计 20 个，且每个根都是单重的。现考虑该多项式的一个扰动其中 是一个非常小的数。这相当于是对（1.1）中x 的系数作一个小的扰动。191.11.21.1）的解对扰动的敏感性。实验内容：为了实现方便，我们先介绍两个 Matlabroots和“其中若变量 a存储 n+1 u 为一个 n a的元素依次为a ,a , ,a ，则输出 u 的各分

2、量是多项式方程12n1的全部根；而函数的输出 b 是一个 n+1维变量，它是以 n 维变量 v 的各分量为根的多项式的系数。可见“roots和“poly”是两个互逆的运算函数。上述简单的 Matlab 程序便得到（1.2）的全部根，程序中的“ess即是（1.2）中的 。（1）选择充分小的 ess，反复进行上述实验，记录结果的变化并分析它们。如果扰动项的系数 很小，我们自然感觉（1.1）和（1.2）的解应当相差很小。计算中你有什么出乎意料的发现表明有些解关于如此的扰动敏感性如何（2）将方程（1.2）中的扰动项改成 x 或其它形式，实验中又有怎样的现18象出现（3）（选作部分）请从理论上分析产生这

3、一问题的根源。注意我们可以将方程（1.2）写成展开的形式，同时将方程的解 x 看成是系数的函数，考察方程的某个解关于的扰动是根关于的变化更敏感 实验过程：程序：a=poly(1:20);rr=roots(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)

4、;ve(n)=ess;a=poly(1:20)+ve;y=poly2sym(a);r=solve(y);-6-ms=max(abs(r-rr)endend数值实验结果及分析：format long0.9230n234560000007800000000000000000000000000000090n000000000000000000000000000000002300000000000000000000000000000000000000400000000000000000056789讨论：其根的扰动敏感性就越明显。实验总结：利用 MATLAB来进行病态问题的实验，虽然其得出的结果是有误差

5、的，但是可项式的根会有一定的扰动的，所以对于这类病态问题可以借助于 MATLAB来进行问题的分析。 实验 2.1（多项式插值的振荡现象）问题提出用的节点越多，插值多项式的次数就越高。 我们自然关心插值多项式的次数增加时，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)及插值多项式函

6、数 L(x)在-1，1上的图象，比较分析实验结果。（2）选择其它的函数，例如定义在区间-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)为插值节点构造上述各函数的拉格朗日插值多项式，比较其结果。实验过程：程序：多项式插值的震荡现象（实验 2.1）for m=1:6subplot(2,3,m)largrang(6*m)if m=1把窗口分割成 2*3 大小的窗口对 largrang函数进行运行t

7、itle(longn=6)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保存为：chazhi.mfunction largrang(longn)mm=input(please input mm(运行第几个函数就输入 mm 为几):mm=)if mm=1%d表示定义域的边界值d=1;elseif mm=2|mm=3d=5;endx0=linspac

8、e(-d,d,longn);if mm=1%x 的节点y0=1./(1.+25.*x0.2);elseif mm=2y0=x0./(1.+x0.4);elseif mm=3y0=atan(x0);endx=sym(x);n=length(x0); s=0.0;for k=1:np=1.0;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 mm=2ezplot(x/(1+x4)elseif mm=3ezplot(atan(x)endhold onezp

9、lot(y,-d,d)hold off保存为：largrang.m数值实验结果及分析：对于第一个函数 f(x)=1/(1+25x )2对于第二个函数 h(x)=x/(1+x )4对于第三个函数g(x)=arctan(x)讨论：通过对三个函数得出的 largrang象，说明了对函数的支点不是越多越好，而是在函数的两端而言支点越多，而largrang插值多项式不是更加靠近被逼近的函数，反而更加远离函数，在函数两端的跳动性更加明显，argrang插值多项式对函数不收敛。利用 MATLAB来进行函数的 largrang 插值多项式问题的实验，虽然其得出的结果是有误差的，但是增加支点的个数进行多次实验，

10、可以找出函数的largrang插值多项式的一般规律，当支点增加时，largrang 插值多项式对函数两端不收敛，不是更加逼近，而是更加远离，跳动性更强。所以对于函数的largrang插值多项式问题可以借助于 MATLAB来进行问题的分析，得到比较准确的实验结规律。 实验 5.1 （主元的选取与算法的稳定性）Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值 Gauss消去法作为数值算法的稳定性呢 Gauss 元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。实验内容：考虑线性方程组 Gauss 消去过程。实验要求：6 1 7 8 6 115 （1 , x n=10

11、8 6 1Ab*T 15 14 8 6计算矩阵的条件数。让程序自动选取主元，结果如何（2）现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模大的元素作为主元，结果又如何分析实验的结果。（3 n=20程中的作用。（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。程序：建立 M 文件：function x=gauss(n,r)请输入矩阵 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-范数：p=)

12、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,:)=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:

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

14、元： gauss请输入矩阵 A的阶数:n=10n =10条件数对应的范数是 p-范数：p=2p =pp=2请输入是否为手动，手动输入 1，自动输入 0：r=1r =1ans=1111111111选取绝对值最小的元素为主元： gauss请输入矩阵 A的阶数:n=10n =10条件数对应的范数是 p-范数：p=2p = 2pp =03请输入是否为手动，手动输入 1，自动输入 0：r=1r =1ans =1.000000000000001.000000000000001.000000000000011.000000000000001.000000000000001.000000000000001.0

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

16、0000000011.000000000000001.000000000000001.000000000000061.000000000000001.00000000000000999891.000000000000231.000000000000901.000000000012731.000000000003521.00000000002910将 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=

17、1p =1pp =请输入是否为手动，手动输入 1，自动输入 0：r=1r =1ans =1.000000000000511.000000002688051.000000000843371.00000000031354 gauss请输入矩阵 A的阶数:n=7n =7条件数对应的范数是 p-范数：p=2p =2pp =请输入是否为手动，手动输入 1，自动输入 0：r=1r = 1ans =1.000000000043371.000000001211431.00000000152825结果产生了较大的误差，条件数越大产生的误差就越大。讨论：在gauss选取绝对值大的元素作为主元比绝对值小的元素作为主

18、元时对结果产生的误差 gauss对用 gauss消去法解线性方程组时，主元的选取与算法的稳定性有密切的联数量级加大，使大数吃掉小数，产生舍入误差。 5.2问题提出：理论上，线性代数方程组 的摄动满足Ax b实，因为计算 A 通常要比求解方程Axb还困难。1实验内容：Matlab 中提供有函数“condest”可以用来估计矩阵的条件数，它给出的是按 范数的条件数。首先构造非奇异矩阵 A 和右端，使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动 A和b，使得 A和b 充分小。实验要求：（1）假设方程 Ax=b 的解为 x，求解方程(A)xbb，以 1-范数，给出xx的计算结果。xx（2c

19、ondest”所需机器时间的差别.很容易给出 cond (A)2的数值。将它与函数“”所得到的结果进行比较。x（3）利用“condest”给出矩阵A1）中的结果给出的x理论估计，并将它与（1）给出的计算结果进行比较，分析所得结果。注意，如果给出了 和 A 的估计，马上就可以给出的估计。A1（4）估计着名的 Hilbert 矩阵的条件数。程序：n=input(please input n:n=)a=fix(100*rand(n)+1x=ones(n,1)输入矩阵的阶数随机生成一个矩阵 a假设知道方程组的解全为 1用矩阵 a和以知解得出矩阵 b随即生成扰动矩阵 data随即生成扰动矩阵 datbb

20、=a*xdata=rand(n)*0.00001datb=rand(n,1)*0.00001A=a+dataB=b+datbxx=geshow(A,B)解扰动后的解xxx0=norm(xx-x,1)/norm(x,1)得出的理论结果xx保存为：fanshu.mfunction 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)*AB(i,i:nb);endendx=zeros(n,1);x(n)=AB(

21、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保存为：geshow.mfunction cond2(A)B=A*A;自定义求二阶条件数V1,D1=eig(B);V2,D2=eig(B(-1);cond2A=sqrt(max(max(D1)*sqrt(max(max(D2)end保存为：cond2.mformat longfor n=10:10:100n=n%n为矩阵的阶A=fix(100*randn(n); 随机生成矩阵 AcondestA=condest(A) 用 condest 求条件数

22、cond2(A)condA2=cond(A,2)pause用自定义的求条件数用 cond求条件数运行一次暂停end保存为：shiyan52.mn=input(please input n:n=)a=fix(100*rand(n)+1;x=ones(n,1);输入矩阵的阶数随机生成一个矩阵 a假设知道方程组的解全为 1用矩阵 a 和以知解得出矩阵 b随即生成扰动矩阵 data随即生成扰动矩阵 datbb=a*x;data=rand(n)*0.00001;datb=rand(n,1)*0.00001;A=a+data;B=b+datb;xx=geshow(A,B);利用第一小问的 geshow.m

23、求出解阵xxx0=norm(xx-x,1)/norm(x,1)得出的理论结果xxx00=cond(A)/(1-norm(inv(A)*norm(xx-x)*(norm(xx-x)/(norm(A)+norm(datb)/nxxorm(B)得出的估计值xxdatx=abs(x0-x00)保存为：sy5_2.m求两者之间的误差format longfor n=4:11n=n%n为矩阵的阶数生成 Hilbert 矩阵Hi=hilb(n);cond1Hi=cond(Hi,1)求 Hilbert 矩阵得三种条件数cond2Hi=cond(Hi,2)condinfHi=cond(Hi,inf)pausee

24、nd数值实验结果及分析： fanshuplease input n:n=6n =6a =1432142340142593405210 10031685881988485029719921323770896016322417227x =b =111111251 410 221 157 218 187data =1.0e-005 *datb =1.0e-005 *A =1.0e+002 *B =1.0e+002 *xx =x0 =xxx的计算结果为：xx（2）NcondestAe+002e+002e+002e+002e+002e+004e+003condA2e+002e+003e+002e+003

25、e+002e+002e+003e+003e+002e+002e+002e+002 sy5_2please input n:n=8n =8xx给出对的估计是：xxxx的理论结果是：xx结果相差：(4)cond1Hie+004e+005e+007e+008e+010e+012e+013e+015cond2Hie+004e+005e+007e+008e+010e+011e+013e+014condinfHie+004e+005e+007e+008e+010e+012e+013e+015n456789讨论：矩阵性质的一个重要的依据，条件数越大，矩阵“病态”性越严重，在解线性代求解。hilbert 矩阵

26、是一个很病态的矩阵，他的条过程中要尽量避免 hilbert矩阵 实验 7.1（迭代法、初始值与收敛性）及初始值与迭代收敛性的关系。问题提出：迭代法是求解非线性方程的基本思想方法，与线性方程的情况一样，其构造方法可以有多种多样，但关键是怎样才能使迭代收敛且有较快的收敛速度。实验内容：考虑一个简单的代数方程针对上述方程，可以构造多种迭代法，如在实轴上取初始值 x ，请分别用迭代（7.1（7.3）作实验，记录各算法的迭代过程。0实验要求：（1）取定某个初始值，分别计算（7.1）（7.3）迭代结果，它们的收敛性如何用 Matlab（2）对三个迭代法中的某个，取不同的初始值进行迭代，结果如何试分析迭代法

27、对不同的初值是否有差异（3）线性方程组迭代法的收敛性是不依赖初始值选取的。比较线性与非线性问题迭代的差异，有何结论和问题。程序：clearclcs=input(s=);clfif s=1a=-1.5;b=2.5;y00=0;x00=input(x00=);elseif s=2a=0.1;b=6.5;y00=0;x00=input(x00=);elseif s=3a=0;b=2; y00=0; x00=input(x00=);endx=linspace(a,b,80);y0=x;y=xy1=zxy7f(x,s);y=f(x)clear y;y=y0;y1;if s=1%plot(x,y,linewidth,1)legend(y=x,y=f1)title(x(n+1)=x(n)2-1)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)title(x(n+1)=sqrtx(n)+1)endhold onplot(a b,0,0,k-,0 0,a b,k-)axis(a,b,a,b)z=;for i=1:15次xt(1)=x00;yt