《数值分析第一次作业.doc》由会员分享,可在线阅读,更多相关《数值分析第一次作业.doc(17页珍藏版)》请在三一文库上搜索。
1、内科大数值分析作业班级:研5班 姓名:王雪萍 学号:201102208 专业:双控-问题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)= M
3、atlab程序代码如下: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*(
4、f(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-
7、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)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 o
9、nr.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-
11、0.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+
12、1); 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.9
15、8-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上从图
19、3-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); %看
22、0 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)0103050
23、80110求运动方程。分析:由所给的数据在坐标纸上描出如图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
24、50 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=
27、5: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)+
28、1;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