求體積的方法

蒙特卡羅方法求冰淇淋錐的體積200820302021趙越一、問題敘述蒙特卡羅方法是一種通過隨機(jī)變量的統(tǒng)計(jì)試驗(yàn)求近似解的方法。在冰淇淋錐體和其外圍六面體中隨機(jī)產(chǎn)生N個(gè)坐標(biāo)點(diǎn)，統(tǒng)計(jì)其中落在錐體中的點(diǎn)和總的點(diǎn)數(shù)的比例就可以近似地求出錐體的體積，在該實(shí)驗(yàn)中求出空間中x(-l,l)，y(-l.l)，z(0.2)區(qū)間中的冰淇淋錐體的體積大小。二、問題分析x2+x2+y2&z<1%%設(shè)置函數(shù)名%設(shè)置循環(huán)次數(shù)最多為10%設(shè)置采樣點(diǎn)的數(shù)量N為10000%循環(huán)L次%隨機(jī)取值函數(shù)N行3列%設(shè)置x在(-1,1)范圍取值%設(shè)置y在(-1,1)范圍取值%設(shè)置z在(0,2)范圍取值%設(shè)置錐體的函數(shù)表達(dá)%設(shè)置錐體的函數(shù)表達(dá)%判斷落在錐體內(nèi)的點(diǎn)的數(shù)量%計(jì)算出錐體的體積大小%返回體積大小和計(jì)算誤差冰淇淋錐含于體積等于8的立方體中，G={(x,y,z)I-1<x<1,-1<y<1,0<z<2}由于rand產(chǎn)生0到1之間的隨機(jī)數(shù)，所以：x=2*rand-1產(chǎn)生-1到1之間的隨機(jī)數(shù)；y=2*rand-1產(chǎn)生-1到1之間的隨機(jī)數(shù)；z=2*rand;產(chǎn)生0到2之間的隨機(jī)數(shù)；N個(gè)點(diǎn)均勻分布于六面體中，錐體中占有m個(gè)，則錐體與六面體體積之比近似為m:N多次反復(fù)試驗(yàn)，可以觀察出其誤差情況三、實(shí)驗(yàn)程序及注釋functiondata=icecream(L)ifnargin==0,L=10;endN=10000;fork=1:LP=rand(N,3);x=2*P(:,1)-1;y=2*P(:,2)-1;z=2*P(:,3);R2=x.A2+y.A2;R=sqrt(R2);II=find(z>=R&z<=1+sqrt(1-R2));m=length(II);q(k)=8*m/N;enddata=[q;q-pi];四、實(shí)驗(yàn)數(shù)據(jù)結(jié)果及分析程序運(yùn)行后，分別列出了取值N=10000的時(shí)候進(jìn)行10次同樣的蒙特卡羅方法計(jì)算出來的錐體體積大小，可以看出小數(shù)點(diǎn)后第2位的數(shù)值已經(jīng)很不穩(wěn)定，因此只能保證兩位數(shù)字的有效數(shù)位，因此修改N的大小，進(jìn)行相同的實(shí)驗(yàn)，當(dāng)N分別取10人3,10人4,10人5,10人6時(shí)可以看出采樣點(diǎn)的多少對(duì)數(shù)據(jù)精度的影響。10八310八410八510飛13.0883.12243.12563.147923.243.13683.16783.144432.963.1683.14933.149443.0723.17683.15053.149753.0243.18483.12263.141763.1043.06323.13193.139273.3043.11923.15453.140482.9683.19843.15893.141593.123.11843.14023.1393103.1123.183.16143.1376在N為10A3時(shí)，數(shù)據(jù)精度很差，為10A4時(shí)只能保證不到2位的有效數(shù)位，隨著N的增加計(jì)算精度也隨著增加，在NT時(shí)，數(shù)據(jù)比較接近3.14。五、 實(shí)驗(yàn)結(jié)論蒙特卡羅方法近似計(jì)算幾何體體積，其隨著采樣點(diǎn)的數(shù)量提高，數(shù)據(jù)精度也逐漸提高，采樣點(diǎn)的增加對(duì)計(jì)算量的影響成正比例關(guān)系增加，而不是成幾何級(jí)數(shù)增加，因此使對(duì)圖形體積的快速計(jì)算成為可能，但數(shù)據(jù)精度較低。六、 注記蒙特卡羅方法可以用于快速計(jì)算函數(shù)曲線所圍圖形面積，即計(jì)算定積分，其優(yōu)點(diǎn)是計(jì)算量小，采樣點(diǎn)的提高對(duì)計(jì)算量提高不快，適合快速計(jì)算，計(jì)算面積的時(shí)候和計(jì)算體積的時(shí)候維度的增加對(duì)計(jì)算量的增加也是線性的，不像數(shù)值方法計(jì)算量提高很快，因此蒙特卡羅方法還適合多維度問題的計(jì)算。對(duì)拋射曲線的改進(jìn)實(shí)驗(yàn)200820302021趙越一，問題敘述在研究榴彈炮等類似的拋射曲線發(fā)射角度和射程的關(guān)系問題中，建立的物理模型為炮彈發(fā)射初速不變，炮彈發(fā)射后僅在重力的作用下作慣性飛行，飛行的軌跡為拋物線，但在實(shí)際問題中，用該模型計(jì)算出的射程會(huì)與實(shí)際的射程發(fā)生較大偏差，理論射程會(huì)大于實(shí)際射程，這主要是由于沒有考慮到炮彈在空氣中受到空氣的助力，速度不斷減小的原因。二，問題分析x=vcosaxt原拋物曲線公式為：o原拋物曲線公式為：<1y=vsinaxt一—gt2〔02考慮到炮彈在空氣中飛行受到空氣助力，設(shè)空氣阻力始終與飛行軌跡的切線方向相反，dV且阻力大小與速度呈正比例關(guān)系，側(cè)有：m叱=-ZV，其中Z為阻力系數(shù)，再根據(jù)初值,

dt當(dāng)t=0時(shí)，V=V0，可得到V=Vb*exp(-Zt),用V代替原拋物線公式中的V0得到：x=ve-ztcosaxt0<1y=ve-ztsinaxt一—gt2〔02%定義求解拋物線零點(diǎn)的函數(shù)先假設(shè)Z=0.01三，實(shí)驗(yàn)程序及注釋%定義求解拋物線零點(diǎn)的函數(shù)v0=515;alpha=45*pi/180;T=fzero('funx',50);v0=515;alpha=45*pi/180;T=fzero('funx',50);forn=1:1:17t=(n-1)*T/16;x(n)=v0*exp(-0.01*t)*t*cos(alpha);y(n)=v0*exp(-0.01*t)*t*sin(alpha)-end%定義初速度和發(fā)射角%求出并定義T為炮彈飛行的時(shí)間%計(jì)算17個(gè)點(diǎn)%定義采樣時(shí)間間隔%定義各點(diǎn)橫坐標(biāo):*t.A2; %定義各點(diǎn)縱坐標(biāo)plot(x,y,x,y,'r*')Xmax=x(17)四，實(shí)驗(yàn)數(shù)據(jù)結(jié)果及分析下圖為未經(jīng)過修正的原函數(shù)曲線

Xmax=1.0652e+004。已經(jīng)比較接近于真實(shí)值11800了。五，實(shí)驗(yàn)結(jié)論實(shí)驗(yàn)描述了通過對(duì)拋物曲線速度函數(shù)的修正，使得拋射曲線更加接近真實(shí)情況的方法，并通過程序化出了較理想的曲線圖，但為了逼近真實(shí)數(shù)據(jù)，還需要對(duì)阻力系數(shù)進(jìn)行精確的求解。六，注記下面嘗試求出助力系數(shù)的較精確的值。11800=v e-ztcosaxt0<1ve-ztsinaxt=—gt得出：一gt2得出：一gt2=11800tana=11800,t=49.07Z=0.0085重新畫曲線得出一functionfx=funx(x)fx=exp(-0.0085*x)-sqrt(2)*4.9*x/515;I0 2

v0=515;alpha=45*pi/180;T=fzero('funx',50);forn=1:1:17t=(n-1)*T/16;x(n)=v0*exp(-0.0085*t)*