蒙特卡洛求定積分及龍哥庫(kù)塔解微分方程

《蒙特卡洛求定積分及龍哥庫(kù)塔解微分方程》由會(huì)員上傳分享，免費(fèi)在線閱讀，更多相關(guān)內(nèi)容在工程資料-天天文庫(kù)。

1、作業(yè)一、蒙特卡洛法求定積分:0=(cosx+2.0)dx整體思路，蒙特卡洛求定積分的主要思想就是通過(guò)在取值范圍內(nèi)大量隨機(jī)數(shù)的隨機(jī)選取對(duì)函數(shù)進(jìn)行求值，進(jìn)而除以所取次數(shù)并乘以其區(qū)間，即為所積值。Step1：FileEditTextGoCellToolsDebugDesktc?o回1.0+■1.1X必癬O1Efunction[y3=f(X)2尸cos(x)+2.0:1?3end:'Editor-Untitled*在執(zhí)行程序前,打開(kāi)matlab,執(zhí)行以下操作打開(kāi)File—>New—>Functioix并點(diǎn)擊Function,彈出Editor窗口，將其中內(nèi)容修改為function[y]=f(x

2、)y=cos(x)+2?0;end(如圖所示)。Step2：在Editor窗口中點(diǎn)擊File->SaveAs,保存至所建立的文件夾內(nèi)，保存名稱必須為f.m。Step3：冋到Matlab程序中，將CurrentFolder改為與剛剛Function函數(shù)定義的保存路徑一致的路徑°Step4：在CommandWindow里輸入以卜?源程序。源程序：?n=10A6;?y=0;%用來(lái)表示精度，表示需要執(zhí)行次數(shù)。%初始化y=0,因?yàn)閒(x)二cos(x)+2.0,下面出現(xiàn)y=y+f(x),故爾y用來(lái)統(tǒng)計(jì)計(jì)算出的所有f(x)值的和。?count=l;?whilecountv二nx=unifrnd(

3、0,4);%從1幵始計(jì)數(shù)，顯然要執(zhí)行n次。%當(dāng)執(zhí)行次數(shù)少于n次時(shí)，繼續(xù)執(zhí)行%在matlab中，unifrnd是在某個(gè)區(qū)間內(nèi)均勻選取實(shí)數(shù)(可為小數(shù)或整數(shù))，表示x取0到4的任意數(shù)。y=y+f(x);count=count+1;%求出到目前為止執(zhí)行出的所有f(x)值的和，并用y表示%count+l表示已執(zhí)行次數(shù)再加一，將來(lái)通過(guò)與n的比對(duì)，判斷是否執(zhí)行下一次endz=y/n*4%如果執(zhí)行次數(shù)已達(dá)n次，那么結(jié)束循環(huán)%用y除以執(zhí)行次數(shù)①那么取得平均值，并用它乘以區(qū)間長(zhǎng)度4,即可得到z的近似值z(mì)=7.2430由于matlab中不能使用&字符，故用z代替，即可得到o=7.2395o在執(zhí)行過(guò)程中，發(fā)

4、現(xiàn)每一次執(zhí)行結(jié)果都會(huì)不一樣，這是因?yàn)槭请S機(jī)選取，所以不同次數(shù)的計(jì)算結(jié)果會(huì)有偏差，但由于執(zhí)行次數(shù)很多，從概率的角度來(lái)講都是與真實(shí)值相近似的值。作業(yè)二、用四階龍格庫(kù)塔法解微分方程:fd2yNew—'Function,并點(diǎn)擊Function,彈出Editor窗口，將其中內(nèi)容修改為將二階函數(shù)化為一階過(guò)程中，定義如下：functiondy=f

5、unc(~,y)dy=zeros(2,1);密初始化，2行一列,即(；器)=C)dy(l)=y(2);%將上述方程組化簡(jiǎn)得來(lái)dy(2)=-4*y(2)-29*y(1);拓將上述方程組化簡(jiǎn)得來(lái)end%定義結(jié)束回到Matlab程序中，將CurrentFolder改為與剛剛Functional數(shù)定義的保存路徑一致的路徑。在CommandWindow里輸入以下源程序。>>[xzy]=ode45(*func1,[012]z[015])繪圖>>y=size(y)y=2292%數(shù)組y包含的元素?cái)?shù)>>reshape(y,458,1)%將含有229*2=458格元素的數(shù)組y化為458行1列的數(shù)組>>p

6、lot(x,y)%繪圖>>title(1數(shù)值圖T%圖名>>holdon%保持當(dāng)前圖形>>xlabel％力Ux坐標(biāo)>>ylabel(fy*)%加丫坐標(biāo)數(shù)值圖15-1%281012方法二2.最常用的R-K公式標(biāo)準(zhǔn)4階R-K公式為:Yi+i=并+石(K」+2K2+2K3+KJKi=f(Xi+yi)"fg雋,刃將Kt)K4=f(Xi+h,yi+hK3)在editor窗口中打開(kāi)一個(gè)新函數(shù)，對(duì)龍格庫(kù)塔函數(shù)定義：function[x,y]=lgkt4j(odefun,xspan,yO,h)%odefun為微分方程的函數(shù)描述，xspan為解區(qū)間[xO,xn],y0為初始條件，h為迭代步長(zhǎng)x=xspa

7、n(1):h:xspan(2);%定義區(qū)為從xspan(1)開(kāi)女臺(tái)至Ijxspan(2)且步長(zhǎng)為hk=l;%初始化k=ly(:zl)=yO(:);%初始化y(:,l)fori=l:length(x)-1為從第一次循壞開(kāi)始，共循壞(length(x)-1)次Kl=feval(ode￡un,x(k),y(:,k));%feva丄(函數(shù)名，xl,x2,...)K2=feval(odefunzx(k)+h/2,y(:rk)+h/2*Kl);K3=feval(ode