《数值分析第一次作业.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; n=5for j=1:1:n-1 h(j)=x(j+1)
