數(shù)值分析自主上機(jī)題

1、2016級數(shù)值分析上機(jī)實(shí)踐報(bào)告機(jī)械工程學(xué)院2016021910吳臻標(biāo)準(zhǔn)題：迭代格式的比較3設(shè)方程 f(x)=x - 3x -=0 有三個(gè)實(shí)根 x1 =1.8793 , x2 =-0.34727 ,x3=-1.53209現(xiàn)采用下面三種不同計(jì)算格式，求f(x)=0的根x1或x23x 121、x = xX312、x =33、x = 33x 1數(shù)學(xué)原理：簡單迭代法是根據(jù)f(x)=0這個(gè)方程，對其進(jìn)行等價(jià)變換化為x=書(x)并由該 式產(chǎn)生逼近解x*的迭代數(shù)列xk，這就是簡單迭代法的思想。對于同一 f(x)=0可 以等價(jià)變換出不同的x=書(x)而且它們各自的收斂性不同。程序設(shè)計(jì)：#in clude<

2、iostream>#in clude<cmath>#in clude<cstdlib> using n amespace std;double f(double i)/外調(diào)函數(shù)f(x),每次更新新的函數(shù)/以第一種迭代方式為例子double k,m,sum;k=3*i+1;m=pow(i,2.0);sum=k/m; return sum;int mai n()double x,x0;int N;最大迭代次數(shù)int k;cout«"輸入初解:"cin> >x0;coutvv"輸入最大迭代次數(shù):"cin>

3、;>N;for(k=1;k<=N;k+)x=f(xO);if(fabs(x-x0)<0.0000001)cout«"迭代次數(shù):"vvkvvendl; cout«"輸出得到的解:"<<x<<endl; system("pause");return 0;else x0=x;coutvv"已達(dá)到最大迭代次數(shù):"<<N<<endl; cout«"輸出得到的解:"<<x<<endl; s

4、ystem("pause");return 0;實(shí)驗(yàn)結(jié)果：揄入最大迭代枕數(shù) 陸代枚數(shù)輸出得到的解= -1-53209 陳袴任意毎靈續(xù)愉人初解：-乩* 輸入嘉大迭代枕數(shù)迓 已達(dá)到最犬迭代次數(shù)応回 輸出得到的解:-輸入初解= -9-3 輸入最大迭代枚數(shù) 迭代初數(shù)泊輸出得至|的解=-0-3472?6 請按任童鍵繼繕輸人初解=15輸入最大迭代次數(shù);5回 迭代枕數(shù):輸出得到的M = -0-347296 請按任盍鍵継繚-輸入初=輸入最大迭代初數(shù);別 迭代枕數(shù):解軸出得到的解= 1-87929 請無任意鍵輸入初i?=2-0驚入最夭迭祀枕數(shù)汚0 逹代枕數(shù)七輸岀得1的解= 1-87939 諳

5、轉(zhuǎn)任竟鍵繼續(xù).程序運(yùn)行結(jié)果討論和分析：對于第一種迭代格式，收斂區(qū)間-8.2 -0.4，在該收斂區(qū)間內(nèi)迭代收斂于-1.53209,只能求得方程的一個(gè)根；對于第二種迭代格式，收斂區(qū)間-1.51.8，在該收斂區(qū)間內(nèi)迭代收斂于-0.34730,同樣只能求得方程的一個(gè)根；對于第三種迭代格式，收斂區(qū)間-0.3 + g），在該收斂區(qū)間內(nèi)迭代收斂于1.87937,只能求得方程的一個(gè)根；由以上結(jié)果很容易發(fā)現(xiàn)，初值的選取對迭代斂散性有很大影響。以第一種迭代格式為例，當(dāng)初值大于等于-0.3時(shí)，迭代格式發(fā)散；當(dāng)初值小于等于-8.3時(shí)，迭代 格式也發(fā)散；只有初值在-0.3和-8.3之間時(shí)，迭代格式才收斂于一1.5320

6、9。其他 迭代格式也有這樣的性質(zhì)，即收斂于某個(gè)數(shù)值區(qū)間，超出這個(gè)區(qū)間迭代格式就是 發(fā)散的，這就是所謂迭代格式的收斂性。自主題：機(jī)械運(yùn)動的數(shù)值仿真背景：描述物理學(xué)里把物體位置的變化叫機(jī)械運(yùn)動。如我們所知，力的作用效果有：改變物體的運(yùn)動狀態(tài) 改變物體的形狀 改變物的運(yùn)動狀態(tài)大多會引起物 體的位置變化，引起機(jī)械運(yùn)動。改變物體的形狀而不改變它的運(yùn)動狀態(tài)就叫是非機(jī)械運(yùn)動中的一種。在工程實(shí)際與生產(chǎn)生活中我們常常需要對特定對象的機(jī)械 運(yùn)動進(jìn)行研究分析得出其機(jī)械運(yùn)動的規(guī)律，然后將該規(guī)律應(yīng)用于對人們有益的方面。實(shí)例：如下圖1 （左）所示，假設(shè)有一煙花火箭，其初始條件為零。將其放在地 方然后點(diǎn)火，該煙花火箭的初始

7、質(zhì)量為 m° 120g，其中粉末燃料占70g。經(jīng)過實(shí) 驗(yàn)得知，燃料的持續(xù)時(shí)間為tc 2.0s。燃料所產(chǎn)生的恒定推力為T 5.2N。這也說 明燃料的消耗率恒定?？諝猱a(chǎn)生的阻力和煙花火箭的速度的平方成正比：2422R kv,k 4.010 （Ns /m）。這里，要求選擇一種數(shù)值方法對其運(yùn)動過程進(jìn)行仿真4并且其截?cái)嗾`差為°（h）或者更高。要求計(jì)算出該煙花火箭的最高高度，同時(shí)計(jì) 算出從燃料消耗到該煙花火箭運(yùn)動到最高點(diǎn)的時(shí)間延遲。4數(shù)學(xué)原理及數(shù)學(xué)模型：該實(shí)際問題要求其截?cái)嗾`差要求大于或等于O（h ），這就使得較為簡單的歐拉法，中點(diǎn)法不適合本例。龍格-庫塔法以其優(yōu)異的數(shù)值特性成為解決本

8、問題的首選。圖1煙花火箭的機(jī)械運(yùn)動的數(shù)值仿真（左）和結(jié)構(gòu)力學(xué)問題的數(shù)值求解（右）很顯然，該問題屬于變質(zhì)量的運(yùn)動學(xué)問題， 在該運(yùn)動過程中，其前兩秒是在驅(qū)動 力和阻力的共同作用下加速上升的， 而后的時(shí)間內(nèi)，該煙花火箭是在空氣的阻力 下減速上升的，同時(shí)注意到空氣的阻力和速度的平方成正比。 為了對該運(yùn)動過程 進(jìn)行數(shù)值仿真，那么必須建立相應(yīng)的微分方程組。 分析該運(yùn)動過程可知，應(yīng)該將 該運(yùn)動過程分為兩部分：加速上升過程和減速上升過程。從而得到相應(yīng)的微分方 程組。加速上升過程：dv11000T1000kv2dt120 35tg120 35tdgvdt1Vi(O)0, h1(0)0減速上升過程:dh2V2圖2

9、煙花火箭的上升過程高度（左）和速度（右）的數(shù)值仿真dtdv21000kv2gdt 120 35 2h2（0） hi（2）V2（0） Vi（2）（2）式中，h為上升的高度，v為上升過程的速度，g為重力加速度。h（2）表示加速 上升過程的最終高度，v表示加速上升過程的最終速度。使用龍格-庫塔法求解如上的微分方程組。該系統(tǒng)的數(shù)值仿真結(jié)果如圖2所示。圖2 （左）是該煙花火箭的上升過程高度的數(shù)值仿真，圖2 （右）是其上升過程速度的數(shù)值仿真。同時(shí)亦可以得到煙花火箭上升的最大高度和問題中所需的時(shí)間 延遲：hmax198.462mt delay6.185（ S）從如上的分析和仿真可知，使用數(shù)值方法進(jìn)行機(jī)械運(yùn)動

10、的數(shù)值仿真，可以簡化本 身復(fù)雜變化的物理運(yùn)動過程。同時(shí)注意到，即使改變機(jī)械運(yùn)動的初始條件或者系 統(tǒng)的某些特征，其對應(yīng)的仿真只需要進(jìn)行簡單的調(diào)整。 而且其仿真精度也可以得 以預(yù)見。編程過程：加速上升過程fx_,y_：=5200/(120-35t)-0.4xA2/(120-35t)-10;gx_,y_：二x；x,y=0,0;h=0.1;t=0.1;xx=Table0,i,1,70;yy=Table0,i,1,70;tt=Table0,i,1,70;Doa=fx,y;xa=x+h (a+fx+h,y+h*a)/2; b=gx,y;ya=y+h (b+gx+h,y+h*b)/2;Prin tk,&qu

11、ot;",t,"",xa,"",ya;t,x,y,xxk,yyk,ttk=t+h,xa,ya,xa,ya,t+h,k,1,20£7B9101112131-4IS lfiIT1320O1 , 20.30,4 aS0.60.7a.eo甲9i.1.1le 21.31- 41-S1. fi1.71. S1_ *2.:a 4 3 =.2 . 00E7.060920.35352mml0®74414,esse2.1515713.64143.621837*7712631314510 474335,"4EE13.E工口740.260

12、317190344-843121.221349 47S-3ZSIOCS4,139230.32ss t eoss3 6 亠 05223.4石號$41.9684S 313T2.591255.130377.01113945fil .2 90*77CK100C65< 397 77EZ34Ttt0.6,1.1,1.6,2.1,2.6,3.1,3.6,4.1,4.6,5.1,5.6,6.1,6.6,7.1,7.6xx17.3174,37.9564,61.2759,85.0712,104.674,113.404,109.842,75.5189,8.73997,-59.2801,-96.8113,-11

13、1.485,-117.825,-121.341,-123.752yy 0.125,8.90869,28.0119,58.7749,101.435,153.898,210.725,265.771,303.655,308.15,278.635,230.354,174.737,115.949,55.4041ListPlotTabletti,xxi,i,1,20,PlotStylePoi ntSize0.02ListPlotTable tti,xxi,i,1,20,PlotStylePointSize0.02180L-60一-40一20*1 1 1 1 1 1 1i1aH1I1111iHi0.51.0

14、1.52.0ListPlotTabletti,xxi,i,1,20,PlotJoi nedTrue減速上升過程fx_,y_：=-0.4xA2心20-35*2)-10;gx_,y_：二x；x,y=85.3977,78.2347;h=0.1;t=2.1;Doa=fx,y;xa=x+h (a+fx+h,y+h*a)/2; b=gx,y;ya=y+h (b+gx+h,y+h*b)/2;Prin tk,"",t,"",xa,"",ya;t,x,y,xx20+k,yy20+k,tt20+k=t+h,xa,ya,xa,ya,t+h,k,1,502.

15、2 72. 13494.百M(fèi)L2.3 67-3S95L01+3062.4 62.751106.-&5E丐58.5956X14+93L54.S44412C.7952-751.433612.2552,843-3132131*4332 汽45.4421.眥，迪3,訛3.14C-31S4145.1023,236,014"145,133.3S5.B5Efi152.9453 + 433.E24215,5353,531.9062159,9233.6MoosesK3-H93*72S,362S146-1333.32-7161-62.3743 + 3£5+li36ii,esi4,23.

16、f55B174,174.122-15717-E.5334.22C,79L4le.7c2313.4433】麹.S4744IS.13994.5W,e?S2154.£1S123456189I口II1213U1$IE1"1 =空20212Z232425107-377189.327Z23293031325334箔3537IS3940鼻42434445箕47454350生首4$i.l14.44534.B13.250712.133511.0197.9217191,S79192,9365.2e .a 421.5193.9335.37-77089194-S7Z5*4.729351954655,333s.e4空適19孔跡S.73.«4SS?157-379&.£L3-.7435.9l.«32L6198.018乞(L 629£»4ISE.iefi5.2fi.36.45.7-l,37O5fi-2.37136-3*37627-4,35513-S*40016.42307-7-45556-B.4994419C.254196,222138-09197«S5£137.525197.09