清華大學(xué)謝金星數(shù)學(xué)實(shí)驗(yàn)作業(yè)

1、實(shí)驗(yàn)二 常微分方程的數(shù)值解土木系 結(jié)23 李會(huì)平【實(shí)驗(yàn)?zāi)康摹?、 掌握用matlab解數(shù)值微分方程2、 了解龍格-庫(kù)塔方法的基本原理3、 用這些手段解決一些實(shí)際的問(wèn)題【實(shí)驗(yàn)內(nèi)容】4-5 核廢料問(wèn)題l 首先列出問(wèn)題的運(yùn)動(dòng)學(xué)方程，由牛頓第二定律，md2s/dt2=G-F-f其中m=G/g, f=kv,由于我們熟悉的單位是公制的，所以在定義函數(shù)的時(shí)候進(jìn)行了單位轉(zhuǎn)換。以下是函數(shù)的定義：這其實(shí)是一個(gè)關(guān)于s的二階常微分方程，需要定義兩個(gè)變量x(1),x(2)將其化為一階微分方程組，此處x(1),x(2)實(shí)際上分別代表速度和深度，相應(yīng)的微分方程如代碼中所示：dx=(G-F-k*x(1)/m;x(1)（如課堂

2、提醒，中間應(yīng)該是分號(hào)，這點(diǎn)容易出錯(cuò)。）l 接寫(xiě)來(lái)進(jìn)行m文件的命令編寫(xiě)，如下所示：（1）執(zhí)行pause之前的代碼，得到的結(jié)果如下圖所示：該圖代表了假設(shè)水無(wú)限深的情況下，物體速度的變化情況，容易看出速度是有一個(gè)上限的，這也符合直觀的感受，因?yàn)樗俣炔豢赡軣o(wú)限大，否則阻力無(wú)限大，物體將無(wú)法繼續(xù)運(yùn)動(dòng)。（2）紅線(xiàn)為題給的速度閾值，可見(jiàn)初步判斷速度是有可能超過(guò)該閾值的，需要進(jìn)一步判斷沉底的時(shí)候是否超過(guò)該速度，于是需要做出s和v的關(guān)系曲線(xiàn)，即x(2),x(1)的關(guān)系曲線(xiàn)，如下所示：從圖像中看出，s趨于時(shí)候，v也趨于極限值，題給的smax=300*0.3048m,在編程的時(shí)候，如果要繪制該水平直線(xiàn)，需要將其轉(zhuǎn)化

3、為數(shù)組才能進(jìn)行繪圖，否則會(huì)出現(xiàn)錯(cuò)誤。（3）但是從這幅圖像中，由于s非常大，紅線(xiàn)觸底，無(wú)法直接看出s,v的相交情況，需要進(jìn)一步限定坐標(biāo)軸范圍進(jìn)行細(xì)化，如下所示：從圖中明顯讀出，當(dāng)s達(dá)到smax時(shí)候，其相應(yīng)的vvmax,所以從題給的情況看，工程師們的說(shuō)法是更有道理的。在這個(gè)例題中，加深了對(duì)微分方程數(shù)值解的理解，同時(shí)在用matlab繪圖時(shí)也有了更多的心得。繪圖時(shí)候尤其要注意為數(shù)組和數(shù)組的對(duì)應(yīng)！另外在通過(guò)查詢(xún)相關(guān)知識(shí)知道，如果要繪制直接過(guò)（a,b）和（c,d）兩點(diǎn)的直線(xiàn)，可以直接用plot(a,c,b,d)命令快速實(shí)現(xiàn)，方便了一些作圖過(guò)程。4-6小船渡河問(wèn)題l 首先給出問(wèn)題的理論解法：建立小船運(yùn)行的數(shù)

4、學(xué)模型，用極坐標(biāo)，以B為原點(diǎn)進(jìn)行計(jì)算。取定BA為極軸。小船位置設(shè)為C，可以得到方程組： ddt=v1sin-v2ddt=v1cos/ 消去td=v1sin-v2v1cosd積分得到：ln=ln1cos12kln1+sin1-sin+Ck=v1v2可以解得：=Lcos(1+sin1-sin)-12kl 接下來(lái)用程序進(jìn)行數(shù)值解和理論解答的分析：（1）用matlab運(yùn)算時(shí)，采用普通坐標(biāo)更加方便，所以以A為原點(diǎn)建立直角坐標(biāo)系，容易得到 dxdt=v1-v2xx2+d-y2dydt=v2(d-y)x2+d-y2 據(jù)此定義函數(shù)m文件如下所示：（2）相應(yīng)的主程序?yàn)椋浩渲衟ause前的代碼表示進(jìn)行v1=1，v

5、2=2時(shí)候的數(shù)值運(yùn)算求解繪圖。pause后的代碼表示分別繪制出v1=0,0.5,1,1.5,2時(shí)候的運(yùn)動(dòng)軌跡并實(shí)現(xiàn)題目（2）中要求的數(shù)值解與理論解的對(duì)比情況，運(yùn)行的結(jié)果如下所示：v1=1m/s,v2=2m/s時(shí)候的軌跡圖：（2）然后是v1=0,0.5m/s,1m/s,1.5m/s,2m/s的理論曲線(xiàn)，并實(shí)現(xiàn)v1=1m/s時(shí)候的理論曲線(xiàn)和實(shí)際曲線(xiàn)的對(duì)比：該圖中有兩點(diǎn)需要注意的地方：1）標(biāo)注均在曲線(xiàn)的右方，而v1=0時(shí)候的曲線(xiàn)實(shí)際上和y軸重合，所以繪出的曲線(xiàn)沒(méi)有顯現(xiàn)出來(lái)；2）注意到v1=1m/s時(shí)候的曲線(xiàn)，實(shí)際上此處還有一條藍(lán)色曲線(xiàn)，圖上可以看到藍(lán)色曲線(xiàn)的痕跡，這是v1=1m/s時(shí)候的數(shù)值解曲線(xiàn)，可以看出二者的重合度是非常高的，幾乎完全重合，這說(shuō)明之前的數(shù)值解的結(jié)果是可信的。另外我們注意到1）當(dāng)v1=2m/s的時(shí)候，船已經(jīng)不可能正好到達(dá)河的正對(duì)岸了，從物理規(guī)律也可以直觀的理解，因?yàn)楫?dāng)v1v2時(shí)候，v2的任何一個(gè)分量都不可能大于v1，這樣一定會(huì)隨著水流向下游走一定的距離，而不可能完全到達(dá)對(duì)岸；2）當(dāng)水流速度逐漸增大的時(shí)候，路徑越來(lái)越凸向右側(cè)，即要到達(dá)正對(duì)岸會(huì)更加費(fèi)力氣，這也是符合生活規(guī)律的。【實(shí)驗(yàn)總結(jié)】本次實(shí)驗(yàn)總的來(lái)說(shuō)思路較為清晰，通過(guò)實(shí)驗(yàn)進(jìn)一步加深了對(duì)matlab的了解，鞏固了課堂所學(xué)的常微分方程的一些知識(shí)。在matlab的實(shí)現(xiàn)過(guò)程中，有很多需要注意的