数值分析实验讲解.docx

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

2、(0.02),end用迭代法求出x=2*x,的值，使得2xXx=1;while x+xx,x=x/2,pause(0.02),end用迭代法求出x=x/2,的最小值，使得2xX3.考虑下面二次代数方程的求解问题 公式是熟知的，与之等价地有，对于，应当如何选择算法。应该用计算，因为b与相近，两个相加减不宜做分母4.函数有幂级数展开利用幂级数计算的MATLAB程序为function s=powersin(x)s=0;t=x;n=1;while s+t=s; s=s+t； t=-x2/（(n+1)*(n+2)）*t； n=n+2；endt1=cputime;pause(10);t2=cputime;

3、t0=t2-t1(a)解释上述程序的终止准则。 当s+t=s，终止循环。(b)对于计算的进度是多少？分别计算多少项？ X=pi/2时，s =1.0000 x=11pi/2时，s=-1.0000 x=21pi/2时，s =0.9999Cputime分别是0.1563 0.0469 0.0156 5.考虑调和级数，它是微积分中的发散级数，在计算机上计算该级数的部分和，会得到怎么样的结果，为什么？function s=fun(n)s=0;t=1/n;for i=1:n s=s+1/i;end 当n=100时s =5.1874 当n=80时s =4.9655 当n=50时，s =4.4992当n=10

4、时，s =2.9290 6.指数函数的级数展开，如果对于，用上述的级数近似计算指数函数的值，这样的算法结果是否会好，为什么？function s=powerexp(x)s=1;n=1;t=1;while s+t=s; t=(xn)/factorial(n); s=s+t; n=n+1; end当x=-1时，s =0.3679 当x=-2时，s =0.1353 当x=-3时，s =0.04987.考虑数列，它的统计平均定义为 ，标准差数学上等价于作为标准差的两种算法，你将如何评价他们的得与失。clc,clearx=randn(1,10000);n=length(x);a=sum(x)/n;y1=

5、sqrt(sum(x-a).2)/(n-1); y2=sqrt(sum(x.2)-n*a2)/(n-1);y1,y2 后面的公式更好 改变m的值求出不同个数x标准差 ，没有好大差别实验二 插值法计算实习题1.已知函数在下列个点的值为0.20.40.60.81.00.980.920.810.640.38试用4次插值多项式及三次样条插值(自然边界条件)对数据进行插值。用图给出，及。 function f=lagfun(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:5 L(i)=1; for j=1:5 if j=i L(i)=

6、L(i)*(x-a(j)/(a(i)-a(j); end endendf=0;for i=1:5 f=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,newton(x0,y0,x),r);三次样条插值：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;y1=interp1(x,y,x1,spl

7、ine)plot(x1,y1)2.在区间-1,1上分别取用两组等距节点对龙格函数作多项式插值及三次样条插值，对每个值，分别画出插值函数及的图形。拉格朗日插值：function y=lagrl(x0,y0,x)n=length(x0);m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s;end再做12和20等距结点插值 for n=12:8:20 x=-1:0.01:1; y=1./(1+2

8、5.*x.2); z=0*x; x0=-1:2/(n-1):1; y0=1./(1+25.*x0.2); y1=lagrl(x0,y0,x);plot(x,z,r,x,y,k:,x,y1,r) gtext(Lagr.,num2str(n) hold on endtitle(Lagrange)legend(Lagr插值,f(x)图象)图象3.下列数据点的插值x01491625364964y012345678可以得到平方根函数的近似值，在区间0,64上作图。（1）用这九个点作8次多项式插值.（2）用三次样条（第一边界条件）程序求从得到结果看在0,64上，哪个插值更精确；在区间0,10上，两种插值哪

9、个更精确？（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,length(x)-1)poly2sym(a)答案-0.0000 0.0000 -0.0000 0.0002 -0.0050 0.0604 -0.3814 1.3257 -0.0000L(x)=-6345667567529957/19342813113834066795298816*x8+5071957851450983/75557863725914323419136*x7-3204839575550849/590295810358705651712*

10、x6+8226197088139413/36893488147419103232*x5-717795609662967/144115188075855872*x4+2177199843684719/36028797018963968*x3-3435436202510413/9007199254740992*x2+5970618836686703/4503599627370496*x-7696702421972085/2475880078570760549798248448画图x0=0,1,4,9,16,25,36,49,64;y0=sqrt(x);plot(x0,y0,r-)hold onx=

11、0,1,4,9,16,25,36,49,64;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=interp1(x,y,x1,spline)plot(x,sqrt(x),r-)hold onplot(x1,y1,b)hold on实验三 函数逼近与快速傅立叶变换 1. 对于给函数在区间-1,1上去，试求3次曲线拟合，试画出拟合曲线并打印方程

12、，与实验二，第二章的计算实习题2的结果比较。首先求出拟合曲线 x=-1:0.2:1;y=1./(25*x.2+1); p=polyfit(x,y,3)得出-1.6723e-016 -5.7518e-001 1.4919e-017 4.8412e-001Y=-0.5752x2+0.4814再输入x1=-1:0.01:1;y1=polyval(p,x1); plot(x,y,-*,x1,y1),grid2. 由实验给出数据表00.10.20.30.50.811.00.410.500.610.912.022.46试求3次，4次多项式的曲线拟合，再根据数据曲线图形，求一个另外函数的拟合曲线，用图示数据

13、曲线及相应的三种拟合曲线。先输入代码x=0.0 0.1 0.2 0.3 0.5 0.8 1.0; y=1.0 0.41 0.5 0.61 0.91 2.02 2.46;p3=polyfit(x,y,3)p4=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+000 9.2659e-0012.8853e+000 -1.2335e+001 1.6275e+001 -5.2987e+000 9.4272e-

14、001实验四 数值微分与数值积分 1. 用不同数值分析方法计算积分，（1） 取不同步长h。分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善？（2） 用龙贝格求积计算完成问题（1）。（3） 用自适应辛普森积分，使得精度达到。复合梯形公式建立m文件 clc,clear;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

15、.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(inline(sqrt(x).*log(x),0.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.99

16、38e-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.6548e009存在一个最小的h使精度不能在改善，且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);

17、while abs00.0001 k=k+1; h=h/2;p=0; for i=1:Mx=a+h*(2*i-1); p=p+subs(f,x); end T(k+1,1)=T(k,1)/2+h*p; M=2*M; for j=1:kT(k+1,j+1)=(4j)*T(k+1,j)-T(k,j)/(4j-1); endabs0=abs(T(k+1,j+1)-T(k,j); end q=T(k+1,k+1);n=k;在命令窗口输入Fx,n=shiyan42(sqrt(x)*log(x),10(-8),1)自适应辛普森积分Fx=quad(inline(sqrt(x).*log(x),0.00000

18、1,1)结果Fx=-0.4444实验五 数值微分 精度及步长的关系实验目的：数值计算中误差是不可避免的，希望同学能通过本实验初步认识数值分析中极端重要的两个概念：截断误差和舍入误差，并认真体会误差对计算结果的影响.问题提出：设一元函数，则在的导数定义为实验内容：计算在的导数值可以用如下的差分给出其近似： （1） （2）这也就给出了计算函数导数的两个简单算法.请给出几个计算高阶导数的近似算法，并完成如下工作：1. 对同样的h，比较（1）和（2）的计算结果.2. 针对计算高阶导数的算法，比较h取不同值时的计算结果.实验要求：选择有代表性的函数（最好选择多个），利用MATLAB提供的绘图工具画出该函

19、数在某个区间的导数曲线.再将数值计算的结果用MATLAB画出来，认真思考实验的结果.同学们还可以围绕这一问题设计其它的实验内容，特别认真的体会导数的结果。这是n阶导数的程序function C,L,dyk,k=ndaolag(X,Y,n)m=length(X);n1=m; L=ones(m,m);for k=1:m v=1; for i=1:m if k=i v=conv(v,poly(X(i)/(X(k)-X(i); endendL1(k,:)=v;l(k,:)=poly2sym(v);endC=Y*L1;L=Y*l;syms x dykfor k=1:nk;dyk=diff(L,x,k)

20、end在命令窗口输X=pi/6,pi/4,pi/3,5*pi/12,pi/2;Y=0.05,0.7071,0.8660,0.9659,1;C,L,dyk,k=ndaolag(X,Y,4)For i=1:4syms i x,dyi=diff(sin(x),x,i)end实验六 数值积分 误差传播与算法稳定性实验目的：体会算法稳定性在选择算法中的地位.问题提出：考虑一个简单的用积分定义的序列求解实验内容：由递推关系，可以得到计算积分序列的两种算法.算法一初始值为 （1）递推公式： （2）算法二初始值为 （3）递推公式： （4）实验要求：分别用算法一，算法二，并分别在计算中采用5位、6位和7位有效数

21、字，请判断哪种算法能给出更精确的结果。算法一：a=input(input youxiao weishu a: );syms n inin=vpa(exp(-1),a)for n=2:10 in=vpa(1-n*in),a)end结果 a=5 .36788 .26424 .20728 .17088 .14560 .12640 .11520 .7840e-1.29440 -1.9440算法二：function in=noibb=input(input youxiao weishu a: );syms n inin=vpa(0,b)for n=10:-1:2 in=vpa(1-in)/n),b)en

22、d结果a=50. .10000 .10000 .11250 .12679 .14554 .17089 .20728 .26424 .36788 实验七 多项式求根 病态问题实验目的：希望同学们通过本次实验对问题的病态性有一个初步的体会.问题提出：考虑一个高次的代数多项式 （1）显然该多项式的全部根为1,2，20,共计20个，且每个根都是单重的（也称为简单的），现在考虑该多项式的一个扰动 （2）其中是一个非常小的数，这相当于考虑（1）中的的系数有一个小的扰动，我们比较（1），（2）根的情况实验内容：为了实现方便，我们先介绍两个MATLAB函数：“roots”和“poly”.对u=roots(a)

23、若其中为维的向量，则该函数的输出为一个维的向量.设则的各分量是多项式方程的全部根.而函数若其中是一个维向量，该函数的输出是一个维向量，它是以的各分量为根的多项式的系数，可见“roots”和“poly”是两个互逆的运算.上述简单的MATLAB程序便得到（2）的全部根，程序中“ess”即使（2）中的 实验要求：选择充分小的ess，重复上述实验，记录结果的变化，还可以选择（2）中扰动项改为或者其它形式，观察会有怎样的现象出现。functiont_charp1_1%数值实验1.1变态问题%输入：020之间的扰动项及小的扰动常数%输出：加扰动都得到的全部根clcresult=inputdlg(请输入扰动

24、项：在020之间的整数：,charpt1-1,1,19);Numb=str2num(char(result);if(Numb20)|(Numb0)errordlg(请输入正确的扰动项:020之间的整数!);return;endresult=inputdlg(请输入（01）间的扰动常数：,charpt1_1,1,0.00001);ess=str2num(char(result);ve=zeros(1,21);ve(21-Numb)=ess;root=roots(poly(1:20)+ve);disp(对扰动项,num2str(Numb),加扰动,num2str(ess),得到的全部根为：);disp(num2str(root);结果分析：请输入扰动项：在020之间的整数：19请输入（01）间的扰动常数：0.00001