数值分析作业第一次.docx

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；分析：已知两端的一阶导数值为第一种边界条件。可写成矩阵： 其中j=， i=，dj=6fxj-1,xj,xj+1， n=1， 0=1对于第一种边界条件d0=（fx0,x1-f0），dn=（fn-fxn-1,xn解：由matlab计算得：xyhdn0.250.5000-5.52000.300.54770.05000.35711.0000-4

2、.31430.390.62450.09000.60000.6429-3.26670.450.67080.06000.42860.4000-2.42860.530.72800.08001.00000.5714-2.1150、由此得矩阵形式的线性方程组为：解得M0=-2.0286 ;M1=-1.4627;M2= -1.0333; M3= -0.8058; M4=-0.6543S(x)= matlab源程序 x=0.25 0.30 0.39 0.45 0.53; y=0.5 0.5477 0.6245 0.6708 0.7280; S0=1; S5=0.6868; for j=1:1:4 h(j)=

3、x(j+1)-x(j); end for j=2:1:4 r(j)=h(j)/(h(j-1)+h(j); end r(1)=1;for j=1:1:3 u(j)=h(j)/(h(j)+h(j+1); end u(4)=1; for j=1:1:4 f(j)=(y(j+1)-y(j)/h(j); end d(1)=6*(f(1)-S0)/h(1); d(5)=6*(S5-f(4)/h(4); for j=2:1:4 d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j); enda=zeros(5,5); for i=1:1:5 a(i,i)=2; end for i=1:1:4 a(

4、i+1,i)=u(i); a(i,i+1)=r(i); end b=inv(a); M=b*d;s=csape(x,y,complete,1 0.6868) fnplt(s,r) xlabel(x) ylabel(y)title(三次样条插值函数)plot(x,y,o,x,y,) s.coefs (2)S(0.25)=S(0.53)=0. 分析：已知两端的二阶导数只为零，可以利用自然边界条件。可写成矩阵：其中j=， i=，dj=6fxj-1,xj,xj+1， n=0=0 d0=dn=0解：由matlab计算得：xyhdn0.250.500000.300.54770.05000.35710.00

5、00-4.31470.390.62450.09000.60000.6429-3.26670.450.67080.06000.42860.4000-2.42860.530.72800.08001.00000.57140由此得矩阵形式的线性方程组为：解得M0=0 ;M1=-1.8925;M2= -0.8234; M3= -1.2108; M4=0.6054;S(x)= matlab程序：x=0.25 0.30 0.39 0.45 0.53;y=0.5 0.5477 0.6245 0.6708 0.7280;for j=1:1:4 h(j)=x(j+1)-x(j);endfor j=1:1:4 r(

6、j)=h(1)/(h(j)+h(1);endfor j=1:1:4 u(j)=1-r(j);endfor j=1:1:4 f(j)=(y(j+1)-y(j)/h(j);endfor j=1:1:4 d(j)=6*(f(1)-f(j)/(h(1)+h(j);end a=zeros(4,4);for j=1:1:4 a(j,j)=2;endfor j=1:1:3 a(j+1,j)=u(j+1); a(j,j+1)=r(j);enda(1,4)=u(1);a(4,1)=r(4);b=inv(a); M=b*d; s=csape(x,y,second,0 0) fnplt(s,r) xlabel(x)

7、 ylabel(y) title(三次样条插值函数) plot(x,y,o,x,y,) s.coefs第二种中情况时的s（x）函数的程序：x=0.25 0.30 0.39 0.45 0.53; y=0.5 0.5477 0.6245 0.6708 0.7280;for j=1:1:4 h(j)=x(j+1)-x(j); end for j=2:1:4 r(j)=h(j)/(h(j-1)+h(j); end r(1)=0; for j=1:1:3 u(j)=h(j)/(h(j)+h(j+1); end u(4)=1; for j=1:1:4 f(j)=(y(j+1)-y(j)/h(j); end

8、 d(1)=0; d(5)=0; for j=2:1:4 d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j); enda=zeros(5,5); for i=1:1:5 a(i,i)=2; end for i=1:1:4 a(i+1,i)=u(i); a(i,i+1)=r(i); end b=inv(a); M=b*d; s=csape(x,y,second,0 0) fnplt(s,r) xlabel(x) ylabel(y) title(三次样条插值函数) plot(x,y,o,x,y,) s.coefs计算实习题1已知函数在下列各点的值为xi0.20.40.6.0.81.0

9、f（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）。分析：由差商的定义及牛顿插值多项式的表示方式可知：fx0,x1，xk= Pn=f(x0)+fx0,x1(x-x0)+ fx0,x1,x2(x-x0) (x-x1)+ fx0,x1，xn(x-x0) (x-xn-1)用自然边界条件。可写成矩阵：其中j=， i=，dj=6fxj-1,xj,xj+1， n=0=0 d0=dn=0解：由matlab计算得：xif(xi)

10、h一阶差商二阶差商三阶差商四阶差商dj0.20.9800.40.920.2-0.30000.5000-3.75000.60.810.2-0.5500-0.625000.5000.500-4.50000.80.640.2-0.8500-0.75000-0.208330.50000.500-6.75001.00.380.2-1.3000-1.12500-0.62500-0.520831.00000.5000P4（x）=0.98-0.3(x-0.2)-0.625(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.

11、6)(x-0.8)解得：M0=0 ;M1=-1.5165;M2= -0.9339; M3= -1.6228; M4=0S(x)=程序：四次牛顿插值程序：x=0.2 0.4 0.6 0.8 1.0;fx=0.9 0.92 0.81 0.64 0.38;%由此函数可得差分表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); endendfor i=1:n fprintf(%4.2f,x(i); for j=1:i

12、fprintf(%10.5f,FF(i,j); end fprintf(n);endfor i=1:1:4; x=0.2 0.28 1 1.08; y(i)=0.98+0.1*(x(i)-0.2)-1.625*(x(i)-0.2)*(x(i)-0.4)+1.45833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-2.60417*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8);endplot(x,y,b,x,y,o);title(牛顿四次插值)三次牛顿插值程序：x=0.2 0.4 0.6 0.8 1.0; y=0.98 0.92 0.

13、81 0.64 0.38; for j=1:1:4 h(j)=x(j+1)-x(j); end for j=2:1:4 r(j)=h(j)/(h(j-1)+h(j); end r(1)=0; for j=1:1:3 u(j)=h(j)/(h(j)+h(j+1);%向前提 end u(4)=1; for j=1:1:4 f(j)=(y(j+1)-y(j)/h(j); end d(1)=0; d(5)=0; for j=2:1:4 d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j); enda=zeros(5,5); for i=1:1:5 a(i,i)=2; end for i=1

14、:1:4 a(i+1,i)=u(i); a(i,i+1)=r(i); end b=inv(a); M=b*d; plot(x,y,o,x,y,) s.coefsx=0.2 0.4 0.6 0.8 1.0; y=0.98 0.92 0.81 0.64 0.38;s=csape(x,y,variation)x1=0.2 ; y=-1.3393*(x1-0.4)*(x1-0.4)*(x1-0.4)-0.2464*(0.2-x1)+0.9800*(x1-0.4);a=y;x1=0.28; y=-1.3393*(x1-0.4)*(x1-0.4)*(x1-0.4)-0.2464*(0.2-x1)+0.98

15、00*(x1-0.4);b=yx1=1;y=-1.3393*(x1-0.4)*(x1-0.4)*(x1-0.4)-0.2464*(0.2-x1)+0.9800*(x1-0.4);c=yx1=1.08; y=-1.3393*(x1-0.4)*(x-0.4)*(x1-0.4)-0.2464*(0.2-x1)+0.9800*(x1-0.4);d=yx2=a b c d m=0.2 0.28 1 1.08h=fnval(m,s)plot(m,h,o)s.coefstitle(三次样条插值函数)xlabel(x)ylabel(y) 图1.13下列数据点的插值x01491625364964y0123456

16、78可以得到平方根函数的近似，在区间0,64上作图。（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

17、 l6(x)+7 l7(x)+8 l8(x)由matlab计算得：xyhdn0001110.25000-2.62504230.37500.7500-2.35429350.41670.6250-2.524316470.43750.5833-0.008425590.45000.5625-4.8037366110.45830.5500-2.7518487130.46430.5417-2.78676481500.53570 由此得矩阵形式的线性方程组为：解得:M0=0;M1=-1.1200;M2=-0.5132;M3=-1.4523;M4=1.0184;M5=-2.5065;M6=-0.4526;M7

18、=-1.2883;M8=0S(x)= 图3-1为0 64的曲线为拉格朗日插值函数与三次样条插值函数如图中所示。由图3-1可以看出，绿色的线条更靠近红色的线条，三次样条插值函数的曲线更接近函数曲线，几乎是重合的；图3-2在0 1区间，是绿色的线几乎和红色的线重合，可能是程序写的不够完美，从图上看三次样条插值的曲线接近函数曲线。由图3.2可以看出在区间0,1上，S(x)更精确。L8（x）matlab 编程0 64上程序：x1=0 1 4 9 16 25 36 49 64;y1=0 1 2 3 4 5 6 7 8; P = polyfit(x1,y1,8);%8表示8次多项式 X=0:1:64;Y

19、= polyval(P,X); plot(x1,y1,r-,X,Y,b-)title(9点8次多项式插值)xlabel(x)ylabel(y)图3-1三次插值在区间0 64的程序：x=0 1 4 9 16 25 36 49 64; y=0 1 2 3 4 5 6 7 8; for j=1:1:8 h(j)=x(j+1)-x(j); end for j=2:1:8 r(j)=h(j)/(h(j-1)+h(j); end r(1)=0;for j=1:1:7 u(j)=h(j)/(h(j)+h(j+1); end u(9)=0; for j=1:1:8 f(j)=(y(j+1)-y(j)/h(j)

20、; end d(1)=0; d(9)=0; for j=2:1:7 d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j); enda=zeros(9,9); for i=1:1:9 a(i,i)=2; end for i=1:1:4 a(i+1,i)=u(i); a(i,i+1)=r(i); end b=inv(a);M=b*d; s=csape(x,y,variational,0 0) fnplt(s,r) X=0:1:64; h=fnval(s,X) plot(x,y,r-,X,h,g*) s.coefs title(三次样条插值函数自然条件) xlabel(x)ylabel(

21、y) 图3-1L8（x）在区间0 1的程序： x1=0 1 4 9 16 25 36 49 64;y1=0 1 2 3 4 5 6 7 8;P = polyfit(x1,y1,8);%8表示8次多项式X=0:0.01:1;Y = polyval(P,X);plot(x1,y1,r-,X,Y,b-)title(9点8次多项式插值) xlabel(x) ylabel(y)三次插值在区间0 1的程序：x=0 1 4 9 16 25 36 49 64; y=0 1 2 3 4 5 6 7 8; for j=1:1:8 h(j)=x(j+1)-x(j); end for j=2:1:8 r(j)=h(j

22、)/(h(j-1)+h(j); end r(1)=0;for j=1:1:7 u(j)=h(j)/(h(j)+h(j+1); end u(9)=0; for j=1:1:8 f(j)=(y(j+1)-y(j)/h(j); end d(1)=0; d(9)=0; for j=2:1:7 d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j); enda=zeros(9,9); for i=1:1:9 a(i,i)=2; end for i=1:1:4 a(i+1,i)=u(i); a(i,i+1)=r(i); end b=inv(a);M=b*d; s=csape(x,y,variat

23、ional,0 0) fnplt(s,r) X=0:0.1:1; h=fnval(s,X) plot(x,y,r-,X,h,g*) s.coefs title(三次样条插值函数自然条件) xlabel(x)ylabel(y) 图3-2 图3-2第三章16观测物体的直线运动，得出下数据时间t/s0 0.9 1.9 3.0 3.9 5.0距离s/m0 10 30 50 80 110 求运动方程。 解：根据所给数据，在坐标纸上标出，见图16-1，从图中看到各点在一条直线附近，故可选择现行函数做拟合曲线，即令y=ax+b.这里m=5，n=1， 0=1， 1=x，故（j， k）=j（xi）k(xi),

24、(f, k)= f(xi) k=dk, k=0,1,n,法方程：Ga=d, a=(a0,a1,an)T, d=(d0, d1,dn)T,G= 由法方程得：=由matlab计算得， =； =，所以运动方程：y=-15,0002x+15.9804Matlab程序：x=0 0.9 1.9 3.0 3.9 5.0;y=0 10 30 50 80 110;plot(x,y,+)%画给出的点的图，看趋于那类型的图形，为后面设方程打基础figurex=0 0.9 1.9 3.0 3.9 5.0;y=0 10 30 50 80 110;A=polyfit(x,y,1);%拟合成二次曲线,A为返回值%提取系数a

25、=A(1);b=A(2);%画原图plot(x,y);hold on;%保存图%画拟合图plot(x,a*x+b,r);hold off;计算的matlab程序：x=0 0.9 1.9 3.0 3.9 5.0;y=0 10 30 50 80 110;t=zeros(2,2);for i=1:1:5;t(1,1)=t(1,1)+1;end for i=1:1:5 t(1,2)=t(1,2)+x(i);end for i=1:1:5 t(2,1)=t(2,1)+x(i);end for i=1:1:5 t(2,2)=t(2,2)+x(i)2;end for i=1:1:5 A=zeros(2,2)

26、; A(1,1)=A(1,1)+y(i); A(2,1)=A(2,1)+y(i)*x(i); m=A(:,1) end d=inv(t)a=d*m 图16-1由求出来的矩阵可知道拟合的方程与原方程的拟合比较得： 图16-218.在某化学反应中，由实验得分解物浓度与时间关系如下：时间t/s0 5 10 15 20 25 30 35 40 45 50 55浓度y/（*10*（-4）0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.62 4.64用最小二乘法算求y=f（t）解：根据所给数据，在坐标纸上标出，见图18-1从图中看到各点在一条直线附近，故可选择现行函数做拟合曲线，它不是线性形式，两边取对数得：y=a-b/x ; 如令Y=y，A=a，则得Y=A-b/x; =1,x.为确定A，b,先将（xi，yi）转化为（xi,Yi）.根据最小二乘法，取，得（j， k）=j（xi）k(xi), (f, k)= f(xi) k=dk, k=0,1,n,法方程：Ga=d, a=(a0,a1,an)T, d=(d0, d1,dn)T,