(完整版)数值分析第一次作业

问题1:20给定数据如下:Xj0.250.300.390.450.53yj0.50000.54770.62450.67080.7280试求三次样条插值S(x)，并满足条件S'(0.25)=1.0000，S'(0.53)=0.6868；S'(0.25)=S'(0.53)=0。(1)是：本问题是已知五个点，由这五个点求一三次样条插值函数。边界条件有两种,已知一阶倒数，(2)是已知自然边界条件。对于第一种边界条件d°=6h°(f[X0,xl]-f0、)，dn=6hn-1(f、n-f'[Xn-l,Xn])20000M0d012100M1d102220M2d200323M3d300042M4d4h；-1hj其中j=，i=dj=6f[xj-1,xj,xj+1]，n=1，0=1hj-1hjhj-1hj对于第一种边界(已知边界的一阶倒数值)，可写出下面的矩阵方程。解：由matlab计算得:Xyhd0.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.1150由此得矩阵形式的线性方程组为:21000M0-5.52000.357120.642900M!-4.314300.600020.40000M2-3.2667000.428620.5714M3-2.428600012M4-2.1150解得Mq=-2.0286;M1=-1.4627;M2=-1.0333;M3=-0.8058;M4=-0.6546S(x)=-6.76209(0.30x)3-4.8758(x-0.25)310.0169(0.30X)10.9662(x0.25)，x[0.25,0.30]3-2.708779(0.39x)-1.9136(x-0.30)36.1075(0.39x)6.9544(x0.30)，x[0.30,0.39]-2.87040(0.45x)3-2.2384(x-0.39)310.4187(0.45X)11.188(x0.39)，x[0.39,0.45]-1.67881(0.53x)-1.3637(3x-0.45)8.3956(0.53X)9.1087(x0.45)，x[0.45,0.53]Matlab程序代码如下: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(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*d'p=zeros(n-1,4);%p矩阵为S(x)函数的系数矩阵forj=1:1:n-1p(j,1)=m(j)/(6*h(j));p(j,2)=m(j+1)/(6*h(j));P(j，3)=(y(j)-m(j)*(h(j)A2/6))/h(j);P(j，4)=(y(j+1)-m(j+1)*(h(j)A2/6))/h(j);end对于第二中边界，已知边界二阶倒数，可写出下面的矩阵：20000M0d012100M1d102220M2d200323M3d300042M4d4hi_ihj其中j=，i=-，dj6f[xj-1,Xj,Xj+1]，n=0=0d0=dn=0hj-ihihii-1hj解：由matlab计算得:xyhdn0.250.500000.300.54770.050.35710-4.31430.390.62450.090.60000.6429-3.26670.450.67080.060.42860.4000-2.42680.530.72800.0800.57140由此得矩阵形式的线性方程组为:20000M000.357120.642900M1-4.314300.600020.40000M2-3.2667000.428620.5714M3-2.428600002M40解得M0=0;M1=-1.8795;M2=-0.8636;M3=-1.0292;M4=0S(x)=30(0.30x)36.2652(x-0.25)10(0.30x)10.9697(x0.25)，x[0.25,0.30]3.4806(0.39、3x)1.5993(x-0.30)36.1137(0.39x)6.9518(x0.30)，x[0.30,0.39]2.3990(0.45\3x)2.8590(x--0.39)310.417(0.45x)11.1903(x0.39)，x[0.39,0.45]-2.1442(0.53x)-0(x0.45)3'8.3987(0.53x)9.1000(x0.45)，x[0.45,0.53]matlab程序代码如下：functiontgsanci(n,s,t)%n代表元素数,x=[0.250.300.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-1a(j+1,j)=u(j+1);a(j,j+1)=r(j);endb=inv(a)k=am=b*d'p=zeros(n-1,4);%p矩阵为S(x)函数的系数矩阵forj=1:1:n-1p(j,1)=m(j)/(6*h(j));p(j,2)=m(j+1)/(6*h(j));P(j，3)=(y(j)-m(j)*(h(j)A2/6))/h(j);P(j，4)=(y(j+1)-m(j+1)*(h(jF2/6))/h(j);endP由下面的一段程序，画出自然边界与第一边界的图形：程序代码如下：x=[0.250.300.390.450.53];y=[0.50.54770.62450.67080.7280];s=csaPe(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图中的三条线几乎重合。yy问题2:1已知函数在下列各点的值为Xi0.20.40.6.0.81.0f(Xi)0.980.920.810.640.38试用4次牛顿插值多项式P4（x）及三次样条函数S（x）（自然边界条件）对数据进行插值。用图给出｛（Xi,yi），xi=0.2+0.08i，i=0，1,11,10｝，P4（x）及S（x）。分析：（1）要用4次牛顿插值多项式处理数据，Pn=f（x0）+f[x0,X1]（X_X0）+f[x0,X1,X2]（X_X0）（X-X1）+…+f[x0,X1，…Xn]（X_X0）…（X-Xn_1）首先我们要知道牛顿插值多项式的系数，即均差表中得部分均差。解：由matlab计算得：Xif(Xi)一阶差商二阶差商三阶差商四阶差商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-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)%由此函数可得差分表n=length(x);fprintf('*****************差分表*****************************\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');end2)用三次样条插值函数由上题分析知,要求各点的M值：有matlab编程计算求出个三次样条函数：20000M000.50020.50000M1-3.750000.50020.5000M2-4.5000000.50020.500M3-6.750000002M40解得：M0=0;M1=-1.6071;M2=-1.0714;M3=-3.1071;M4=0[0.4,0.6][0.6,0.8](00.4x)31.3393(x0.2)4.90(00.4x)4.6536(x0.2)，x[0.2,0.4]-1.339(30.6x)30.8929(x-0.4)34.653(60.6x)4.085(7x0.4)，x-0.8929(0.8x)3-2.589(3x-0.6)34.085(70.8x)3.303(6x0.6)，x-2.589(31.0x)3-(0x0.8)33.3036(1.0x)1.9(x0.8)，x[0.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)+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*d'p=zeros(n-1,4);%p矩阵为S(x)函数的系数矩阵forj=1:1:n-1p(j,1)=m(j)/(6*h(j));p(j,2)=m(j+1)/(6*h(j));P(j，3)=(y(j)-m(j)*(h(j)A2/6))/h(j);P(j，4)=(y(j+1)-m(j+1)*(h(j)A2/6))/h(j);endP下面是画牛顿插值以及三次样条插值图形的程序：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=[011011]x0=0.2+0.08*kfori=1:1:4y0(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)endplot(xO,yO,'o',xO,yO)holdony1=spline(x,y,xO)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)及P4(x)图形如下所示：问题3:3下列数据点的插值x01491625364964y012345678可以得到平方根函数的近似，在区间［0,64］上作图。（1）用这9各点作8次多项式插值L8（x）.（2）用三次样条（自然边界条件）程序求S（x）。从结果看在［0,64］上，那个插值更精确；在区间［0,1］上，两种哪个更精确？n分析：L8（X）可由公式Ln（x）=yJk（Xj）得出。k0三次样条可以利用自然边界条件。写成矩阵:L8(x)=11(x)+2l2(x)+3l3(x)+4l4(x)+5l5(x)+6l6(x)+7l7(x)+818(x)求三次样条插值函数由matlab计算：可得矩阵形式的线性方程组为：2.000000000000M000.25002.00000.7500000000M11.000.37502.00000.625000000M20.1000.41672.00000.58330000M30.02860000.43752.00000.5625000M40.011900000.45002.00000.550000m50.0061000000.45832.00000.54170M60.00350000000.46342.00000.5357M70.0022000000002.0000M80解得：其中解：lo（x）=|l（x）=|2（X）=l3（X）=l4（X）=l5（X）=l6（X）=l7（X）=l8（X）=20000M。d012100M1d102220M2d200323M3d300042m4d4hj-ihjdj=6f[xj-1,xj,xj+1]，0=0d0=dn=0jhj-ihj，i=hj-1hj，n=(X-1)(x4)(X9)(x16)(x25)(x36)(x49)(X64)(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)(x1)(x9)(x16)(x25)(x36)(x49)(X64)(40)(41)(49)(416)(425)(436)(449)(464)(x-0)(x1)(x4)(x16)(x25)(x36)(x49)(X64)(90)(91)(94)(916)(925)(936)(949)(964)(X-0)(x1)(x4)(X9)(x:25)(x36)(x49)(x64)(160)(161)(164)(169)(1625)(1636)(1649)(16l54)(x-0)(x-1)(x4)(x9)(x16)(x36)(x49)(x64)(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)(x4)(X9)(x16)(x25)(x36)(x64)(490)(491)(494)(499)(4916)(4925)(4936)(4964)(x-0)(x-1)(x4)(x9)(x16)(x25)(x36)(x49)(640)(641)(644)(649)(6416)(6425)(6436)(6449)Mo=O;Mi=-O.52O9;M2=0.0558;M3=-0.0261;M4=0.0006;M5=-0.0029;M6=-0.0008;M7=--0.0009;M8=0-20-0.0006(16x)30(x9)30.4590(16x)-0.5708(x-9)，x[9,16]0(25x)3-0.0001(x16)3-0.4436(25x)0.5600(x-16)，x[16,25]0(36x)3-0(x25)30.4599(36x)0.5470(x-25)，x[25,36]0(48x)3-0(x36)30.4633(48x)0.5404(x-36)，x[36,48]0(64x)30(x48)30.4689(64x)0.5333(x-48)，x[48,64]0(1X)30.0868(x0)30(1x)1.0868(x-0),x[0,1]-0.0289(4x)3-0.0031(x-1)30.5938(4x)0.6388(x1),x[1,4]0.0019(9x)30.0009(x-4)30.3535(9x)0.627(x4),x[4,9]S(x)=0拉格朗日插值函数与三次样条插值函数如图中所示，绿色实线条为三次样条插值曲线，蓝色虚线条为x=y2的曲线，另外一条红色线条为拉格朗日插值曲线。图3-1为［01］的曲线,图3-2为［064］区间上的曲线。由图3-1可以看出，红色的线条与蓝色虚线条几乎重合，所以可知拉格朗日插值函数的曲线更接近开平方根的函数曲线，在［01］朗格朗日插值更精确。而在区间［064］上从图3-2中可以看出蓝色虚线条和绿色线条是几乎重合的，而红色线条在［3070］之间有很大的振荡，所在在区间［064］三次样条插值更精确写。图3-110080604020原图像102030拉格朗日插值三-次样条插值40506070图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-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*d't=ap=zeros(n-1,4);%p矩阵为S(x)函数的系数矩阵forj=1:1:n-1p(j,1)=m(j)/(6*h(j));p(j,2)=m(j+1)/(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)A2/6))/h(j);endP%画图形比较那个插值更精确的函数：x0=[01491625364964];yO=[O12345678];x=0:64;y=lagr1(x0,y0,x);plot(xO,yO,'o')holdonplot(x,y,'r');holdon;pp=csape(xO,yO,'variational')fnplt(pp,'g');holdon;%看［01］区间的图形时加上这条指令plot(xO,yO,':b');holdon%axis([0201]);gtext('三次样条插值')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-xO(j))/(xO(k)-xO(j));endends=p*yO(k)+s;endy(i)=s;end问题：16.观测物体的直线运动，得出以下数据:时间（t/s）00.91.93.03.95.0距离（s/m）010305080110求运动方程。分析：由所给的数据在坐标纸上描出如图16-1所示。从图中可以看到各点在一条直线的附近，故可以选择线性函数作拟合曲线，即令s（t）=a+bt,这里m=5,n=1。法方程矩阵为下面的形式：0,0D0,nDa0f,0Dn,0Dn,nDanf,nD0,1,的线性无关性，知道该方程存在唯一解Axb(j,k)jiki)i1(j,jiyii1解：由上面的分析以及通过法方程矩阵如下：matlab计算得:614.7a28014.753.63b1078解之得：a=-7.8550,b=22.2538。由此可得运动方程为：S(t)=22.2538t-7.8550.在matlab中拟合的曲线如图16-2所示：源数据点最小二乘拟合曲线120120a100100--80r"y?q80)60卜/\60Q4040Q20r/r200f/r0Cr-20Iarott图16-2图16-1Matlab源程序代码：计算部分的程序代码:x=[00.91.93.03.95.0];y=[010305080110];r=zeros(2,2);fori=1:1:6;r(1,1)=r(1,1)+1;endfori=1:1:6r(1,2)=r(1,2)+x(i);endfori=1:1:6r(2,1)=r(2,1)+x(i);endfori=1:1:6r(2,2)=r(2,2)+x(i)A2;enda=zeros(2,1);fori=1:1:6a(1,1)=a(1,1)+y(i);a(2,1)=a(2,1)+y(i)*x(i)endk=rt=ad=inv(r);a=d*a画图的程序代码：x=[00.91.93.03.95.0];y=[010305080110];xx=0:0.02:5.0;给定点0102030405060tpp=polyfit(x,y,1);yy=polyval(pp,xx);plot(x,y,'o',xx,yy)问题：18在某化学反应中，由试验得分解物浓度与时间的关系如下:时间(t/s)0510152025303540455055浓度y/(10-4)01.272.162.863.443.874.154.374.514.584.624.64用最小二乘法求y=f(t)。分析：首先要选择拟合曲线，作出给定数据的散点图如下:54.543.532.521.510.50图18-1给定数据的散点图通过对散点图的分析可以看出，数据点的分布为一条单调上升的曲线。具有这种特性的曲线很多，我们选取如下三种函数来拟合：多项式(x)=a0a1X…amxm,m为适当选取的正整数；x(X)axbb(x)aex(a>0,b<0)。在matlab中分别用上述拟合函数拟合曲线，本题应采用倒指数拟合，即y(t)=a*exp(bt-1)拟合曲线。对y(t)=a*exp(bt-1)两边取对数有Iny=Ina+b/t;如令Y=Iny,A=Ina,则得Y=A+b/t;为确定A,b,先将(ti,yi)转化为(ti,Yi).根据最小二乘法，取0(t)1，1(t)1/t。用lsqcurvefit函数拟合曲线，程序代码如下：创建M文件：functionF=mf(a,t)F=a(1)*exp(a(2).*t.A(-1));在命令窗口中敲入如下代码：t=[5:5:55]y=[1.272.162.863.443.874.154.374.514.584.624.64];a0=[0.5,0.5];a=lsqcurvefit('pf',a0,t,y)回车后可得：a(1)=5.5031,a(2)=-8.7872下面的计算问题用matlab编程实现：x=[510152025303540455055];y0=[1.272.162.863.443.874.154.374.514.584.624.64];y=log(y0)y(1)=0r=zeros(2,2);fori=1:1:12;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)+1/x(i)A2;enda=zeros(2,1);fori=2:1:12a(1,1)=a(1,1)+y(i);a(2,1)=a(2,1)+y(i)*(1/x(i));endk=rt=ad=inv(r);m=d*a-4.8922/to由matlab软件计算得：a=3.9864,b=-4.8922。所以y(t)=3.9864e有下面的语句给出拟合后的图形：x=[0510152025303540455055];y0=[01.272.162.863.443.874.154.374.514.584.624.64];fplot('5.5031*exp((-8.7872)*(x).A(-1))',[0,55]);holdonplot（x,yO,'o',x,yO,'r'）holdongtext（'倒指数拟合曲线'）gtext（'原图像'）5图18-2