數(shù)值分析報(bào)告作業(yè)第一次

1、第二章 習(xí)題20.給定數(shù)據(jù)如下表：xj0.250.300.390.450.53yj0.50000.54770.62450.67080.7280試求三次樣條插值S(x)，并滿足條件(1)S(0.25)=1.0000，S(0.53)=0.6868；分析：已知兩端的一階導(dǎo)數(shù)值為第一種邊界條件?？蓪懗删仃嚕?其中j=，i=，dj=6fxj-1,xj,xj+1， n=1，0=1對于第一種邊界條件d0=（fx0,x1-f0），dn=（fn-fxn-1,xn解：由matlab計(jì)算得：xyhdn0.250.5000-5.52000.300.54770.05000.35711.0000-4.31430.390.

2、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)=x(j+1)-x(j

3、); 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(i+1,i)=u(i

4、); 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(三次樣條插值函數(shù))plot(x,y,o,x,y,) s.coefs (2)S(0.25)=S(0.53)=0. 分析：已知兩端的二階導(dǎo)數(shù)只為零，可以利用自然邊界條件?？蓪懗删仃嚕浩渲衘=，i=，dj=6fxj-1,xj,xj+1，n=0=0 d0=dn=0解：由matlab計(jì)算得：xyhdn0.250.500000.300.54770.05000.35710.0000-4.31470.3

5、90.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(j)=h(1)/(h(j

6、)+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) ylabel(y) t

7、itle(三次樣條插值函數(shù)) plot(x,y,o,x,y,) s.coefs第二種中情況時(shí)的s（x）函數(shù)的程序：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 d(1)=0; d(5

8、)=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(三次樣條插值函數(shù)) plot(x,y,o,x,y,) s.coefs計(jì)算實(shí)習(xí)題1已知函數(shù)在下列各點(diǎn)的值為xi0.20.40.6.0.81.0f（xi）0.980.9

9、20.810.640.38試用4次牛頓插值多項(xiàng)式P4（x）及三次樣條函數(shù)S（x）（自然邊界條件）對數(shù)據(jù)進(jìn)行插值。用圖給出（xi，yi），xi=0.2+0.08i，i=0，1, 11, 10，P4（x）及S（x）。分析：由差商的定義及牛頓插值多項(xiàng)式的表示方式可知：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)用自然邊界條件?？蓪懗删仃嚕浩渲衘=，i=，dj=6fxj-1,xj,xj+1，n=0=0 d0=dn=0解：由matlab計(jì)算得：xif(xi)h一階差商二階差商三階差商四階

10、差商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.6)(x-0.8)解得：M0=

11、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;%由此函數(shù)可得差分表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 fprintf(%10.5f

12、,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.81 0.64 0.38;

13、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:1:4 a(i+1,i)=

14、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.9800*(x1-0.4);b=

15、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(三次樣條插值函數(shù))xlabel(x)ylabel(y) 圖1.13下列數(shù)據(jù)點(diǎn)的插值x01491625364964y012345678可以得到平方根函數(shù)的近似

16、，在區(qū)間0,64上作圖。（1）用這9各點(diǎn)作8次多項(xiàng)式插值L8(x).（2）用三次樣條（自然邊界條件）程序求S（x）。從結(jié)果看在0,64上，那個(gè)插值更精確；在區(qū)間0,1上，兩種哪個(gè)更精確？分析：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

17、l8(x)由matlab計(jì)算得：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=-1.2883;M8=0S(x)

18、= 圖3-1為0 64的曲線為拉格朗日插值函數(shù)與三次樣條插值函數(shù)如圖中所示。由圖3-1可以看出，綠色的線條更靠近紅色的線條，三次樣條插值函數(shù)的曲線更接近函數(shù)曲線，幾乎是重合的；圖3-2在0 1區(qū)間，是綠色的線幾乎和紅色的線重合，可能是程序?qū)懙牟粔蛲昝溃瑥膱D上看三次樣條插值的曲線接近函數(shù)曲線。由圖3.2可以看出在區(qū)間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次多項(xiàng)式 X=0:1:64;Y = polyval(P,X); p

19、lot(x1,y1,r-,X,Y,b-)title(9點(diǎn)8次多項(xiàng)式插值)xlabel(x)ylabel(y)圖3-1三次插值在區(qū)間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); end d(1)=0; d(9

20、)=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(三次樣條插值函數(shù)自然條件) xlabel(x)ylabel(y) 圖3-1L8（x）在區(qū)間0

21、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次多項(xiàng)式X=0:0.01:1;Y = polyval(P,X);plot(x1,y1,r-,X,Y,b-)title(9點(diǎn)8次多項(xiàng)式插值) xlabel(x) ylabel(y)三次插值在區(qū)間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)/(h(j-1)+h(j); e

22、nd 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,variational,0 0) fnplt(

23、s,r) X=0:0.1:1; h=fnval(s,X) plot(x,y,r-,X,h,g*) s.coefs title(三次樣條插值函數(shù)自然條件) xlabel(x)ylabel(y) 圖3-2 圖3-2第三章16觀測物體的直線運(yùn)動，得出下數(shù)據(jù)時(shí)間t/s0 0.9 1.9 3.0 3.9 5.0距離s/m0 10 30 50 80 110 求運(yùn)動方程。 解：根據(jù)所給數(shù)據(jù)，在坐標(biāo)紙上標(biāo)出，見圖16-1，從圖中看到各點(diǎn)在一條直線附近，故可選擇現(xiàn)行函數(shù)做擬合曲線，即令y=ax+b.這里m=5，n=1，0=1，1=x，故（j，k）=j（xi）k(xi), (f, k)= f(xi) k=dk,

24、k=0,1,n,法方程：Ga=d, a=(a0,a1,an)T, d=(d0, d1,dn)T,G= 由法方程得：=由matlab計(jì)算得，=；=，所以運(yùn)動方程：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,+)%畫給出的點(diǎn)的圖，看趨于那類型的圖形，為后面設(shè)方程打基礎(chǔ)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為返回值%提取系數(shù)a=A(1);b=A(2);%畫原圖plot(

25、x,y);hold on;%保存圖%畫擬合圖plot(x,a*x+b,r);hold off;計(jì)算的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); 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.在某化學(xué)反應(yīng)中，由實(shí)驗(yàn)得分解物濃度與時(shí)間關(guān)系如下：時(shí)間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）解：根據(jù)所給數(shù)據(jù)，在坐標(biāo)紙上標(biāo)出，見圖18-