數(shù)值分析第一次作業(yè)

1、問題1：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；（2）S(0.25)=S(0.53)=0。分析：本問題是已知五個點，由這五個點求一三次樣條插值函數(shù)。邊界條件有兩種，（1）是已知一階倒數(shù)，（2）是已知自然邊界條件。對于第一種邊界(已知邊界的一階倒數(shù)值)，可寫出下面的矩陣方程。其中j=，i=，dj=6fxj-1,xj,xj+1， n=1，0=1對于第一種邊界條件d0=（fx0,x1-f0），dn=（fn-fx

2、n-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)= Matlab程序代碼如下：function tgsanci(n,s,t) %n代表元素

3、數(shù)，s,t代表端點的一階導(dǎo)。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(j)-f(j-1)/(h(j-1)+h(j);end d(1)=6*(f(1)-

4、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*d'p=zeros(n-1,4); %p矩陣為S(x)函數(shù)的系數(shù)矩陣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(

5、j)2/6)/h(j);end對于第二中邊界，已知邊界二階倒數(shù)，可寫出下面的矩陣：其中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= -0.8636; M3= -1.0292; M4=0S(x)= matl

6、ab程序代碼如下：function tgsanci(n,s,t) %n代表元素數(shù)， 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 d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j);en

7、d 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*d'p=zeros(n-1,4); %p矩陣為S(x)函數(shù)的系數(shù)矩陣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);e

8、nd 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

9、',1.0 0.6868)fnplt(r,'b')gtext('三次樣條第一邊界')hold onr.coefs圖中的三條線幾乎重合。 圖20-1 在一小段區(qū)間內(nèi)的圖形 圖20-2 在整個給出的區(qū)間的圖形問題2：1已知函數(shù)在下列各點的值為xi0.20.40.6.0.81.0f（xi）0.980.920.810.640.38試用4次牛頓插值多項式P4（x）及三次樣條函數(shù)S（x）（自然邊界條件）對數(shù)據(jù)進行插值。用圖給出（xi，yi），xi=0.2+0.08i，i=0，1, 11, 10，P4（x）及S（x）。分析：（1）要用4次牛頓插值多項式處理數(shù)據(jù)，Pn=

10、f(x0)+fx0,x1(x-x0)+ fx0,x1,x2(x-x0) (x-x1)+···+ fx0,x1，···xn(x-x0) ···(x-xn-1)首先我們要知道牛頓插值多項式的系數(shù)，即均差表中得部分均差。解：由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-

11、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.4)(x-0.6)（x-0.8）計算牛頓插值中多項式系數(shù)的程序如下：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)%由

12、此函數(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',FF(i,j); end fprintf('n');end（2）用三次樣條插值函數(shù)由上題分析知,要求各點的M值： 有matlab編程計算求出

13、個三次樣條函數(shù)：解得：M0=0 ;M1= -1.6071;M2= -1.0714; M3= -3.1071; M4=0三次樣條插值函數(shù)計算的程序如下：function tgsanci(n,s,t) %n代表元素數(shù)，s,t代表端點的一階導(dǎo)。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)-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

14、)-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*d'p=zeros(n-1,4); %p矩陣為S(x)函數(shù)的系數(shù)矩陣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

15、)-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-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

16、)*(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)-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,&

17、#39;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下列數(shù)據(jù)點的插值x01491625364964y012345678可以得到平方根函數(shù)的近似，在區(qū)間0,64上作圖。（1）用這9各點作8次多項式插值L8(x).（2）用三次樣條（自然邊界條件）程序求S（x）。從結(jié)果看在0,64

18、上，那個插值更精確；在區(qū)間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)求三次樣條插值函數(shù)由matlab計算：可得矩陣形式的線性方程組為：解得:M0=0;M1=-0.5209;M2=0.0558

19、;M3=-0.0261;M4=0.0006;M5=-0.0029;M6=-0.0008;M7=-0.0009;M8=0S(x)= 拉格朗日插值函數(shù)與三次樣條插值函數(shù)如圖中所示，綠色實線條為三次樣條插值曲線，藍色虛線條為x=y2的曲線，另外一條紅色線條為拉格朗日插值曲線。圖3-1為0 1的曲線，圖3-2為0 64區(qū)間上的曲線。由圖3-1可以看出，紅色的線條與藍色虛線條幾乎重合，所以可知拉格朗日插值函數(shù)的曲線更接近開平方根的函數(shù)曲線，在0 1朗格朗日插值更精確。而在區(qū)間0 64上從圖3-2中可以看出藍色虛線條和綠色線條是幾乎重合的，而紅色線條在30 70之間有很大的振蕩，所在在區(qū)間0 64三次樣條

20、插值更精確寫。 圖3-1 圖3-2Matlab程序代碼如下：求三次樣條插值函數(shù)的程序：function tgsanci(n,s,t) %n代表元素數(shù)，s,t代表端點的一階導(dǎo)。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 j=1:1:n-1 f(j)=(y(j+1)-y(j)/h(j);endfor j=2:1:n-1 d(j)=6*

21、(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*d't=ap=zeros(n-1,4); %p矩陣為S(x)函數(shù)的系數(shù)矩陣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(

22、j+1)-m(j+1)*(h(j)2/6)/h(j);end p%畫圖形比較那個插值更精確的函數(shù)：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 1

23、區(qū)間的圖形時加上這條指令gtext('三次樣條插值')gtext('原圖像')gtext('拉格朗日插值')%下面是求拉格朗日插值的函數(shù)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.觀測物體的直線運動，得出以下數(shù)據(jù)：時間（t/s）00.

24、91.93.03.95.0距離（s/m）010305080110求運動方程。分析：由所給的數(shù)據(jù)在坐標紙上描出如圖16-1所示。從圖中可以看到各點在一條直線的附近，故可以選擇線性函數(shù)作擬合曲線，即令s(t)=a+bt ,這里m=5,n=1。法方程矩陣為下面的形式：的線性無關(guān)性，知道該方程存在唯一解解：由上面的分析以及通過matlab計算得： 法方程矩陣如下：解之得：a=-7.8550 ,b=22.2538 。由此可得運動方程為：S(t)=22.2538t-7.8550 .在matlab中擬合的曲線如圖16-2所示： 圖16-1 圖16-2Matlab源程序代碼：計算部分的程序代碼：x=0 0.9

25、 1.9 3.0 3.9 5.0;y=0 10 30 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

26、 0.9 1.9 3.0 3.9 5.0;y=0 10 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在某化學(xué)反應(yīng)中，由試驗得分解物濃度與時間的關(guān)系如下：時間(t/s)0510152025303540455055濃度y/(10-4)01.272.162.863.443.874.154.374.514.584.624.64用最小二乘法求y=f(t)。分析：首先要選擇擬合曲線，作出給定數(shù)據(jù)的散點圖如下： 圖18-1 給定數(shù)據(jù)的散點圖通過對散點圖的分析可以看

27、出，數(shù)據(jù)點的分布為一條單調(diào)上升的曲線。具有這種特性的曲線很多，我們選取如下三種函數(shù)來擬合： 多項式 j (x) = a0 + a1x + + amxm，m為適當選取的正整數(shù)； ； （a >0, b <0）。在matlab中分別用上述擬合函數(shù)擬合曲線，本題應(yīng)采用倒指數(shù)擬合，即y(t)=a*exp(bt-1)擬合曲線。對y(t)=a*exp(bt-1)兩邊取對數(shù)有 y=a+b/t ; 如令Y=y，A=a，則得Y=A+b/t;為確定A，b,先將（ti，yi）轉(zhuǎn)化為（ti,Yi）.根據(jù)最小二乘法，取，。用lsqcurvefit函數(shù)擬合曲線，程序代碼如下：創(chuàng)建M文件：function F=mf(a,t)