数值分析第一次作业.doc

《数值分析第一次作业.doc》由会员分享，可在线阅读，更多相关《数值分析第一次作业.doc（14页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、【精品文档】如有侵权，请联系网站删除，仅供学习与交流数值分析第一次作业.精品文档.问题1：20.给定数据如下表：xj0.250.300.390.450.53yj0.50000.54770.62450.67080.7280试求三次样条插值S(x)，并满足条件(1)S(0.25)=1.0000，S(0.53)=0.6868；（2）S(0.25)=S(0.53)=0。分析：本问题是已知五个点，由这五个点求一三次样条插值函数。边界条件有两种，（1）是已知一阶倒数，（2）是已知自然边界条件。对于第一种边界(已知边界的一阶倒数值)，可写出下面的矩阵方程。其中j=，i=，dj=6fxj-1,xj,xj+1，

2、 n=1，0=1对于第一种边界条件d0=（fx0,x1-f0），dn=（fn-fxn-1,xn）解：由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由此得矩阵形式的线性方程组为：解得 M0=-2.0286;M1=-1.4627;M2= -1.0333; M3= -0.8058; M4=-0.6546S(x)= Ma

3、tlab程序代码如下：function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。x=0.25 0.30 0.39 0.45 0.53;y=0.5 0.5477 0.6245 0.6708 0.7280;n=5,s=1.0,t=0.6868for j=1:1:n-1 h(j)=x(j+1)-x(j);endfor j=2:1:n-1 r(j)=h(j)/(h(j)+h(j-1);endfor j=1:1:n-1 u(j)=1-r(j);endfor j=1:1:n-1 f(j)=(y(j+1)-y(j)/h(j);endfor j=2:1:n-1 d(j)=6*(f

4、(j)-f(j-1)/(h(j-1)+h(j);end d(1)=6*(f(1)-s)/h(1) d(n)=6*(t-f(n-1)/h(n-1) a=zeros(n,n);for j=1:1:n a(j,j)=2;end r(1)=1; u(n)=1;for j=1:1:n-1 a(j+1,j)=u(j+1); a(j,j+1)=r(j);endb=inv(a)m=b*dp=zeros(n-1,4); %p矩阵为S(x)函数的系数矩阵for j=1:1:n-1 p(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)

5、2/6)/h(j); p(j,4)=(y(j+1)-m(j+1)*(h(j)2/6)/h(j);end对于第二中边界，已知边界二阶倒数，可写出下面的矩阵：其中j=，i=，dj=6fxj-1,xj,xj+1，n=0=0 d0=dn=0解：由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由此得矩阵形式的线性方程组为：解得M0=0 ;M1= -1.8795;M2=

6、-0.8636; M3= -1.0292; M4=0S(x)= matlab程序代码如下：function tgsanci(n,s,t) %n代表元素数， x=0.25 0.30 0.39 0.45 0.53; y=0.5 0.5477 0.6245 0.6708 0.7280; n=5for j=1:1:n-1 h(j)=x(j+1)-x(j);endfor j=2:1:n-1 r(j)=h(j)/(h(j)+h(j-1);endfor j=1:1:n-1 u(j)=1-r(j);endfor j=1:1:n-1 f(j)=(y(j+1)-y(j)/h(j);endfor j=2:1:n-1

7、 d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j);end d(1)=0 d(n)=0 a=zeros(n,n);for j=1:1:n a(j,j)=2;end r(1)=0; u(n)=0;for j=1:1:n-1 a(j+1,j)=u(j+1); a(j,j+1)=r(j);endb=inv(a)k=am=b*dp=zeros(n-1,4); %p矩阵为S(x)函数的系数矩阵for j=1:1:n-1 p(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)2/6)/h(j); p(j,4)=

8、(y(j+1)-m(j+1)*(h(j)2/6)/h(j);end p由下面的一段程序，画出自然边界与第一边界的图形：程序代码如下：x=0.25 0.30 0.39 0.45 0.53; y=0.5 0.5477 0.6245 0.6708 0.7280; s=csape(x,y,variational) fnplt(s,r)hold onxlabel(x)ylabel(y)gtext(三次样条自然边界)plot(x,y,o,x,y,:g)hold ons.coefsr=csape(x,y,complete,1.0 0.6868)fnplt(r,b)gtext(三次样条第一边界)hold on

9、r.coefs图中的三条线几乎重合。 图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，yi），xi=0.2+0.08i，i=0，1, 11, 10，P4（x）及S（x）。分析：（1）要用4次牛顿插值多项式处理数据，Pn=f(x0)+fx0,x1(x-x0)+ fx0,x1,x2(x-x0) (x-x1)+ fx0,x1，xn(x-x0) (x-xn-1)首

10、先我们要知道牛顿插值多项式的系数，即均差表中得部分均差。解：由matlab计算得：xi f(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

11、.4)(x-0.6)（x-0.8）计算牛顿插值中多项式系数的程序如下：function varargout=newtonliu(varargin)clear,clcx=0.2 0.4 0.6 0.8 1.0;fx=0.98 0.92 0.81 0.64 0.38;newtonchzh(x,fx);function newtonchzh(x,fx)%由此函数可得差分表n=length(x);fprintf(*差分表*n);FF=ones(n,n);FF(:,1)=fx;for i=2:n for j=i:n FF(j,i)=(FF(j,i-1)-FF(j-1,i-1)/(x(j)-x(j-i+1

12、); endendfor i=1:n fprintf(%4.2f,x(i); for j=1:i fprintf(%10.5f,FF(i,j); end fprintf(n);end（2）用三次样条插值函数由上题分析知,要求各点的M值： 有matlab编程计算求出个三次样条函数：解得：M0=0 ;M1= -1.6071;M2= -1.0714; M3= -3.1071; M4=0三次样条插值函数计算的程序如下：function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。x=0.2 0.4 0.6 0.8 1.0;y=0.98 0.92 0.81 0.64 0.38;

13、 n=5for j=1:1:n-1 h(j)=x(j+1)-x(j);endfor j=2:1:n-1 r(j)=h(j)/(h(j)+h(j-1);endfor j=1:1:n-1 u(j)=1-r(j);endfor j=1:1:n-1 f(j)=(y(j+1)-y(j)/h(j);endfor j=2:1:n-1 d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j);end d(1)=0 d(n)=0 a=zeros(n,n);for j=1:1:n a(j,j)=2;end r(1)=0; u(n)=0;for j=1:1:n-1 a(j+1,j)=u(j+1); a(j,

14、j+1)=r(j);endb=inv(a)m=b*dp=zeros(n-1,4); %p矩阵为S(x)函数的系数矩阵for j=1:1:n-1 p(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)2/6)/h(j); p(j,4)=(y(j+1)-m(j+1)*(h(j)2/6)/h(j);end p下面是画牛顿插值以及三次样条插值图形的程序：x=0.2 0.4 0.6 0.8 1.0;y=0.98 0.92 0.81 0.64 0.38;plot(x,y)hold on for i=1:1:5y(i)=0.98

15、-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=0 1 10 11x0=0.2+0.08*kfor i=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)-

16、0.4)*(x(i)-0.6)*(x(i)-0.8)endplot( x0,y0,o,x0,y0 )hold ony1=spline(x,y,x0)plot(x0,y1,o)hold ons=csape(x,y,variational) fnplt(s,r)hold ongtext(三次样条自然边界)gtext(原图像)gtext(4次牛顿插值)运行上述程序可知：给出的（xi，yi），xi=0.2+0.08i，i=0，1, 11, 10点，S（x）及P4（x）图形如下所示：问题3 ：3下列数据点的插值x01491625364964y012345678可以得到平方根函数的近似，在区间0,64上作

17、图。（1）用这9各点作8次多项式插值L8(x).（2）用三次样条（自然边界条件）程序求S（x）。从结果看在0,64上，那个插值更精确；在区间0,1上，两种哪个更精确？分析：L8(x)可由公式Ln(x)=得出。 三次样条可以利用自然边界条件。写成矩阵：其中j=，i=，dj=6fxj-1,xj,xj+1，n=0=0 d0=dn=0解：l0(x)=l1(x)= l2(x)= l3(x)= l4(x)= l5(x)= l6(x)= l7(x)= l8(x)= L8(x)= l1(x)+2 l2(x)+3 l3(x)+4 l4(x)+5 l5(x)+6 l6(x)+7 l7(x)+8 l8(x)求三次样

18、条插值函数由matlab计算：可得矩阵形式的线性方程组为：解得:M0=0;M1=-0.5209;M2=0.0558;M3=-0.0261;M4=0.0006;M5=-0.0029;M6=-0.0008;M7=-0.0009;M8=0S(x)= 拉格朗日插值函数与三次样条插值函数如图中所示，绿色实线条为三次样条插值曲线，蓝色虚线条为x=y2的曲线，另外一条红色线条为拉格朗日插值曲线。图3-1为0 1的曲线，图3-2为0 64区间上的曲线。由图3-1可以看出，红色的线条与蓝色虚线条几乎重合，所以可知拉格朗日插值函数的曲线更接近开平方根的函数曲线，在0 1朗格朗日插值更精确。而在区间0 64上从图3

19、-2中可以看出蓝色虚线条和绿色线条是几乎重合的，而红色线条在30 70之间有很大的振荡，所在在区间0 64三次样条插值更精确写。 图3-1 图3-2Matlab程序代码如下：求三次样条插值函数的程序：function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。y=0 1 2 3 4 5 6 7 8;x=0 1 4 9 16 25 36 49 64; n=9for j=1:1:n-1 h(j)=x(j+1)-x(j);endfor j=2:1:n-1 r(j)=h(j)/(h(j)+h(j-1);endfor j=1:1:n-1 u(j)=1-r(j);endfor

20、j=1:1:n-1 f(j)=(y(j+1)-y(j)/h(j);endfor j=2:1:n-1 d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j);end d(1)=0 d(n)=0 a=zeros(n,n);for j=1:1:n a(j,j)=2;end r(1)=0; u(n)=0;for j=1:1:n-1 a(j+1,j)=u(j+1); a(j,j+1)=r(j);endb=inv(a)m=b*dt=ap=zeros(n-1,4); %p矩阵为S(x)函数的系数矩阵for j=1:1:n-1 p(j,1)=m(j)/(6*h(j); p(j,2)=m(j+1)/(

21、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);end p%画图形比较那个插值更精确的函数：x0=0 1 4 9 16 25 36 49 64;y0=0 1 2 3 4 5 6 7 8;x=0:64;y=lagr1(x0,y0,x);plot(x0,y0,o)hold onplot(x,y,r);hold on;pp=csape(x0,y0,variational)fnplt(pp,g);hold on;plot(x0,y0,:b);hold on%axis(0 2 0 1); %看0

22、 1区间的图形时加上这条指令gtext(三次样条插值)gtext(原图像)gtext(拉格朗日插值)%下面是求拉格朗日插值的函数function y=lagr1(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问题：16.观测物体的直线运动，得出以下数据：时间（t/s）00.91.93.03.95.0距离（s/m）01030508

23、0110求运动方程。分析：由所给的数据在坐标纸上描出如图16-1所示。从图中可以看到各点在一条直线的附近，故可以选择线性函数作拟合曲线，即令s(t)=a+bt ,这里m=5,n=1。法方程矩阵为下面的形式：的线性无关性，知道该方程存在唯一解解：由上面的分析以及通过matlab计算得： 法方程矩阵如下：解之得：a=-7.8550 ,b=22.2538 。由此可得运动方程为：S(t)=22.2538t-7.8550 .在matlab中拟合的曲线如图16-2所示： 图16-1 图16-2Matlab源程序代码：计算部分的程序代码：x=0 0.9 1.9 3.0 3.9 5.0;y=0 10 30 5

24、0 80 110;r=zeros(2,2);for i=1:1:6; r(1,1)=r(1,1)+1;end for i=1:1:6 r(1,2)=r(1,2)+x(i); end for i=1:1:6 r(2,1)=r(2,1)+x(i); end for i=1:1:6 r(2,2)=r(2,2)+x(i)2; end a=zeros(2,1); for i=1:1:6 a(1,1)=a(1,1)+y(i); a(2,1)=a(2,1)+y(i)*x(i) end k=r t=a d=inv(r);a=d*a 画图的程序代码：x=0 0.9 1.9 3.0 3.9 5.0;y=0 10

25、30 50 80 110;xx=0:0.02:5.0;pp=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)。分析：首先要选择拟合曲线，作出给定数据的散点图如下： 图18-1 给定数据的散点图通过对散点图的分析可以看出，数据点的分布为一条单调上升的曲线。具有这种特性的曲线很多，我们选取如下三

26、种函数来拟合： 多项式 j (x) = a0 + a1x + + amxm，m为适当选取的正整数； ； （a 0, b 0）。在matlab中分别用上述拟合函数拟合曲线，本题应采用倒指数拟合，即y(t)=a*exp(bt-1)拟合曲线。对y(t)=a*exp(bt-1)两边取对数有 y=a+b/t ; 如令Y=y，A=a，则得Y=A+b/t;为确定A，b,先将（ti，yi）转化为（ti,Yi）.根据最小二乘法，取，。用lsqcurvefit函数拟合曲线，程序代码如下：创建M文件：function F=mf(a,t)F=a(1)*exp(a(2).*t.(-1);在命令窗口中敲入如下代码：t=5

27、:5:55y= 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.62 4.64;a0=0.5,0.5;a=lsqcurvefit(pf,a0,t,y)回车后可得：a(1)= 5.5031 ,a(2)= -8.7872下面的计算问题用matlab编程实现：x=5 10 15 20 25 30 35 40 45 50 55 ;y0=1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.62 4.64;y=log(y0)y(1)=0r=zeros(2,2);for i=1:1:12; r(1,1)=r(1,1)+1

28、;end for i=2:1:12 r(1,2)=r(1,2)+1/x(i); end for i=2:1:12 r(2,1)=r(2,1)+1/x(i); end for i=2:1:12 r(2,2)=r(2,2)+1/x(i)2; end a=zeros(2,1); for i=2:1:12 a(1,1)=a(1,1)+y(i); a(2,1)=a(2,1)+y(i)*(1/x(i) ; end k=r t=a d=inv(r);m=d*a由matlab软件计算得：a= 3.9864,b= -4.8922。所以y(t)= 3.9864e-4.8922/t。有下面的语句给出拟合后的图形：x=0 5 10 15 20 25 30 35 40 45 50 55 ;y0=0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.62 4.64;fplot(5.5031*exp(-8.7872)*(x).(-1),0,55);hold onplot(x,y0,o,x,y0,r)hold ongtext(倒指数拟合曲线)gtext(原图像) 图18-2