数值分析第一次作业

1、精品文档(1)是已知一阶倒数, 对于第种边界(已知边界的一阶倒数值),口与出卜面的矩阵方程。一 2九0000M 0 1-d 02九100M 1d 10N 22九 20M 2=d 200N 32九3M 3d 3一°00N 42-M 4 1:d 4(2)是已知自然边界条件。,hj-1其中j=hj-1 hji=hj-1hj hjdj=6fx j-1 凶,xj+1 ,对于第一种边界条件d0=9 h。(fx0,x1-f0、), dn=h(f n-f'Xn-i,Xn)n-1问题1:20.给定数据如下表:xj0.250.300.390.450.53yj0.50000.54770.62450

2、.67080.7280试求三次样条插值S(x),并满足条件(1)S'(0.25)=1.0000,S'(0.53)=0.6868;(2)S'(0.25)=S'(0.53)=0。分析：本问题是已知五个点，由这五个点求一三次样条插值函数。边界条件有两种,精品文档解：由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.000

3、0.5714-2.1150由此得矩阵形式的线性方程组为：一21000M0I-5.5200-0.357120.642900M1-4.314300.600020.40000M2-3.2667000.428620.5714M3-2.4286-00012M41i-2.1150解得M0=-2.0286;M1=-1.4627;M2=-1.0333;M3=-0.8058；M4=-0.6546S(x)=,x 0.25,0.30,x 0.30,0.39x 0.39,0.45x 0.45,0.53-6.76209(0.30-x)3-4.8758(x-0.25)3+10.0169(0.30x)+10.9662(x0

4、.25)-2.708779(0.39-x)3-1.9136(x-0.30)3+6.1075(0.39-x)+6.9544(x0.30)- 2.87040 (0.45 -x)-1.67881 (0.53 -x)-2.2384(x-0.39)3+10.4187(0.45x)+11.188(x0.39),3-1.3637(x-0.45)+8.3956(0.53x)+9.1087(x0.45),Matlab程序代码如下:functiontgsanci(n,s,t)%n代表元素数，s,t代表端点的一阶导。x=0.250.300.390.450.53;y=0.50.54770.62450.67080.72

5、80;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

6、=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);PO,3)=(yO)-m(j)*(hO)A2/6)/h(j);PO,4)=(yO+1)-mO+1)*(h(j)A2/6)/h(j);end对于第二中边界，已知边界二阶倒数，可写出下面的矩阵：'2九0000M0d°1得2%00M1d10出2九20M2二d200%2%M3d3000P-42M41d

7、4hj-1.hi一一其中出=,九i=,dj6fxj-1,xj,xj+1,Hn=九0=0d0=dn=0hj-1+hjhj-1+hj解：由matlab计算得:xyhy九dn0.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由此得矩阵形式的线性方程组为：一20000"M000.357120.642900M1-4.314300.600020.40000M2-3.2667000.428620.

8、5714M3-2.4286-00002M4110解得Mo=0;M1=-1.8795;M2=-0.8636;M3=-1.0292;M4=0S(x)=0.25,0.300.30) , x 0.30,0.39-0.39) , x 0.39,0.45x 0.45,0.53r,、3,、3,、0(0.30-x)-6.2652(x-0.25)+10(0.30-x)+10.9697(x-0.25),xgj3.4806(0.39-x)3-1.5993(x-0.30)3+6.1137(0.39x)+6.9518(x-2.3990(0.45-x)3-2.8590(x-0.39)3+10.417(0.45-x)+11

9、.1903(x1-2.1442(0.53-x)3-0(x-0.45)3+8.3987(0.53-x)+9.1000(x-0.45),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);end

10、forj=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);PO,3)=(y(j)-mO)*(hO)A2/6)/h(j);PO,4)=(y(

11、j+1)-mO+1)*(hO)A2/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,

12、'complete',1.00.6868)fnplt(r,'b')gtext('三次样条第一边界')holdonr.coefs图中的三条线几乎重合。y0.0.5450.0.5350.0.5250.0.5150.0.505图20-1在一小段区间内的图形图20-2在整个给出的区间的图形问题2:1已知函数在下列各点的值为xi0.20.40.6.0.81.0f(xi)0.980.920.810.640.38试用4次牛顿插值多项式P4(x)及三次样条函数S(x)(自然边界条件)对数据进行插值。用图给出(Xi,y。，Xi=0.2+0.08i,i=0,1,11

13、,10,P4(x)及S(x)。分析：(1)要用4次牛顿插值多项式处理数据，Pn=f(x0)+fx0,x1(x-x0)+fx0,x1,x2(x-x0)(x-x1)+fx0,x，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所以有四次插值

14、牛顿多项式为：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;仅=0.980.920.810.640.38;newtonchzh(x,fx);functionnewtonchzh(x,fx)%由此函数可得差分表n=length(x);fprintf('*差分表

15、*n');FF=ones(n,n);FF(:,1)=fx'fori=2:nforj=i:nFF0,i)=(FF0,i-1)-FF0-1,i-1)/(x0)-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一MoI10.50020.50000M1-3.750000.50020.5000m2-4.50000

16、00.50020.500M3-6.750000002_M4一110-解得：Mo=0；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(x0.2),x0.2,0.4j-1.3393(0.6x)30.8929(x-0.4)3+4.6536(0.6x)+4.0857(x0.4),x=0.4,0.6-0.8929(0.8x)3-2.5893(x-0.6)3+4.0857(0.8x)+3.3036(x0.6),xw0.6,0.8-2.5893(1.0-x)3-0(x0.8)3+3.3036(1

17、.0x)+1.9(x0.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)+

18、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);PO,3)=(y(j)-mO)*(hO)A2/6)/h(j);PO,4)=(yO+1)-mO+1)*(hO)A2/6)/h(j);endp下面是画牛顿插值以及三次样条

19、插值图形的程序：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*kfo门=1:1:4y0(i)=0.98-0.3*(x(i)-0.2)-0.62500*(x(i)-0.2)*(

20、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,yO,'o',x0,y0)holdony1=spline(x,y,x0)plot(x0,y1,'o')holdons=csape(x,y,'variational')fnplt(s,'r')holdongtext('三次样条自然边界)gtext('原图像')gtext('4次牛

21、顿插值)运行上述程序可知：给出的(xi,%),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)=£ydk(xj)得出。k=0一2九0000一M0-fd出27M00M1d02八20M2二d00%2，-3M3d10002M4一1d写成

22、矩阵:01234三次样条可以利用自然边界条件。其中解：l0(x)=l1(x)=l2(x)=l3(x)=l4(x)=l5(x)=l6(x)=l7(x)=l8(x)=,hj-ihjNj=,九i=,dj=6fxji,Xj,Xj+i,Nn=儿0=0do=dn=0hj-ihjhj-ihj(X-1)(x4)(x9)(x16)(x25)(x36)(x49)(x64)(01)(。4)(09)(016)(025)(036)(049)(064)(x-0)(x4)(x9)(x16)(x25)(x36)(x49)(x64)(1-0)(1-4)(1-9)(1-16)(1-25)(1-36)(1-49)(1-64)(x-

23、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)(x-1)(x-4)(x-16)(x-25)(x-36)(x-49)(x-64)(9-0)(9-1)(9-4)(9-16)(9-25)(9-36)(9-49)(9-64)(x-0)(x-1)(x-4)(x-9)(x-25)(x-36)(x-49)(x-64)(160)(161)(164)(169)(1625)(1636)(1649)(1664)(x-0)(x-1)(x4)(x9)(x16)(x36)(x49)(x

24、64)(25-0)(25-1)(25-4)(25-9)(25-16)(25-36)(25-49)(25-64)(x-0)(x1)(x4)(x9)(x16)(x25)(x49)(x64)(36-0)(36-1)(36-4)(36-9)(36-16)(36-25)(36-49)(36-64)(x-0)(x-1)(x-4)(x-9)(x-16)(x-25)(x-36)(x-64)(490)(491)(494)(499)(4916)(4925)(4936)(4964)(x-0)(x-1)(x4)(x9)(x16)(x25)(x36)(x49)(64-0)(64-1)(64-4)(64-9)(64-16

25、)(64-25)(64-36)(64-49)L8(X)=11(x)+2l2(x)+3l3(x)+4l4(x)+5l5(x)+6l6(x)+7l7(x)+8l8(X)求三次样条插值函数由matlab计算：可得矩阵形式的线性方程组为：2.000000000000I-m000.25002.00000.7500000000M1-1.000.37502.00000.625000000M2-0.1000.41672.00000.58330000M3-0.02860000.43752.00000.5625000M4=-0.011900000.45002.00000.550000M5-0.0061000000

26、.45832.00000.54170M6-0.00350000000.46342.00000.5357M7-0.0022000000002.0000_M80解得：M0=0;M1=-0.5209;M2=0.0558；M3=-0.0261;M4=0.0006;M5=-0.0029;M6=-0.0008；M7=-0.0009;M8=010(1x)3+0.0868(x0)3+0(1x)+1.0868(x-0),xw0,1-0.0289(4x)3-0.0031(x-1)3+0.5938(4x)+0.6388(x1),xw1,40.0019(9x)30.0009(x-4)3+0.3535(9-x)+0.6

27、271(x-4),x4,9S(x)=x 9,16x 16,25GJ-0.0006(16-x)3+0(x9)3+0.4590(16x)-0.5708(x-9),0(25-x)3-0.0001(x16)3-0.4436(25x)+0.5600(x-16),0(36x)3-0(x25)3+0.4599(36x)+0.5470(x-25),x25,360(48x)3-0(x36)3+0.4633(48x)+0.5404(x-36),xw36,48064_x)3+0(x48)3+0.4689(64_x)+0.5333(x-48),xw48,64拉格朗日插值函数与三次样条插值函数如图中所示，绿色实线条为三

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

29、1625364964;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);

30、endb=inv(a)m=b*d't=ap=zeros(n-1,4);%p矩阵为S(x)函数的系数矩阵forj=1:1:n-1p(j,1)=mO)/(6*h(j);pQ,2)=mQ+1)/(6*hO);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%画图形比较那个插值更精确的函数：x0=01491625364964;y0=012345678;x=0:64;y=lagr1(x0,y0,x);plot(x0,y0,'o')holdonplot(x,y,'r

31、9;);holdon;pp=csape(x0,y0,'variational')fnplt(pp,'g');holdon;plot(x0,y0,':b');holdon%axis(0201);%看01区间的图形时加上这条指令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;

32、forj=1:nifj=kp=p*(z-xO(j)/(x0(k)-x0(j);endends=p*y0(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。法方程矩阵为下面的形式：(九,%、9，平。卜3+.3I3鼻dn,%b(邛nnb/laj(fWn)D邛0,叼，邛n的线性无关性，知道该方程存在唯一解Ax=bm/

33、中,cp=ycp.cpL./<P.(j,k)jiki,(ji=1解：由上面的分析以及通过matlab计算得：法方程矩阵如下：614.7a_28014.753.63b|J078解之得：a=-7.8550,b=22.2538。由此可得运动方程为：S(t)=22.2538t-7.8550.在matlab中拟合的曲线如图16-2所示：源数据点源数据点120.10080X6040200J10246tmf)八 iyii =1图 16-2图16-1Matlab源程序代码：计算部分的程序代码：x=00.91.93.03.95.0;y=010305080110;r=zeros(2,2);fo门=1:1:6

34、;r(1,1)=r(1,1)+1;endfo门=1:1:6r(1,2)=r(1,2)+x(i);endfo门=1:1:6r(2,1)=r(2,1)+x(i);endfo门=1:1:6r(2,2)=r(2,2)+x(i)A2;enda=zeros(2,1);fo门=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;pp=polyfit(x,y,1);yy=polyval(pp,xx);给 定 点4.53.52.51.50.5102030405060t图18-1给定数据的散点图通过对散点图的分析可以看出，数据点的分布为一条单调上升的曲线。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。分析：首先要选择拟合曲线，作出给定数据的散点图如下: