数值分析第一次作业

1、XJ0.250.300.390450.53yj05000054770.624506708C7280问题1：20.给定数据如下表:试求三次样条插值S(x),并满足条件(l)S'(025)=10000,S'(053)=06868：(2)S,(025)=S>>(053)=0o分析：本问题是已知五个点，由这五个点求一三次样条插值函数。边界条件有两种，（1）是已知一阶倒数，（2）是已知自然边界条件。对于第一种边界（已知边界的一阶倒数值），可写出卜.面的矩阵方程。2OOO42从2OOO42“3O111111其中尸匚，一k+hjhj.+hj4n=l，%于第一种边界条件&=

2、9-(fIxo,Xi-6),*一(fn-f'xn.1,xn>hohn.l解:由matlab计算得:XyhAAd0.250.5000-5.52000.3005477005000.35711-4.314303906245009000.600006429-3.266704506708006000.428604000-2.428605307280008001.00005714-2.1150由此得矩阵形式的线性方程组为:21000Mo-5.520003571o0642900-4.3143006000040000-32667000428605714m3-2.428600012-2.1150解

3、得0286,Mi=l46”,M尸-10333,My=-08058,岫=06546S(x>-676209030-x)3-48758x-025)3-F100169030-x)+109662x-025),xe025,0,30-2.708779039-x)3-19136(x-0.30)3+6.1075C0.39-x)4-69544x-030),xe0.30.039-2.8704Q0,45-x)3-2.2384x-0.39)3+104187045-x)+l1188x-0.39),xe0,39,045L6788Q0.53-x)3-L3637cx0.45)3+8395d0.53-x)+9,1087x-

4、045),xe0.45.0.53Matlab程序代码如卜.：functiontgsaiici(n,s,t)%n代表元素数，s.t代表端点的阶存。x=0.25030039045053,y=0505477062450670807280,n=5,s=l0,t=06868forj=lLn-1h(jKG+l)-xG),endforj=2:ln-1。冲6)他。)+蚂-1),endforj=ll:n-lendfbrj=lln-1endforj=2:l:n-lendd(l)=6*(f(l)-s)/h(l)d(n)=6*(t-f(n-l)/h(n-l)a=zeros(n,n),forj=lInendKA1,u(

5、n)=l;fora(j+lj)=uG+l),aj+lE)，endb=inv(a)m=b*d'p=zeros(n-l,4),%p矩阵为S(x)函数的系数矩阵forp(j,l)=m®/(6*h(j),P(j3Hy0-m(j)*0i(jr2/6)/h0；pG，4HyG+l)-mG+D*(hGT2/6)/hG),end对于第二中边界，己知边界二阶倒数，可写出卜面的矩阵:2%000Modo%00Midi02当0M=*003o一当M3d30002.m4.A.其中4J=，4尸dj=6fXj.i,Xj不j+J4n=70=0&)=nW%+hj忆+hj解：由matlab计算得:Xyh4A

6、dn0.250.5000003005477005035710-43143039062450090600006429-32667045067080060428604000-24268053072800080057140由此得矩阵形式的线性方程组为:035710000O106000000064292042860000004000020571402Mo必m3m3m40-43143-3.2667-242860解得Mq=0,Mi=-18795,M尸-08636,M3=-10292,4=0S(x>0(0.30-x)3-6.2651x-0.25)3+10(0.30-X)4-1O.9697(X-0.25

7、),xG0.25,030一3.48000.39-x)3-1.5999x0.30>+6.11330.39x)+6.951Ix0.30),xe0.30,039-2.399tt0.45-x)3-2.859ax-0.39)3+10.417(0.45-x)+11.1901x-0.39),xG0.39,0.452.14420.53-x)3-0(x-0.45)3+83987(0.53-x)+9.1000(x-0.45),xe0.45,0.53matlab程序代码如下：functiontgsanci(n,s,t)%n代表元素数,x=0.250300.390.45053,y=050547706245067

8、0807280),n=5forj=l1n-1MA加+DM),endforj=2:ln-13曲颂冲。-1),endforj=l:ln-1uGAi-B);endforendforj=2:l:n-lendd(l>0d(n>0a=zeros(h,n),fbrj=llnendr(l)=0,u(n>0,forj=l1n-1endb=inv(a)k=am=b*d'p=zeros(n-l,4),%p矩阵为S(x)函数的系数矩阵forj=l:l:n-lpG，gnG)/(6*hG);endP由卜面的一段程序，画出自然边界与第一边界的图形：程序代码如下：x=0.250.3003904505

9、3,y=0505477062450670807280,ssapeCx.y.Variational')娜1心尸)holdonxlabelCx')ylabel(y)gtextC三次样条白然边界i)plotCx,y/o'.x.y；g')holdons.coefsr=csape(xytomplete1.006868)fhpltgb1)gtextC三次样条第一边界i)holdonrcoefs图中的三条线几乎重合。图20-1在一小段区间内的图形图20-2在整个给出的区间的图形问题2；1已知函数在卜列各点的值为X10.204060810f(Xx)0.980920810.640

10、38试用4次牛顿插值多项式P4(x)及三次样条函数S(x)(白然边界条件)对数据进行插值。用图给出(知力),Xf=02+008i,1=0,1,11,10,P4(x)及S(x).分析：(1)要用4次牛顿插值多项式处理数据，Pn=f(XQ)+fIXo,Xi(X-Xo)+fXo,XbX2(X-Xo)(X-XjH+4河,，一4伏-沏)(X-X»i)首先我们要知道牛顿插值多项式的系数，即均差表中得部分均差。解：由matlab计算得:片fg一阶差商二阶差商三阶差商四阶差商0.2009800.400.920-0.300000.600.810-0.55000-0.625000800640-08500

11、0-0.75000-0.208331.000380-1.3U000-1.12500-062500-0.52083所以有四次插值牛顿多项式为:P4(x)=098-03(x-02)-0.62500(x-02)(x-04)-020833(x-0,2)(x-04)(x-0.6)-052083(x-0.2)(x-04)(x-0.6)(x-0.8)计算牛顿插值中多项式系数的程序如下：fanctionvarargout=newtonhu(yarargin)cleai;clcx=020,4060810,fc=0980.920810640.38,newtonchzh(x,fe),fanctionnewtonch

12、zh(x,fx)%由此函数可得差分表n=length(X).OrintfC*"*差分表*"),FF=ones(n,n),FFC.l>fe',for1=2:nforj=i:nendfori=lnfprintf(%42f,x(i),forj=l:i(2)用二次样条插值函数由上题分析知，要求各点的M值;有matlab编程计算求出个三次样条函数：0.5000.5000.500020.5000.500020.500M。MiM,M；M4-3.7500-4.5000-6.7500解得：Mo=O,M-1,6071,M尸-10714;M卡-3,1071;M«=00(0

13、.4-x)3-1.339Xx-0.2)+4.90(X04-x)+4.6536(x-0.2),xg0.2,0.4-1.33950.6-x)3-0.8929(x-0.4)3+4.65360.6-x)+4.0857(x-0.4),xe0,4,0.6-0.89290.8-x)3-2.5893：x-0.6)3+4.085X0.8-x)+3.3036x-0.6),xg0.6,0.8-2.589J1.0-x)3-0(x-0.8)3+3.3036a.0-x)+1.9(x-0.8),xgO.8JL.O三次样条插值函数计算的程序如下：functiontgsanci(n,s,t)%n代表元素数，s,t代表端点的一阶

14、导。x=0204060.81.0,y=0,980920810.64038,n=5fbrj=l1n-1endforj=2:l:n-lr(j)=hG)/(gh(H);endforU(J>1-r(j),endforj=l1n-1fGRG+i)-yG)/hG),endforj=2:l:n-lendd(D=0d(n)=0a=zei-os(n,n),forj=l1nendr(l)=0,u(n>0,fbrj=l1n-1aG+igG+D,aGj+lK0,endb=inv(a)rn=b*d'p=zeros(n-l,4),%p矩阵为S(x)函数的系数矩阵fbrj=l:ln-1pGJFiG)/(

15、6*h。)，pG，3Xy(j)-mG)忡92/6)颂);P&4)=(yG+i)-mG+i)*SGA2/6)/hG),endP卜面是画牛顿插住以及三次样条插值图形的程序：x=0204060810,y=(098092081064038,plot(x,y)holdonfori=l15y(i>098-03*(x(i)-02)-062500*(x(i)-02)*(x(i)-04)920833*(x(1)-02尸(x(i)-0.4)*-06)-052083*(x(i)92)*(x(i)-04)*(x(i)96)*(X(i)98)endk=(oiionx0=0.2+008*kfor1=1,14

16、y0(i)=098-0.3*(x(i)-0.2)-062500*(x(i)-0.2)*(x(i)-04)-020833*(X(0-0.2)*(x®-0.4)*(x(i)-0.6)-052083*(x(i)-02)*(x(i)-04)*(x(i)-06)*(x(i)-0.8)endplot(xOyOJo'KOyO)holdonyl=sphne(x,y,xO)ploKOyl/o1)holdons=csape(x,y,1Vanational1)foplKs.'rOholdongtextC三次样条自然边界,)gtextf原图像)gtext(,4次牛顿插值)运行上述程序可知:给

17、出的(%),xx=02+0.081,1=0,1,11,10点,S(x)及P4(x)图形如下所示：问题3：3下列数据产色插手xI0I1I4I9I16I25I36I49I64y01345678可以得到平方根函数的近似，在区间0,64上作图。(1)用这9各点作8次多项式插值Ls(x).(2)用三次样条(自然边界条件)程序求S(X)。从结果看在0,64匕那个插值更精确：在区间0上，两种哪个更精确？分析：Ls(x)可由公式Ln(x)=ZyJk(Xj)得出。k=0三次样条可以利用自然边界条件。写成矩阵:4】000042心0M。M|M,M；M4其中"产一匚，心k+hjk3=U(x>17(X&

18、gt;ls(x>»句=61耳1盾,01»4n=2广。诙Xn=。k+hj解：七尸(xl)(x-4)(x9)(x-10(x230(x-64)(0-1)(0-4)(0-9)(0-16)(0-25)(0-36)(0-49)(0-64)h(x尸(x-0)(x-4)(x-9)(x-16)(x-25)(x-36)(x-49)(x-64)(1-0)(1-4)(1-9)(1-16)(L-25)(1-36)(1-49)(1-64)匕年尸(x0)(x-D炽一9)(x-16(x-25)(。-36(x-4,(。-64)"(4-0)(4-1)(4-9)(4-16)(4-25)(4-36

19、)(4-4)(4-64)卜口尸(x0)(x-l)(x-4)(x-16)(x-25)(x-36)(x-49)(x-64)3(9-0)(9-1)(9-4)(9-16)(9-25)(9-36)(9-49)(9-64)1u(x-0)(x-l)(x-4)(x-9)(x-25)(x-36)(x-49)(x-64)(16-0)6-1)(16-4)(16-9)(16-25)(16-36)(16-49)(16-64)(x-0)(x-l)(x-4)(x-9)(x-16)(x-36)(x-49)(x-64)(25-0)(25-1)(25-4)(25-9)(25-16)(25-36)(25-49)(25-64)(x-

20、O)(x-l)(x-4)(x-9)(x-16)(x-25)(x-49)(x-64)(36-0)(36-1)(36-4)(36-9)(36-16)(36-25)(36-49)(36-64)(x-0)(x-l)(x-4)(x-9)(x-16)(x-25)(x-36)(x-64)(49-0)(49-1)(49-4)(49-9)(49-16)(49-25)(49-36)(49-64)(x-O)(x-l)(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)Ls(x)=h(x)+2h(x

21、)+3b(xH4l4(x)+5b(x)+6%开？17(x)+81s(x)求三次样条插值函数由matlab计算：可得矩阵形式的线性方程组为：?00000000000o-Mo-00.2500200000.7500000000M】-1.0003750200000625000000-01000416720000058330000M3-00286000043752.000005625000-0.0119000004500200000550000M5-00061000000458320000054170乂6-00035000000046342.000005357M7-0.00220000000020000

22、一M一0,解得：Mo=0,Mi=05209,M2=00558,M3=-002614必W0006,Ms=00029,ky。0008,M-00009,0(1-x)3+-0.0868(x-0)3+0(l-x)+1.0868(x?0),xeOJ-0.02894-x)3-0.0031x-I)3+0.59384-x)+0.6388x-l),xe1,40.00199-x)3-0.0009(x-4)3+0.35359-x)+0.6271x-4),xe49S(x)=0-0.000(516-x)3+0(x-9)3+0.459016-x)-0.5708(9),xe9J6|0(25-x)3-O.OOOlx-16)3-

23、0.443625-x)+0.5600(16)xe16,250(36-x)3-0(x-25)3+0.459936-x)+0.5470(25>xe25360(48-x)3-0(x-36)3+0.463348-x)+0.5404(x.36>xe36,480(64-x)3+0(x-48)3+0.468964-x)+0.5333(48xe48,64拉格朗日插值函数与三次样条插值函数如图中所示，绿色实线条为三次样条插值曲线,蓝色虚线条为H的曲线,另外一条红色线条为拉格朗日插值曲线。图3-1为01的曲线,佟13-2为064区间上的曲线。由图3-1可以看出，红色的线条与蓝色虚线条几乎重合，所以可知

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

25、j=2:l:n-lendd(l>0d(n>0a=zeros(n,n),forj=l:l:naG，j=2,endr(l)=0,u(n)=0;forj=lln-1aG+lj)=uG+l),endb=inv(a)m=b*d,t=aperosCn-1,4),%p矩阵为S(x)函数的系数矩阵forj=l:l:n-lp(j,l)=m(j)/(6*h(j),p(j,2)=m(j+l)/(6*h0);pG,3XV(j)m忡皿2/。颂),P(j,4)y(j+l)-m(j+l)*(h(jr2/6)/h0,endp%画图形比较那个插值更精确的函数：x0=01491625364964,yO=(012345

26、678;x=064,y=lagrl(x0,y0,x),plotCxO7O/o1)holdonplot(x,y,f),holdon,pp=csape(xO,yO,ManationaF)fiipltCpp/g1),holdon,ploKxOyOJbholdon%axis(0201),%看01区间的图形时加上这条指令gtextC三次样条插值。gtextC原图像)gtextC拉格朗口插值，)%卜面是求拉格朗口插值的函数functiony=lagrl(x0>y0,x)n=lengtli(xO),m=length)，fori=lniz=x(i),s=00,fork=lnP=10,forj=lnifP

27、=P*(z-x0(j)/(x0(k)-x0(j),endends=p*yO(k)+s,endy(0=s,end问题：16观测物体的直线运动,得出以下数据:时间(t/s)009193.03950距离（s/m）010305080110求运动方程。分析：由所给的数据在坐标纸上描出如图16-1所示。从图中可以看到各点在一条直线的附近，故可以选择线性函数作拟合曲线，即令s（t）=a+bt,这里m=5,n=l。法方程矩阵为下面的形式：“（%0O）D,、（外，0O）D(f40)D).(%)dY0、,(%，%)d人a”用，%,一%的线性无关性，知道该方程存在唯一解Ax=bmm（名，外）=外外，（？，f）=N%

28、y1=11=1解：由上面的分析以及通过matlab计算得：法方程矩阵如下：614.728014.753.6378解之得：a=-78550,b=222538。由此可得运动方程为：S（t）=22.2538t-78550在matlab中拟合的曲线如图16-2所示：120100804020图16-1图16-2Matlab源程序代码：计算部分的程序代码：x=00919303950,y=010305080110,iNeros(2,2),for1=116,r(l,l)=r(l,l)+l,endfor1=116i<l,2)=r(l,2)+x(i),endfori=l161-(2,1X2,endfori=

29、ll:6r(2,2)=r(2,2)+x(ir2,enda=zeros(2,l),for1=11:6a(l,l)=a(l,l>4-y(i),a(2,l)=a(2,lHy(i)*x(i)endk=rt=ad=uiv(r),a=d*a画图的程序代码：x=00919303950,y=(010305080110,xx=0002:5.0,pp耳。lyfit(x,y,D,yy=polyvaiep,xx),ploKxyb'xx.yy)问题：18在某化学反应中，由试验得分解物浓度与时间的关系如下:时间(t/s)0510152025303540455055浓度y/Cio-4)01272.162.863443874154374514.58462464用最小二乘法求y=f(t)©分析：首先要选择拟合曲线，作出给定数据的散点图如F：G点3.323图18-1给定数据的散点图通过对散点图的分析可以看出，数据点的分布为一条单调卜.升的曲线。具有这种特性的曲线很多，我们选取如下三种函数来拟合：多项式(X)=ao+31X+.+3.m为适当选取的正整数：双