数值分析第一次作业

《数值分析第一次作业》由会员分享，可在线阅读，更多相关《数值分析第一次作业（18页珍藏版）》请在装配图网上搜索。

1、问题1：20.给定数据如下表：xj0.250.300.390.450.53yj0.50000.54770.62450.67080.7280试求三次样条插值S(x),并满足条件(l)S(0.25)=1.0000,S(0.53)=0.6868；2)S(0.25)=S(0.53)=0。分析：本问题是已知五个点,由这五个点求一三次样条插值函数。边界条件有两种已知一阶倒数,（2）是已知自然边界条件。对于第一种边界（已知边界的一阶倒数值）,可写出下面的矩阵方程。其中卩j=21000M，d，00Md11M=dM2d233Md44九00002九00H2九002H22九3300H24九0=1hj-1-h+hj-

2、1j入i=hdj=6fxj-1,xj,xj+1,Hn=1,j-1j6对于第一种边界条件d0=-0h6(fx0，X1-f0)，dn=hn-1解：由matlab计算得：xyh九d0.250.5000-5.52000.300.54770.05000.35711-4.31430.390.62450.09000.60000.6429-3.26670.450.67080.06000.42860.4000-2.42860.530.72800.08001.0000.5714-2.11500由此得矩阵形式的线性方程组为：21000一M，0-5.5200，0.357120.642900M1-4.31430.600

3、020.40000M=-3.266700.428620.5714M2-2.4286000123M-2.11504解得M0=-2.0286;M1=1.4627;M2=-1.0333;M3=-0.8058;M4=-0.6546S(x)=-6.76209(0.30-x)3-4.8758(x-0.25)3+10.0169(0.30-x)+10.9662(x-0.25)-2.708779(0.39-x)3-1.9136(x-0.30)3+6.1075(0.39-x)+6.9544(x-0.30)-2.87040(0.45-x)3-2.2384(x-0.39)3+10.4187(0.45-x)+11.18

4、8(x-0.39),-1.67881(0.53-x)3-1.3637(x-0.45)+8.3956(0.53-x)+9.1087(x-0.45)，xe0.25,0.30，xe0.30,0.39xe0.39，0.45xe0.45,0.53Matlab程序代码如下：functiontgsanci(n,s,t)%n代表元素数，s,t代表端点的一阶导。x=0.250.300.390.450.53;y=0.50.54770.62450.67080.7280;n=5,s=1.0,t=0.6868forj=1:1:n-1h(j)=x(j+1)-x(j);endforj=2:1:n-1r(j)=h(j)/(h

5、(j)+h(j-1);endforj=1:1:n-1u(j)=1-r(j);endforj=1:1:n-1f(j)=(y(j+1)-y(j)/h(j);endforj=2:1:n-1d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j);endd(1)=6*(f(1)-s)/h(1)d(n)=6*(t-f(n-1)/h(n-1)a=zeros(n,n);forj=1:1:na(j,j)=2;endr(1)=1;u(n)=1;forj=1:1:n-1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a)m=b*dp=zeros(n-l,4);%p矩阵为S(x

6、)函数的系数矩阵forj=1:1:n-1p(j,l)=m(j)/(6*h(j);p(j,2)=m(j+l)/(6*h(j);P(j,3)=(y(j)-m(j)*(h(jF2/6)/h(j);p(j,4)=(y(j+1)-m(j+1)*(h(j)人2/6)/h(j);end对于第二中边界，已知边界二阶倒数，可写出下面的矩阵：h其中卩.=U-Jh+hj-1jj-1jdj=6fxj-1,xj,xj+1，卩n=九0=0d0叫=0Md00Md11M=dM2d233Md442九0000卩2九00110卩2九0002卩22九33000卩24解：由matlab计算得：xyh九dn0.250.500000.30

7、0.54770.050.35710-4.31430.390.62450.090.60000.6429-3.26670.450.67080.060.42860.4000-2.42680.530.72800.0800.5714000.642920.42860由此得矩阵形式的线性方程组为：20.3571000020.600000M，0,0M1-4.3143M=-3.2667M2-2.42863M40_00000.4000020.571402解得M0=0;M=-1.8795;M2=-0.8636;M3=-1.0292;M4=0S(x)=0(0.30-x)3-6.2652(x-0.25)3+10(0.3

8、0-x)+10.9697(x-0.25),xe0.25,0.30-3.4806(0.39-x)3-1.5993(x-0.30)+6.1137(0.39-x)+6.9518(x-0.30),xe0.30,0.39-2.3990(0.45-x)3-2.8590(x-0.39)+10.417(0.45-x)+11.1903(x-0.39),xe0.39,0.45-2.1442(0.53-x)3-0(x-0.45)+8.3987(0.53-x)+9.1000(x-0.45),xe0.45,0.53matlab程序代码如下：functiontgsanci(n,s,t)%n代表元素数，x=0.250.30

9、0.390.450.53;y=0.50.54770.62450.67080.7280;n=5forj=1:1:n-1h(j)=x(j+1)-x(j);endforj=2:1:n-1r(j)=h(j)/(h(j)+h(j-1);endforj=1:1:n-1u(j)=1-r(j);endforj=1:1:n-1f(j)=(y(j+1)-y(j)/h(j);endforj=2:1:n-1d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j);endd(1)=0d(n)=0a=zeros(n,n);forj=1:1:na(j,j)=2;endr(1)=0;u(n)=0;forj=1:1:n

10、-1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a)k=am=b*dp=zeros(n-l,4);%p矩阵为S(x)函数的系数矩阵forj=1:1:n-1p(j,l)=m(j)/(6*h(j);p(j,2)=m(j+l)/(6*h(j);p(j,3)=(y(j)-m(j)*(h(j)人2/6)/h(j);p(j,4)=(y(j+1)-m(j+1)*(h(j)人2/6)/h(j);endp由下面的一段程序，画出自然边界与第一边界的图形：程序代码如下：x=0.250.300.390.450.53;y=0.50.54770.62450.67080.7280;s=cs

11、ape(x,y,variational)fnplt(s,r)holdonxlabel(x)ylabel(y)gtext(三次样条自然边界)plot(x,y,o,x,y,:g)holdons.coefsr=csape(x,y,complete,1.00.6868)fnplt(r,b)gtext(三次样条第一边界)holdonr.coefs图中的三条线几乎重合。x图20-2在整个给出的区间的图形问题2：1已知函数在下列各点的值为x.0.20.40.6.0.81.0f(x.)0.980.920.810.640.38试用4次牛顿插值多项式P4(x)及三次样条函数S(x)(自然边界条件)对数据进行插值。

12、用图给出(xi，yi)，令=0.2+0.08九i=0,1,11,10,P4(x)及S(x)。分析：(1)要用4次牛顿插值多项式处理数据，Pn=f(x0)+fx0,x1(x-x0)+fx0,X,x2(x-x0)(x-x1)+fx0,x1，xn(x-x0)(x-xn-1)首先我们要知道牛顿插值多项式的系数，即均差表中得部分均差。解：由matlab计算得：x.f(xj一阶差商一阶差商三阶差商四阶差商0.200.9800.400.920-0.300000.600.810-0.55000-0.625000.800.640-0.85000-0.75000-0.208331.000.380-1.30000-

13、1.12500-0.62500-0.52083所以有四次插值牛顿多项式为：P4(x)=0.98-0.3(x-0.2)-0.62500(x-0.2)(x-0.4)-0.20833(x-0.2)(x-0.4)(x-0.6)-0.52083(x-0.2)(x-0.4)(x-0.6)(x-0.8)计算牛顿插值中多项式系数的程序如下：functionvarargout=newtonliu(varargin)clear,clcx=0.20.40.60.81.0;fx=0.980.920.810.640.38;newtonchzh(x,fx);functionnewtonchzh(x,fx)%由此函数可得差

14、分表n=length(x);*j/f/vI:tfprintf(*分*n)；FF=ones(n,n);FF(:,1)=fx;fori=2:nforj=i:nFF(j,i)=(FF(j,i-1)-FF(j-1,i-1)/(x(j)-x(j-i+1);endendfori=1:nfprintf(%4.2f,x(i);forj=1:ifprintf(%10.5f,FF(i,j);endfprintf(n);end(2)用三次样条插值函数由上题分析知，要求各点的M值:有matlab编程计算求出个三次样条函数：20000，M，00，0.50020.50000M1-3.75000.50020.5000M=-

15、4.500000.50020.500M2-6.7500000023M04解得：M0=O;M1=-1.6071;M2=-1.0714;M3=-3.1071;M4=00(0.4-x)3-1.3393(x-0.2)+4.900(0.4-x)+4.6536(x-0.2)，x0.2,0.4-1.339(0.6-x)3-0.8929(x-0.4)+4.653(0.6-x)+4.0857(x-0.4)，x0.4,0.6-0.892(0.8-x)3-2.5893(x-0.6)+4.0857(0.8-x)+3.3036(x-0.6)，x0.6,0.8-2.589(1.0-x)-0(x-0.8)+3.3036(1

16、.0-x)+1.9(x-0.8)，x0.8,1.0三次样条插值函数计算的程序如下：functiontgsanci(n,s,t)%n代表元素数，s,t代表端点的一阶导。x=0.20.40.60.81.0;y=0.980.920.810.640.38;n=5forj=1:1:n-1h(j)=x(j+1)-x(j);endforj=2:1:n-1r(j)=h(j)/(h(j)+h(j-1);endforj=1:1:n-1u(j)=1-r(j);endforj=1:1:n-1f(j)=(y(j+1)-y(j)/h(j);endforj=2:1:n-1d(j)=6*(f(j)-f(j-1)/(h(j-1

17、)+h(j);endd(1)=0d(n)=0a=zeros(n,n);forj=1:1:na(j,j)=2;endr(1)=0;u(n)=0;forj=1:1:n-1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a)m=b*dp=zeros(n-l,4);%p矩阵为S(x)函数的系数矩阵forj=1:1:n-1p(j,l)=m(j)/(6*h(j);p(j,2)=m(j+l)/(6*h(j);p(j,3)=(y(j)-m(j)*(h(j)人2/6)/h(j);p(j,4)=(y(j+1)-m(j+1)*(h(j)人2/6)/h(j);endp下面是画牛顿插值以及

18、三次样条插值图形的程序：x=0.20.40.60.81.0;y=0.980.920.810.640.38;plot(x,y)holdonfori=1:1:5y(i)=0.98-0.3*(x(i)-0.2)-0.62500*(x(i)-0.2)*(x(i)-0.4)-0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0.52083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endk=011011x0=0.2+0.08*kfori=1:1:4y0(i)=0.98-0.3*(x(i)-0.2)-0.62500*(x(i)-0

19、.2)*(x(i)-0.4)-0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0.52083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endplot(x0,y0,o,x0,y0)holdony1=spline(x,y,x0)plot(x0,y1,o)holdons=csape(x,y,variational)fnplt(s,r)holdongtext(三次样条自然边界)gtext(原图像)gtext(4次牛顿插值)运行上述程序可知：给出的(xi，yi),Xi=0.2+0.08i,i=0,1,11,10点，S(x)及

20、P4(x)图形如下所示：问题3：3下列数据点的插值x01491625364964y012345678可以得到平方根函数的近似，在区间0,64上作图。(1) 用这9各点作8次多项式插值L8(x).(2) 用三次样条(自然边界条件)程序求S(x)。从结果看在0,64上，那个插值更精确；在区间0,1上，两种哪个更精确？分析：L8(x)可由公式Ln(x)=y1(x)得出。8nkkjk=0三次样条可以利用自然边界条件。写成矩阵：h其中卩=口-Jh+hj-1jh.j，1=h+hj-1jdj=6fxj-1,xj,xj+1，卩n=九0=0d0叫=0解：l(x)_(x_1)(x_4)(x_9)(x_16)(x_

21、25)(x_36)(x_49)(x_64)0(01)(04)(09)(016)(025)(036)(049)(064)(x_0)(x4)(x9)(x16)(x25)(x36)(x49)(x64)(10)(14)(19)(116)(125)(136)(149)(164)(x_0)(x-1)(x-9)(x-16)(x-25)(x-36)(x-49)(x-64)(4-0)(4-1)(4-9)(4-16)(4-25)(4-36)(4-49)(4-64)(x_0)(x1)(x4)(x16)(x25)(x36)(x49)(x64)(90)(91)(94)(916)(925)(936)(949)(964)(

22、x_0)(x1)(x4)(x9)(x25)(x36)(x49)(x64)(160)(161)(164)(169)(1625)(1636)(1649)(1664)(x-0)(x_1)(x-4)(x-9)(x-16)(x-36)(x-49)(x-64)(250)(251)(254)(259)(2516)(2536)(2549)(2564)(x-0)(x1)(x4)(x9)(x16)(x25)(x49)(x64)(360)(361)(364)(369)(3616)(3625)(3649)(3664)(x_0)(x_1)(x-4)(x-9)(x-16)(x-25)(x-36)(x-64)(490)(4

23、91)(494)(499)(4916)(4925)(4936)(4964)(x_0)(x_1)(x-4)(x-9)(x-16)(x-25)(x-36)(x-49)(64-0)(64-1)(64-4)(64-9)(64-16)(64-25)(64-36)(64-49)l1(x)=l2(x)=l3(x)=l4(x)=l5(x)=l6(x)=l7(x)=l8(x)=L8(x)=l1(x)+2l2(x)+3l3(x)+4l4(x)+5l5(x)+6l6(x)+7l7(x)+8l8(x)求三次样条插值函数由matlab计算：可得矩阵形式的线性方程组为：_2九0000_M0d012九100M1d102九0

24、Md002322九3M23d2300042M4d42.000000000000，M一0，0.25002.00000.75000000000M1M1.000.37502.00000.6250000000.1000.41672.00000.58330000M230.02860000.43752.00000.5625000M4=0.011900000.45002.00000.550000M50.0061000000.45832.00000.54170M60.00350000000.46342.00000.5357M70.0022000000002.0000_M80_解得:M0=0;M1=-0.520

25、9;M2=0.0558;M3=-0.0261;M4=0.0006;M5=-0.0029;M6=-0.0008;M7=-0.0009;01234567M8=00(1x)3+0.0868(x0)3+0(1X)+1.0868(x-0),xg0,1-0.0289(4x)3-0.0031(x-1)3+0.5938(4x)+0.6388(x1),xg1,40.0019(9x)30.0009(x-4)+0.3535(9x)+0.6271(x4),xg4,9S(x)=-0.0006(16x)3+0(x9)+0.4590(16x)-0.5708(x-9),xg9,16丿0s0(25x)3-0.0001(x16)

26、3-0.4436(25x)+0.5600(x-16),xg16,250(36x)-0(x25)3+0.4599(36x)+0.5470(x-25),xg25,360(48x)-0(x36)3+0.4633(48x)+0.5404(x-36),xg36,480(64x)3+0(x48)3+0.4689(64x)+0.5333(x-48),xg48,64拉格朗日插值函数与三次样条插值函数如图中所示,绿色实线条为三次样条插值曲线,蓝色虚线条为x=y2的曲线，另外一条红色线条为拉格朗日插值曲线。图3-1为01的曲线,图3-2为064区间上的曲线。由图3-1可以看出,红色的线条与蓝色虚线条几乎重合,所以

27、可知拉格朗日插值函数的曲线更接近开平方根的函数曲线,在01朗格朗日插值更精确。而在区间064上从图3-2中可以看出蓝色虚线条和绿色线条是几乎重合的,而红色线条在3070之间有很大的振荡,所在在区间064三次样条插值更精确写。图3-1图3-2Matlab程序代码如下：求三次样条插值函数的程序：functiontgsanci(n,s,t)%n代表元素数，s,t代表端点的一阶导。y=012345678;x=01491625364964;n=9forj=1:1:n-1h(j)=x(j+1)-x(j);endforj=2:1:n-1r(j)=h(j)/(h(j)+h(j-1);endforj=1:1:n

28、-1u(j)=1-r(j);endforj=1:1:n-1f(j)=(y(j+1)-y(j)/h(j);endforj=2:1:n-1d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j);endd(1)=0d(n)=0a=zeros(n,n);forj=1:1:na(j,j)=2;endr(1)=0;u(n)=0;forj=1:1:n-1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a)m=b*dt=ap=zeros(n-l,4);%p矩阵为S(x)函数的系数矩阵forj=1:1:n-1p(j,l)=m(j)/(6*h(j);p(j,2)=m(j+l

29、)/(6*h(j);p(j,3)=(y(j)-m(j)*(h(jF2/6)/h(j);p(j,4)=(y(j+1)-m(j+1)*(h(j)人2/6)/h(j);endp%画图形比较那个插值更精确的函数：x0=01491625364964;y0=012345678;x=0:64;y=lagr1(x0,y0,x);plot(x0,y0,o)holdonplot(x,y,r);holdon;pp=csape(x0,y0,variational)fnplt(pp,g);holdon;plot(x0,y0,:b);holdon%axis(0201);%看01区间的图形时加上这条指令gtext(三次样条

30、插值)gtext(原图像)gtext(拉格朗日插值)%下面是求拉格朗日插值的函数functiony=lagr1(x0,y0,x)n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj=kp=p*(z-x0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;（,）G,）厂a、（f,）00D0nD00D,（,）（,）丿a丿Sf,）0,b0)。在matlab中分别用上述拟合函数拟合曲线，本题应采用倒指数拟合，即y(t)=a*exp(bt-1)拟合曲线。对y(t)=a*exp(bt-i)两边取对数

31、有Iny=lna+b/t;如令Y=Iny,A=lna,则得Y=A+b/t;为确定A,b,先将(ti,y/转化为(t/Yj).根据最小二乘法，取(t)，1,(t)，1/1。01用lsqcurvefit函数拟合曲线，程序代码如下：创建M文件：functionF=mf(a,t)F=a(l)*exp(a(2).*t.A(-l);在命令窗口中敲入如下代码：t=5:5:55y=l.272.l62.863.443.874.l54.374.5l4.584.624.64;a0=0.5,0.5;a=lsqcurvefit(pf,a0,t,y)回车后可得：a(l)=5.503l,a(2)=-8.7872下面的计算问

32、题用matlab编程实现：x=5l0l52025303540455055;y0=l.272.l62.863.443.874.l54.374.5l4.584.624.64;y=log(y0)y(l)=0r=zeros(2,2);fori=l:l:l2;r(1,1)=r(1,1)+1;endfori=2:1:12r(1,2)=r(1,2)+1/x(i);endfori=2:1:12r(2,1)=r(2,1)+1/x(i);endfori=2:1:12r(2,2)=r(2,2)+l/x(i)人2;enda=zeros(2,l);fori=2:l:l2a(l,l)=a(l,l)+y(i);a(2,l)

33、=a(2,l)+y(i)*(l/x(i);endk=rt=ad=inv(r);m=d*a由matlab软件计算得：a=3.9864,b=-4.8922。所以y(t)=3.9864e-4.8922/t。有下面的语句给出拟合后的图形：x=05l0l52025303540455055;y0=0l.272.l62.863.443.874.l54.374.5l4.584.624.64;fplot(5.5031*exp(-8.7872)*(x).人(-1),0,55);holdonplot(x,y0,o,x,y0,r)holdongtext(倒指数拟合曲线)gtext(原图像)原图像倒指数拟合曲线图18-2