第四講 bvp4c函數(shù)應(yīng)用

第四講bvp4c函數(shù)的應(yīng)用,bvp4c針對(duì)兩點(diǎn)邊界值問(wèn)題進(jìn)行數(shù)值求解，實(shí)現(xiàn)對(duì)以下n個(gè)一階常微分方程的數(shù)值解：x的取值范圍為axb區(qū)間內(nèi)，初始條件,bvp4c調(diào)用格式,FunctionName函數(shù)的接口,functiondydx=FunctionName(x,y,p1,p2,.)其中x是一標(biāo)量,y為yj列向量,p1,p2,.等為fj的參數(shù)，參數(shù)值已知，輸出dydx是列向量，對(duì)應(yīng)于fj列向量,BCFunction函數(shù)的創(chuàng)建,functionRes=BCFunction(ya,yb,p1,p2,.)ya是表示yj(a)的列向量,yb是表示yj(b)的列向量，即使不要求邊界條件，已知參數(shù)p1,p2等也必須出現(xiàn)在接口定義語(yǔ)句中，輸出Res為列向量。,調(diào)用bvp4c函數(shù)所涉及到的另外兩個(gè)函數(shù),變量solinit為一結(jié)構(gòu)？，可由函數(shù)bvpinit得到：solinit=bvpinit(x,y)，x向量為初始網(wǎng)格點(diǎn)的估計(jì)值，向量y為每個(gè)yj的估計(jì)值。向量x和y的長(zhǎng)度互不相關(guān)。bvp4c的輸出sol是一結(jié)構(gòu)，為特定數(shù)量點(diǎn)對(duì)應(yīng)的解。為使曲線變得更光滑，需要在中間插入一些點(diǎn)，使用：sxint=deval(sol,xint),xint為點(diǎn)向量，函數(shù)deval根據(jù)這些點(diǎn)向量求解。sol為函數(shù)bvp4c的輸出。,bvp4c函數(shù)應(yīng)用舉例,求解下面這個(gè)單變量二階常微分方程將上式重寫(xiě)為一個(gè)一階方程組：,相關(guān)函數(shù)的創(chuàng)建,下面創(chuàng)建函數(shù)和子函數(shù)，實(shí)現(xiàn)一階常微分方程的子函數(shù)為OdeBvp，記錄邊界條件的函數(shù)為OdeBC,初始估計(jì)值由bvpinit給出，在0和1之間選擇5個(gè)點(diǎn)，假定y1有一個(gè)解-0.05，y2有一個(gè)解0.1,設(shè)k=2,則程序?yàn)?functionbvpExamplek=2;solinit=bvpinit(linspace(0,1,5),-0.05,0.1)exmpsol=bvp4c(OdeBvp,OdeBC,solinit,k);x=linspace(0,1,50);y=deval(exmpsol,x);plot(x,y(1,:)functiondydx=OdeBvp(x,y,k)dydx=y(2);x-k*y(1);functionres=OdeBC(ya,yb,k)res=ya(1);yb(1);,求解結(jié)果顯示,solinit和exmpsol結(jié)構(gòu),solinit=x:00.250000000000000.500000000000000.750000000000001y:2x5doubleexmpsol=x:1x7doubley:2x7doubleyp:2x7doublesolver:bvp4c0,邊值問(wèn)題舉例工程梁,材料力學(xué)解法回顧,根據(jù)材料力學(xué)中關(guān)于軸力，剪力和彎矩符號(hào)的規(guī)定：軸力，拉力為正；剪力，繞隔離體順時(shí)針?lè)较蜣D(zhuǎn)動(dòng)者為正；彎矩，使梁的下側(cè)纖維受拉者為正。畫(huà)出該懸臂梁的剪力圖和彎矩圖。,求解過(guò)程,用積分法，可推出梁的斜率（轉(zhuǎn)角）和撓度表達(dá)式為,工程梁的微分控制方程推導(dǎo),工程梁的無(wú)量綱方程,如果定截面梁橫向位移為w(x)，長(zhǎng)為L(zhǎng)，彈性模量是E，橫截面的慣性矩是I，單位長(zhǎng)度載荷是p0q(x),則無(wú)量綱方程如下：,梁的微分控制方程轉(zhuǎn)化為一階常微分方程組,用bvp4c函數(shù)求解受均布載荷的懸臂梁,懸臂梁的bvp4c函數(shù)求解,functionCantileverBeamsolinit=bvpinit(linspace(0,1,10),0.5,0.5,0.5,0.5);beamsol=bvp4c(BeamODE,BeamCantileverBC,solinit,);eta=linspace(0,1,50);y=deval(beamsol,eta);plot(eta,y(1,:)functiondydx=BeamODE(x,y)dydx=y(2);y(3);y(4);1;functionbc=BeamCantileverBC(y0,y1)bc=y0(1);y0(2);y1(3);y1(4);,懸臂梁的bvp4c函數(shù)求解結(jié)果,課堂練習(xí)完成簡(jiǎn)支梁(一端施加一額外力矩Mr=0.8）的bvp4c函數(shù)程序,functionSimpleBeamMr=0.8;solinit=bvpinit(linspace(0,1,10),0.5,0.5,0.5,0.5);beamsol=bvp4c(BeamODE,SimpleBeamBC,solinit,Mr);eta=linspace(0,1,50);y=deval(beamsol,eta);plot(eta,y(1,:)functiondydx=BeamODE(x,y,Mr)dydx=y(2);y(3);y(4);1;functionbc=SimpleBeamBC(y0,y1,Mr)bc=y0(1);y0(3);y1(1);y1(3)-Mr;,簡(jiǎn)支梁(一端施加一額外力矩Mr=0.8）的bvp4c求解結(jié)果,思考：用ode45函數(shù)和fsolve函數(shù)能用來(lái)求解梁的問(wèn)題嗎？,對(duì)于梁而言，本質(zhì)是一個(gè)邊值問(wèn)題，梁的兩端各為兩個(gè)邊界條件，ode45只能求解初值問(wèn)題，但對(duì)于梁的任一端，可以先假設(shè)另兩個(gè)初值已知，這樣用ode45函數(shù)求可以求解控制微