數(shù)值分析計算實習(xí)題_第1頁
數(shù)值分析計算實習(xí)題_第2頁
數(shù)值分析計算實習(xí)題_第3頁
數(shù)值分析計算實習(xí)題_第4頁
數(shù)值分析計算實習(xí)題_第5頁
已閱讀5頁,還剩21頁未讀, 繼續(xù)免費閱讀

下載本文檔

版權(quán)說明:本文檔由用戶提供并上傳,收益歸屬內(nèi)容提供方,若內(nèi)容存在侵權(quán),請進行舉報或認(rèn)領(lǐng)

文檔簡介

1、插值法下列數(shù)據(jù)點的插值x01491625364964y012345678可以得到平方根函數(shù)的近似,在區(qū)間0,64上作圖.用這9個點作8次多項式插值Ls(x).用三次樣條(第一邊界條件)程序求S(x).從得到結(jié)果看在0,64上,哪個插值更精確;在區(qū)間0,1上,兩種插值哪個更精確?解:(1)拉格朗日插值多項式,求解程序如下symsxl;x1=01491625364964;y1=012345678;n=length(x1);Ls=sym(0);fori=1:nl=sym(y1(i);fork=1:i-1l=l*(x-x1(k)/(x1(i)-x1(k);endfork=i+1:nl=l*(x-x1(

2、k)/(x1(i)-x1(k);endLs=Ls+l;endLs=simplify(Ls)%為所求插值多項式Ls(x).輸出結(jié)果為Ls=-24221063/63504000*x2+95549/72072*x-1/3048192000*x8-2168879/435456000*x4+19/283046400*x7+657859/10886400*x3+33983/152409600*x5T3003/2395008000*x飛(2)三次樣條插值,程序如下x1=01491625364964;y1=012345678;x2=0:1:64;y3=spline(x1,y1,x2);p=polyfit(x2

3、,y3,3);%得到三次樣條擬合函數(shù)S=p(l)+p(2)*x+p(3)*x2+p(4)*x3%得到S(x)輸出結(jié)果為:S=2288075067923491/73786976294838206464-2399112304472833/576460752303423488*x+4552380473376713/18014398509481984*x2+999337332656867/1125899906842624*x3在區(qū)間0,64上,分別對這兩種插值和標(biāo)準(zhǔn)函數(shù)作圖,plot(x2,sqrt(x2),b,x2,y2,r,x2,y3,y)藍色曲線為y=Vx函數(shù)曲線,紅色曲線為拉格朗日插值函數(shù)曲線

4、,黃色曲線為三次樣條插值曲線100806040200100806040200可以看到藍色曲線與黃色曲線幾乎重合,因此在區(qū)間0,64三次樣條插值更精確。在0,1區(qū)間上由上圖看不出差別,不妨代入幾組數(shù)據(jù)進行比較,取x4二0:0.2:1x4=0:0.2:1;sqrt(x4)subs(Ls,x,x4)spline(x1,y1,x4)運行結(jié)果為%準(zhǔn)確值%拉格朗日插值%三次樣條插值ans=00.44720.63250.77460.89441.0000ans=00.25040.47300.67060.84551.0000ans=00.24290.46300.66170.84031.0000從這幾組數(shù)值上可以

5、看出在0,1區(qū)間上,拉格朗日插值更精確。數(shù)據(jù)擬合和最佳平方逼近有實驗給出數(shù)據(jù)表x0.00.10.20.30.50.81.0y1.00.410.500.610.912.022.46試求3次、4次多項式的曲線擬合,再根據(jù)數(shù)據(jù)曲線形狀,求一個另外函數(shù)的擬合曲線,用圖示數(shù)據(jù)曲線及相應(yīng)的三種擬合曲線。解:(1)三次擬合,程序如下symx;x0=0.00.10.20.30.50.81.0;y0=1.00.410.500.610.912.022.46;cc=polyfit(x0,y0,3);S3二cc(l)+cc(2)*x+cc(3)*x2+cc(4)*x3%三次擬合多項式xx=x0(1):0.1:x0(l

6、ength(x0);yy=polyval(cc,xx);plot(xx,yy,-);holdon;plot(x0,y0,x);xlabel(x);ylabel(y);運行結(jié)果S3=-7455778416425083/1125899906842624+1803512222945437/140737488355328*x-655705280524945/140737488355328*x2+4172976892199509/4503599627370496*x3圖像如下(2)4次多項式擬合symx;x0=0.00.10.20.30.50.81.0;y0=1.00.410.500.610.912.0

7、22.46;cc=polyfit(x0,y0,4);S3二cc(l)+cc(2)*x+cc(3)*x2+cc(4)*x3+cc(5)*x4xx=x0(1):0.1:x0(length(x0);yy=polyval(cc,xx);plot(xx,yy,r);holdon;plot(x0,y0,x);xlabel(x);ylabel(y);運行結(jié)果S3=3248542900396215/1125899906842624-3471944732519151/281474976710656*x+4580931990070637/281474976710656*x2-5965836931688425/n2

8、5899906842624*x3+8491275464650307/9007199254740992*x4圖像如下2.521.5y2.521.5y10.5x(3)另一個擬合曲線,新建一個M-file程序如下:functionC,L=lagran(x,y)w=length(x);n=w-1;L=zeros(w,w);fork=1:n+1V=1;forj=1:n+1ifk=jV=conv(V,poly(x(j)/(x(k)-x(j);endendL(k,:)=V;endC=y*L在命令窗口中輸入以下的命令:x=0.00.10.20.30.50.81.0;y=1.00.410.500.610.912

9、.022.46;cc=polyfit(x,y,4);xx=x(1):0.1:x(length(x);yy=polyval(cc,xx);plot(xx,yy,r);holdon;plot(x,y,x);xlabel(x);ylabel(y);x=0.00.10.20.30.50.81.0;y=0.940.580.470.521.002.002.46;%y中的值是根據(jù)上面兩種擬合曲線而得到的估計數(shù)據(jù),不是真實數(shù)據(jù)C,L=lagran(x,y);xx=0:0.01:1.0;yy=polyval(C,xx);holdon;plot(xx,yy,b,x,y,.);圖像如下x解線性方程組的直接解法線性方

10、程組Ax=b的A及b為1078732A=75651078732A=7565,b=2386109337591031則解1111.用MATLAB內(nèi)部函數(shù)求detA及A的所有特征值和cond(A)2若令1078.17.2A+6A=708q504()65,求解(A+SA)(x+6x)=b,輸出向量x和85.989.8996.99599.98|6x|2從理論結(jié)果和實際計算兩方面分析線性方程組Ax=b解得相對誤差l|6x|2/|x|2及a的相對誤差|6a|2/|a|2的關(guān)系.解:(1)程序如下clear;A=10787;7565;86109;75910;det(A)cond(A,2)eig(A)輸出結(jié)果為

11、ans=1ans=2.9841e+003ans=0.01020.84313.858130.2887(2)程序如下A=1078.17.2;7.085.0465;85.989.899;6.99599.98;b=32233331;x0=1111;x=Ab%擾動后方程組的解x1=x-x0norm(x1,2)x1=x-x0norm(x1,2)%8x的值%5x的2-范數(shù)運行結(jié)果為x=-9.586318.3741-3.22583.5240 x1=-10.586317.3741-4.22582.5240ans=20.9322程序如下A=1078.17.2;7.085.0465;85.989.899;6.995

12、99.98;A0=10787;7565;86109;75910;b=32233331;x0=1111;x=Ab;x1=x-x0;norm(x1,2);A1=A-AO;%SA的值norm(xl,2)/norm(x0,2)%|x|2/|x|2的值norm(Al,2)/norm(A0,2)%|A|2/|A|2的值輸出結(jié)果為ans=10.4661ans=0.0076可見A相對誤差只為0.0076,而得到的結(jié)果x的相對誤差就達到了10.4661,該方程組是病態(tài)的,A的條件數(shù)為2984.1遠遠大于1,當(dāng)A只有很小的誤差就會給結(jié)果帶來很大的影響。非線性方程數(shù)值解法求下列方程的實根x八2-3x+2-e八x=0

13、;x八3+2x八2+10 x-20=0.要求:(1)設(shè)計一種不動點迭代法,要使迭代序列收斂,然后再用斯特芬森加速迭代,計算到|x(k)-x(k-1)|10,-8)為止。(2)用牛頓迭代,同樣計算到|x(k)-x(k-1)|10,-8)。輸出迭代初值及各次迭代值和迭代次數(shù)k,比較方法的優(yōu)劣。解:(1)先用畫圖的方法估計根的范圍ezplot(,x“2_3*x+2_exp(x),);gridon;可以估計到方程的根在區(qū)間(0,1);選取迭代初值為x0=0.5;構(gòu)造不動點迭代公式x(k+1)=(x(k廠2+2-e八x(k)/3;W(x)二(x八2+2-e八x)/3;當(dāng)0 x1時,0W(x)1;曠(x)

14、、1;因此迭代公式收斂。0.257534913615250.25753491361525程序如下:formatlong;f二inline(x2+2-exp(x)/3)disp(x=);x=feval(f,0.5);disp(x);Eps=1E-8;i=1;while1x0=x;i=i+1;x=feval(f,x);disp(x);ifabs(x-x0)Epsbreak;endendi,x運行結(jié)果為f=Inlinefunction:f(x)=(x2+2-exp(x)/3x=0.200426243099960.272749065098370.253607156584130.258550376264

15、940.257265636335090.257598985162190.257512454514830.257529084167960.257530285439860.257530285439860.257530597238330.257530204510460.257530306445640.257530279987670.25753028685501i=14x=0.25753028685501斯特芬森加速法,程序如下:formatlong;f二inline(x-(x2+2-exp(x)/3-x)2/(x2+2-exp(x)/3)2+2-exp(x2+2-exp(x)/3)/3-2*(x2+

16、2-exp(x)/3+x);disp(x=);x=feval(f,0.5);disp(x);Eps=1E-8;i=1;while1x0=x;i=i+1;x=feval(f,x);disp(x);ifabs(x-x0)Epsbreak;endendi,x運行結(jié)果為x=0.258684427565790.257530317719810.25753028543986i=4i=40.25753028543986i=40.25753028543986i=4x=0.25753028543986牛頓迭代法,程序如下:formatlong;x=sym(x);f二sym(,x“2-3*x+2-exp(x);df

17、=diff(f,x);FX=x-f/df;Fx=inline(FX);disp(,x=,);x1=0.5;disp(x1);Eps=1E-8;i=0;while1x0=x1;i=i+1;x1=feval(Fx,x1);disp(x1);ifabs(x1-x0)1E10break;endifabs(x-x0)Epsbreak;endendi,x運行結(jié)果:f=Inlinefunction:f(x)=(-2*x2-10*x+20廠1/3x=0.166666666666676.09259259259259-38.38843164151806-8.478196837919431e+002-4.76366

18、0785374071e+005-1.512815059604763e+011i=6x=-1.512815059604763e+011迭代6次后x的值大得令人吃驚,表明構(gòu)造的式子并不收斂.也無法構(gòu)造出收斂的不動點公式牛頓迭代法,程序如下:formatlong;x=sym(x);f=sym(x3+2*x2+10*x-20);df=diff(f,x);FX=x-f/df;Fx=inline(FX);disp(x=);x1=0.5;disp(x1);Eps=1E-8;i=0;while1x0=x1;i=i+1;x1=feval(Fx,x1);disp(x1);ifabs(x1-x0)Epsbreak;

19、endendi,x1運行結(jié)果:x=1.500000000000001.373626373626371.368814819623961.368808107834411.36880810782137i=4x1=1.36880810782137比較三種方法,牛頓法的收斂性比較好,相比不動點迭代法要構(gòu)造出收斂的公式比較難,牛頓法迭代次數(shù)也較少,收斂速度快,只是對初值的要求很高,幾種方法各有利弊,具體采用哪種也需因題而異。常微分方程初值問題數(shù)值解法給定初值問題y=-50y+50 x八2+2x,0 x1;y(0)=1/3;用經(jīng)典的四階R-K方法解該問題,步長分別取h=0.1,0.025,0.01,計算并打

20、印x=0.1i(i=0,l,,10)各點的值,與準(zhǔn)確值y(x)=l/3e八(-50 x)+x八2比較。解:取步長h=0.1,程序如下:%經(jīng)典的四階R-K方法clear;F二-50*y+50*x2+2*x;a=0;b=1;h=0.1;n=(b-a)/0.1;X=a:0.1:b;Y=zeros(1,n+1);Y(1)=1/3;fori=1:nx=X(i);y=Y(i);K1=h*eval(F);x=x+h/2;y=y+K1/2;K2=h*eval(F);y=Y(i)+K2/2;K3=h*eval(F);x=X(i)+h;y=Y(i)+K3;K4=h*eval(F);Y(i+1)=Y(i)+(K1+

21、2*K2+2*K3+K4)/6;end%準(zhǔn)確值temp=;f二dsolve(Dy=-50*y+50*x2+2*x,y(0)=l/3,x);df=zeros(1,n+1);fori=1:n+1temp=subs(f,x,X(i);df(i)=double(vpa(temp);enddisp(步長四階經(jīng)典R-K法準(zhǔn)確值)disp(X,Y,df);運行結(jié)果:步長四階經(jīng)典R-K法準(zhǔn)確值1.0e+010*00.000000000010000.0000000000200000.000000000010000.000000000020000.000000000030000.000000000040000.0

22、00000000050000.000000000060000.000000000070000.000000000080000.000000000090000.00000000010000%畫圖觀察結(jié)果figure;plot(X,df,k*,X,Y,gridon;0.000000000033330.000000000460550.000000006306250.000000086404940.000001184363000.000016235451100.000222560671340.003050935427780.041823239217400.573326903478097.8593563

23、00837710.000000000033330.000000000001220.000000000004000.000000000009000.000000000016000.000000000025000.000000000036000.000000000049000.000000000064000.000000000081000.00000000010000titletitle(四階經(jīng)典R-K法解常微分方程);legend(準(zhǔn)確值,四階經(jīng)典R-K法);當(dāng)X值接近1的時候,偏離準(zhǔn)確值太大。當(dāng)步長h=0.025時,將上面程序中的h改為0.025即可,運行結(jié)果:步長00.10000000000

24、000步長00.100000000000000.200000000000000.300000000000000.400000000000000.500000000000000.600000000000000.700000000000000.800000000000000.900000000000001.00000000000000圖像如下:四階經(jīng)典R-K法0.333333333333330.103135240342880.044285273275990.051967957555070.093957311494390.160345314357570.248085701305570.356241884727580.484525906616270.6

溫馨提示

  • 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
  • 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯(lián)系上傳者。文件的所有權(quán)益歸上傳用戶所有。
  • 3. 本站RAR壓縮包中若帶圖紙,網(wǎng)頁內(nèi)容里面會有圖紙預(yù)覽,若沒有圖紙預(yù)覽就沒有圖紙。
  • 4. 未經(jīng)權(quán)益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
  • 5. 人人文庫網(wǎng)僅提供信息存儲空間,僅對用戶上傳內(nèi)容的表現(xiàn)方式做保護處理,對用戶上傳分享的文檔內(nèi)容本身不做任何修改或編輯,并不能對任何下載內(nèi)容負(fù)責(zé)。
  • 6. 下載文件中如有侵權(quán)或不適當(dāng)內(nèi)容,請與我們聯(lián)系,我們立即糾正。
  • 7. 本站不保證下載資源的準(zhǔn)確性、安全性和完整性, 同時也不承擔(dān)用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。

評論

0/150

提交評論