數(shù)值分析實(shí)驗(yàn)報(bào)告二2匯總

1、數(shù)學(xué)與信息工程學(xué)院實(shí)驗(yàn)報(bào)告課程名稱： 計(jì)算方法實(shí)驗(yàn)室：5206實(shí)驗(yàn)臺(tái)號(hào)：24班 級(jí)：姓 名：實(shí)驗(yàn)日期： 2014 年 3 月26 日實(shí)驗(yàn)名稱非線性方程求解實(shí)驗(yàn)?zāi)康暮鸵髮?shí)驗(yàn)?zāi)康模?.計(jì)算機(jī)實(shí)現(xiàn)最小二乘法2. 最小二擬合與插值法比較3. 計(jì)算機(jī)實(shí)現(xiàn)數(shù)值積分實(shí)驗(yàn)要求：1.每種算法要求達(dá)到給定的精度，輸岀數(shù)值積分結(jié)果及所需分半次數(shù)；2.利用逐次分半遞推公式；實(shí)驗(yàn)內(nèi)容和步驟：實(shí)驗(yàn)內(nèi)容：1、分別用最小二乘法，拉格朗日插值，牛頓插值，求滿足f(1)=-2,f(2)=3,f(3)=4,f(4)=-1,的三次擬合多項(xiàng)式，并在圖中畫出它們的圖像。2、地球衛(wèi)星軌道是一個(gè)橢圓，橢圓周長的計(jì)算方式是兀cs=4af屮-

2、(一)2sin2e d。，這里a是橢圓的半長軸，c是地球中心與軌道中va心(橢圓中心)的距離，記h為近地點(diǎn)距離，H為遠(yuǎn)地點(diǎn)距離，R=6371 (km為地球半徑，則a=(2R + H +h)/2,c = (H -h)/2。我國第一顆人造地球衛(wèi)星近地點(diǎn)距離h=439 ( km，遠(yuǎn)地點(diǎn)距離H=2384( k"，試求衛(wèi)星軌道的周長。第一題:最小二乘法實(shí)驗(yàn)步驟：1) 實(shí)驗(yàn)編程 建立m文件：fun cti on c d l=leasts(x,y, n)% in put x y m維% output n次多項(xiàng)式m=le ngth(x)d=;A=zeros(m, n+1); b=zeros(m,1)

3、;for k=1:n+1A(:,k)=x.A(k-1);endb=A'*y'A=A'*Ac=Abfor i=1:md(i)=c(m+1-i);enddI=poly2sym(d)運(yùn)行下列語句：x=1 2 3 4;y=-2 3 4 -1;n=3;leasts(x,y, n)2) 運(yùn)行結(jié)果d =-0.3333 -0.0000 7.3333 -9.0000l =-3002399751579999/9007199254740992*xA3-311/1125899906842624*xA2+4128299658423301/562949953421312*x-25332747903

4、96013/281474976710656拉格朗日插值 實(shí)驗(yàn)步驟：1)實(shí)驗(yàn)編程建立M文件：fun cti onp,pO=Lagra nge1(x,y,x0)syms t;n=len gth(x);p=0;for i=1:nL=1;for j=1:i-1L=L*(t-x(j)/(x(i)-x(j);en d;for j=i+1:nL=L*(t-x(j)/(x(i)-x(j);en d;L=L*y(i);p=p+L;endp=simplify(p);p0=subs(p, 't',x0)p=vpa(p)運(yùn)行下列語句：x=1 2 3 4;y=-2 3 4 -1;S=Lagra nge1

5、(x,y,1)運(yùn)行結(jié)果：p0 =-2p =-.33333333333333333333333333333333*L3+7.3333333333333333333333333333333*t-9S =-.33333333333333333333333333333333*tA3+7.3333333333333333333333333333333*t-9牛頓插值法: 實(shí)驗(yàn)步驟：1)實(shí)驗(yàn)編程建立M文件：fun ctiony,R,A,C,L=newt on (X,Y,x,M)n=le ngth(X);m=le ngth(x);for t=1:mz=x(t);A=zeros (n,n);A(:,1)=Y&

6、#39;s=0.0;p=1.0;q1=1.0;c1=1.0;for j=2: nfor i=j:nA(i,j)=(A(i,j-1)- A(i-1,j-1)/(X(i)-X(i-j+1);endq1=abs(q1*(z-X(j-1);c1=c1*j;endC=A( n,n );q仁abs(q1*(z-X( n);for k=(n-1):-1:1C=co nv(C,poly(X(k);d=le ngth(C);C(d)=C(d)+A(k,k);end y(k)= polyval(C, z);endR=M*q1/c1;L(k,:)=poly2sym(C);L=vpa(L)運(yùn)行下列語句：syms M,

7、X=1 2 3 4;Y =-2 3 4 -1; x=1;y,R,A,C,P=newt on (X,Y,x,M)2)運(yùn)行結(jié)果：P=-2.6666666666666666666666666666667*xA3+14.*xA2-18.333333333333333333333333333333*x+5.畫圖：運(yùn)行下列命令：>> x1=1 2 3 4;y仁-2 3 4 -1;x=0:0.001:3;t=0:0.001:3;y=-3002399751579999/9007199254740992*x.A3-311/1125899906842624*x42+412829965 8423301/

8、562949953421312*x-2533274790396013/281474976710656;Pl=-.33333333333333333333333333333333*t.A3+7.3333333333333333333333333333333*t -9.;P2=-.33333333333333333333333333333333*x.A3-.55511151231257827021181583404541e- 16*x.A2+7.3333333333333333333333333333333*x-9.;subplot(3,1,1)plot(x1,y1, 'r*',x

9、,y, 'b-')legend('數(shù)據(jù)點(diǎn)(xi,yi)','最小二乘法擬合曲線y=f(x)' );xlabel( 'x' );ylabel( 'y');title('數(shù)據(jù)點(diǎn)(xi,yi)和最小二乘法擬合曲線y=f(x)的圖形')subplot(3,1,2)plot(x1,y1, 'r*',t,p1,'k-')legend('數(shù)據(jù)點(diǎn)(xi,yi)','拉格朗日插值曲線 y=f(x)' );xlabel( 'x' );yla

10、bel( 'y');title('數(shù)據(jù)點(diǎn)(xi,yi)和拉格朗日插值曲線y=f(x)的圖形')subplot(3,1,3)plot(x1,y1, 'r*',x,P2, 'y-')legend('數(shù)據(jù)點(diǎn)(xi,yi)','牛頓插值曲線 y=f(x)' );xlabel( 'x' );ylabel( 'y'); title('數(shù)據(jù)點(diǎn)(xi,yi)和牛頓插值曲線y=f(x)的圖形')運(yùn)行結(jié)果：數(shù)據(jù)點(diǎn)腳.yi)和最小二乘法擬合曲線y=f(x)的團(tuán)形數(shù)擄點(diǎn)(u再

11、和拉格朗日插倩曲線y詔兇的圖形數(shù)據(jù)點(diǎn)浹觀和牛頓插值曲線嚴(yán)地)的團(tuán)形WlnrLi1111牛頓插值曲線y=f(x) >1J111 1 100.511 522.533.54實(shí)驗(yàn)結(jié)果分析：最小二乘法擬合的曲線誤差最小。也可以得到三圖合一的圖像：在以上命令的基礎(chǔ)上運(yùn)行命令 plot(x1,y1,'r*',x,y,'b-',t,p1,'k-',x,P2,'y-')得到| IIIILLLL* 00.611.S22.533 5第二題：逐次分半梯形求積公式建立f.m的m文件fun cti on f=f(x)R=6371.0;H=2384.0;

12、h=439.0;a=(2*R+H+h)/2;c=(H-h)/2;f=4*a*sqrt(1-(c/a)A2*(si n(x)A2);再建立 compT.m 的 m 文件fun cti on n ,t2=compT(a,b,eps)%利用復(fù)化梯形遞推公式求 fun在a,b的定積分; %輸入積分下限a;輸入積分上限b;精度eps； %輸出分半次數(shù)n和積分值t2 ;具體問題1的值s;h=b-a;%節(jié)點(diǎn)步長t1=(f(a)+f(b)*h/2;%梯形值 T0 :h=h/2;t2=t1/2+h*f(a+h);%經(jīng)1次分半復(fù)合梯形值err=t2-t1;淞差n=1;while abs(err)>epssu

13、m=0;h=h/2;t1=t2;n=n+1;k=1;while k<2Ansum=sum+f(a+k*h);k=k+2;endt2=t1/2+h*sumerr=t2-t1;end在程序窗口輸入：n ,t2=compT(0,pi/2,0.00001)得到結(jié)果為t2 =4.8707e+0042t2 =4.8707e+004由此我們可以得到：S =4a o'2 .仁(c)2sinJdv - 4.8707 104復(fù)合辛普森(Simpson)數(shù)值積分建立f.m的m文件fun cti on f=f(x)R=6371.0;H=2384.0;h=439.0;a=(2*R+H+h)/2;c=(H-

14、h)/2;f=4*a*sqrt(1-(c/a)A2*(si n(x)A2);建立 simps on .m 的 m文件fun cti on simps on=simps on (x0,x n, eps)S1= (x n-x0)/6*(f(x0)+4*f(x0+x n)/2)+f(x n);S2=S1;d=0;n=1;m=1;while abs(S2-S1)>eps|d=0d=1;S仁 S2;for k=0:2*m-1S2=S2+2*(-x0+x n)/(2A n*3)*f(x0+(2*k+1)*(x n-x0)/2A( n+1);endfor j=0:m-1S2=S2-(-x0+x n)/

15、(2A n*3)*f(x0+(2*j+1)*(x n-x0)/2A( n);endS2=S2-1/2*S1;n=n+1;m=2*m;S2-S1;end cs=nS2在程序窗口輸入:clearxO=O;xn=pi/2;eps=0.01;sin pus in( x0,x n, eps)運(yùn)行程序得到:cs =S2 =4.8707e+004 while err>=eps由此我們可以得到:S =4a二/21 -(£)2 sin2 日d日4= 4.8707 10龍貝格積分求解 建立 Romberg.m 的 m 文件 fun cti on R,k,T=Romberg(a,b,eps)% f積

16、分函數(shù)% a/b :積分上下限% tol :積分誤差% R : Romberg積分值% k :二分次數(shù)k=1;h=b-a;%第一步T(k,1)=h/2*(f(a)+f(b); err=1;T(k,k)= T(k,1);h=h/2;%第二步求梯形值T0temp=0;i=1;while i<2Ak temp=temp+f(a+i*h);i=i+2;endT(k+1,k)=T(k,k)/2+h*temp;%第三步計(jì)算加速值T(k+1,k+1)=(4Ak*T(k+1,k)-T(k,k)/(4Ak-1);%第四步精度控制err=abs(T(k+1,k+1)-T(k,k);k=k+1;endk=k-1;