数值分析实验2014讲解

1、数值分析实验(2014,9,1610,28)信计1201班，人数34人数学系机房数值分析计算实习报告册专业学号姓名20142015年第一学期实验一数值计算的工具Matlab1. 解释下MATLABS序的输出结果程序：t=0.1n=1:10e=n/10-n*te 的结果：00 -5.5511e-01700-1.1102e-016 -1.1102e-0160 002. 下面MATLABS序的的功能是什么？程序：x=1;while 1+x1,x=x/2,pause(0.02),e nd用迭代法求出x=x/2,的最小值x=1;while x+xx,x=2*x,pause(0.02),e nd用迭代法求

2、出x=2*x,的值，使得2xXx=1;while x+xx,x=x/2,pause(0.02),e nd用迭代法求出x=x/2,的最小值，使得2xX3. 考虑下面二次代数方程的求解问题2ax bx c = 02公式b 如是熟知的，与之等价地有X： 2C，对于2a_b 、b2_4aca =1,b =100000000, c =1，应当如何选择算法。应该用X = b_上2 _4ac计算，因为b与、b2】4ac相近，两个相加减不宜2a2做分母3574. 函数sin(x)有幕级数展开sin-.3! 5!7!利用幕级数计算sinx的MATLAB程序为fun cti on s=powers in(x)s=

3、0;t=x;n=1;while s+t=s;s=s+t ；t=-xA2/( (n+1)*(n+2) )*t ；n=n+2 ；endt仁cputime;pause(10);t2=cputime;t0=t2-t1(a) 解释上述程序的终止准则。当s+t=s，终止循环。(b) 对于x二二/2,11二/2,21二/2计算的进度是多少？分别计算多少项？X=pi/2 时，s =1.0000 x=11pi/2时，s=-1.0000 x=21pi/2 时，s =0.9999Cputime 分别是 0.15630.04690.01565. 考虑调和级数a丄，它是微积分中的发散级数，在计算机上计算该级数的部n =

4、1 n分和，会得到怎么样的结果，为什么？fun ctio ns=fu n(n)s=0;t=1/ n;for i=1: ns=s+1/i;end当 n=80 时 s =4.9655 当 n=50 时，s =4.4992当 n=100 时 s =5.1874当 n=10 时，s =2.9290236.指数函数的级数展开e1 x .,如果对于x ：0，用上述的级数近 2!3!似计算指数函数的值，这样的算法结果是否会好，为什么？function s=powerexp(x)s=1;n=1;t=1;while s+t=s;t=(xA n) /factorial( n);s=s+t;n=n+1;end当 x

5、=-1 时，s =0.3679 当 x=-2 时，s =0.1353 当 x=-3 时，s =0.04987.考虑数列xi,i =1,2,.n，它的统计平均定义为一、 -1 n1 n 2-2 I2标准差o- = I 无(人-x)2数学上等价于=(E xi -nx )作为标准差屮 -1 i 二 n -1 7的两种算法，你将如何评价他们的得与失。clc,clearx=randn (1,10000);n=len gth(x);a=sum(x)/n;y1= sqrt(sum(x-a).A2)/( n-1);y2=sqrt(sum(x.A2)-n*aA2)/( n-1);y1,y2后面的公式更好改变m的

6、值求出不同个数x标准差，没有好大差别实验二插值法计算实习题1.已知函数在下列个点的值为xi0.20.40.60.81.0f (Xi)0.980.920.810.640.38试用4次插值多项式p4(x)及三次样条插值S(x)(自然边界条件)对数据进行插值。用图给出(Xi，y)Xi =0.2 0.081,0,1,11,10，p4(x)及S(x)。fun cti on f=lagf un(x)a=0.2,0.4,0.6,0.8,1.0;b=0.98,0.92,0.81,0.64,0.38;for i=1:5L(i)=1;for j=1:5if j=iL(i)=L(i)*(x-a(j)/(a(i)-a

7、(j);endendendf=0;for i=1:5f=f+L(i)*b(i);end执行文件x0=0.2,0.4,0.6,0.8,1.0;y0=0.98,0.92,0.81,0.64,0.38;plot(x0,y0,o)hold ongrid onfplot(lagfun,0,1);hold onx=0:0.1:1;plot(x, newt on (x0,y0,x),r);00.1020.30.40.6060.7 DB 091三次样条插值：x=0.2,0.4,0.6,0.8,1.0;y=0.98,0.92,0.81,0.64,0.38;x1=0.2:0.08:1;y仁in terp1(x,y

8、,x1,spli ne)plot(x1,y1)70.9 “0.8 -0.7 f0.6 -0.5 -0.20.30.40.50.60.70.80.92.在区间-1,1上分别取n =12,20用两组等距节点对龙格函数f(x)21 +25x2 多项式插值及三次样条插值，对每个 n值，分别画出插值函数及f (x)的图形。 拉格朗日插值：fun cti on y=lagrl(xO,yO,x)n=len gth(x0);m=le ngth(x);for i=1:mz=x(i);s=0.0;for k=1: np=1.0;for j=1: nif j=kp=p*(z-x0(j)/(x0(k)-x0(j);e

9、ndend s=p*y0(k)+s;endy(i)=s;end再做 12和 20等距结点插值for n=12:8:20 x=-1:0.01:1;y=1./(1+25.*x.A2);z=0*x;x0=-1:2/(n-1):1; y0=1./(1+25.*x0.A2);y1=lagrl(x0,y0,x);plot(x,z,r,x,y,k:,x,y1,r) gtext(Lagr.,num2str(n)hold onendtitle(Lagrange) legend(Lagr 插值 ,f(x)图象)图象Lagrange3.下列数据点的插值x01491625364964x=0,1,4,9,16,25,3

10、6,49,64;11y012345678可以得到平方根函数的近似值，在区间0,64上作图。(1) 用这九个点作8次多项式插值L8(x).(2) 用三次样条(第一边界条件)程序求从得到结果看在0,64上，哪个插值更精确；在区间0,10上，两种插值哪个更 精确？(1) 多项式插值x=0,1,4,9,16,25,36,49,64;y=0,1,2,3,4,5,6,7,8;a=polyfit(x,y,le ngth(x)-1)poly2sym(a)答案-0.00000.0000-0.00000.0002-0.00500.0604-0.3814 1.3257 -0.0000L(x)=-6345667567

11、529957/19342813113834066795298816*xA8+5071957851450983/75557863725914323419136*xA7-3204839575550849/590295810358705651712* xA6+8226197088139413/36893488147419103232*xA5-717795609662967/1441151 88075855872*xA4+2177199843684719/36028797018963968*xA3-34354362025104 13/9007199254740992*xA2+5970618836686

12、703/4503599627370496*x-769670242 1972085/2475880078570760549798248448画图 x0=0,1,4,9,16,25,36,49,64;y0=sqrt(x);plot(x0,y0,r-)hold on y=0,1,2,3,4,5,6,7,8;a=polyfit(x,y,8); xx=0:0.1:64; yy=polyval(a,xx); plot(xx,yy,b) hold on(3)三次样条x=0,1,4,9,16,25,36,49,64;y=0,1,2,3,4,5,6,7,8;x1=0:0.1:64;y1= in terp1(x,

13、y,x1,spli ne)plot(x,sqrt(x),r-) hold onplot(x1,y1,b) hold on实验三函数逼近与快速傅立叶变换11.对于给函数f(x)2在区间-1,1上去x-l 0.2i(i =1,2,川,10)，试1 +25x求3次曲线拟合，试画出拟合曲线并打印方程，与实验二，第二章的计算实习题2的结果比较。首先求出拟合曲线x=-1:0.2:1;y=1./(25*x.A2+1);p=polyfit(x,y,3)Y=-0.5752xA2+0.4814再输入x1=-1:0.01:1;y 仁polyval(p,x1);2. 由实验给出数据表x00.10.20.30.50.8

14、1y1.00.410.500.610.912.022.46试求3次，4次多项式的曲线拟合，再根据数据曲线图形，求一个另外函数的拟合曲线，用图示数据曲线及相应的三种拟合曲线。先输入代码x=0.0 0.1 0.2 0.3 0.5 0.8 1.0;p3=polyfit(x,y,3)y=1.0 0.41 0.5 0.61 0.91 2.02 2.46;15p4=polyfit(x,y,4)z=0.0:0.1:1.0;y3=polyval(p3,z);y4=polyval(p4,z);plot(x,y,b,z,y3,r,z,y4,g)结果-6.6221e+000 1.2815e+001 -4.6591e

15、+000 9.2659e-0012.8853e+000 -1.2335e+001 1.6275e+001 -5.2987e+000 9.4272e-0012.521.510.500.10.20.30.40.50.60.70.80.91实验四数值微分与数值积分1.用不同数值分析方法计算积分 ixlnxdx二4 ,09(1)取不同步长h。分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否 存在一个最小的h，使得精度不能再被改善？(2) 用龙贝格求积计算完成问题(1)(3) 用自适应辛普森积分，使得精度达到10J。复合梯形公式建立m文件clc,cl

16、ear;x=0.0000001:0.1:1;y=sqrt(x). *log(x);Fx=trapz(x,y)error=Fx+4/9结果 Fx=-0.4123error=0.0321h=0.01，Fx=-0.4431，error=-0.0014；h=0.001，Fx=-0.4444，error=-5.4875e-005； h=0.0001，Fx=-0.4444，error=-2.0338e-006； h=0.00001，Fx=-0.4444，error=-6.4496e-008.没有一个最小的h使精度不能再改变复合辛普森求积Fx=quad( in li ne(sqrt(x)*log(x),0.

17、000001,1,0.0001)error=Fx+4/9结果 Fx=-0.4440error=4.3273e-004h=0.00001， Fx=-0.4444 error=-6.3537e-005；h=0.000001， Fx=-0.4444，error=-2.9938e-006；h=0.0000001， Fx=-0.4444， error=-3.0657e-007 ；h=0.0000000000001， Fx=-0.4444error=-9.6548e-009;h=0.00000000000001， Fx=-0.4444error=-9.6548e- 009存在一个最小的h使精度不能在改善，

18、且 h在0.0000000000001与 0.00000000000001 之间。龙贝格求积function q,n=shiyan42(f,a,b)M=1;abs0=10;k=0;T=zeros(1,1);h=b-a;T(1,1)=(h/2)*(subs(f,a)+subs(f,b);while abs00.0001k=k+1;h=h/2;p=0;for i=1:Mx=a+h*(2*i-1);p=p+subs(f,x);endT(k+1,1)=T(k,1)/2+h*p;M=2*M;for j=1:kT(k+1,j+1)=(4Aj)*T(k+1,j)-T(k,j)/(4Aj-1);endabs0

19、=abs(T(k+1,j+1)-T(k,j);endq=T(k+1,k+1);n=k;在命令窗口输入Fx,n=shiyan42(sqrt(x)*log(x),10A(-8),1) 自适应辛普森积分Fx=quad(inline(sqrt(x).*log(x),0.000001,1)结果 Fx=-0.4444实验五 数值微分 精度及步长的关系实验目的：数值计算中误差是不可避免的， 希望同学能通过本实验初步认识数值分析中 极端重要的两个概念： 截断误差和舍入误差， 并认真体会误差对计算结果的影响问题提出：设一元函数f ：R R，贝U f在xo的导数定义为f (X。)二 hm0f (Xo h) 一 f

20、 (Xo)h23实验内容:计算在x0的导数值可以用如下的差分给出其近似:f (Xo)f (Xo h) f (Xo)h(1)f (Xo)f (xo h) - f (xo - h)2h(2)这也就给出了计算函数导数的两个简单算法.请给出几个计算高阶导数的近似算法，并完成如下工作:1.对同样的h,比较(1)和(2)的计算结果.2.针对计算高阶导数的算法，比较 h取不同值时的计算结果实验要求:选择有代表性的函数f(x)(最好选择多个)，利用MATLAB提供的绘图 工具画出该函数在某个区间的导数曲线f(x).再将数值计算的结果用 MATLAB画出来，认真思考实验的结果.同学们还可以围绕这一问题设计其它的

21、 实验内容，特别认真的体会导数的结果。这是 n 阶导数的程序 functionC,L,dyk,k=ndaolag(X,Y,n)m=le ngth(X); n1=m; L=on es(m,m);for k=1:mv=1;for i=1:mif k=iv=co nv(v,poly(X(i)/(X(k)-X(i);endendL1(k, :)=v;l (k, :)=poly2sym(v);endC=Y*L1;L=Y*I;syms x dykfor k=1: nk；dyk=diff(L,x,k)end在命令窗输X=pi/6,pi/4,pi/3,5*pi/12,pi/2;Y=0.05,0.7071,0.

22、8660,0.9659,1;C,L ,dyk,k=ndaolag(X,Y,4)For i=1:4syms i x,dyi=diff(si n(x),x,i)end实验六数值积分误差传播与算法稳定性实验目的：体会算法稳定性在选择算法中的地位.问题提出：考虑一个简单的用积分定义的序列1En 二 0 xnexdx, n 二 0,1,2,|1|,求解 En, n =0,1,2,3,HI实验内容：由递推关系En=1 - nEn4, n =1,2,3川I，可以得到计算 巳=xnexdx, n =0,1,2,111,积分序列 洽的两种算法.算法一初始值为E=1/e（ 1）递推公式：En=1-En4, n -

23、1,2,IH算法初始值为En =0( 3)递推公式：EnJ -, n 二N,N_1,N_2,|(,3,2,1(4)n实验要求：分别用算法一，算法二，并分别在计算中采用5位、6位和7位有效数字,请判断哪种算法能给出更精确的结果。算法一：a=in put( in put youxiao weishu a:);syms n inin=vpa(exp(-1),a)for n=2:10in=vpa(1- n*in ),a)end结果a=5.36788 .26424 .20728 .17088 .14560 .12640 .11520 .7840e-1.29440 -1.9440算法二：fun cti o

24、n in=noibb=in put( in put youxiao weishu a:);syms n inin=vpa(0,b)for n=10:-1:2in=vpa(1-i n)/n),b)end结果a=50. .10000 .10000 .11250 .12679 .14554 .17089 .20728.26424 .36788实验七多项式求根病态问题实验目的:希望同学们通过本次实验对问题的病态性有一个初步的体会.问题提出：考虑一个高次的代数多项式20P(x) =(x-1)(x-2)|(x-20)H (x-k)( 1)kA显然该多项式的全部根为1,2,20,共计20个，且每个根都是单重的(也 称为简单的)，现在考虑该多项式的一个扰动P(x)亠 7.x19 =0其中；是一个非常小的数，这相当于考虑(1)中的x19的系数有一个小的扰动, 我们比较(1), (2)根的情况实验内容：为了实现方便，我们先介绍两个 MATLAB函数：“roots”