数值分析第一次作业.doc

上传人:豆**** 文档编号:17610300 上传时间:2022-05-25 格式:DOC 页数:14 大小:308KB
返回 下载 相关 举报
数值分析第一次作业.doc_第1页
第1页 / 共14页
数值分析第一次作业.doc_第2页
第2页 / 共14页
点击查看更多>>
资源描述

《数值分析第一次作业.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

展开阅读全文
相关资源
相关搜索

当前位置:首页 > 教育专区 > 小学资料

本站为文档C TO C交易模式,本站只提供存储空间、用户上传的文档直接被用户下载,本站只是中间服务平台,本站所有文档下载所得的收益归上传人(含作者)所有。本站仅对用户上传内容的表现方式做保护处理,对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私,请立即通知得利文库网,我们立即给予删除!客服QQ:136780468 微信:18945177775 电话:18904686070

工信部备案号:黑ICP备15003705号-8 |  经营许可证:黑B2-20190332号 |   黑公网安备:91230400333293403D

© 2020-2023 www.deliwenku.com 得利文库. All Rights Reserved 黑龙江转换宝科技有限公司 

黑龙江省互联网违法和不良信息举报
举报电话:0468-3380021 邮箱:hgswwxb@163.com