求體積的方法

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

dt當t=0時，V=V0，可得到V=Vb*exp(-Zt),用V代替原拋物線公式中的V0得到：x=ve-ztcosaxt0<1y=ve-ztsinaxt一—gt2〔02%定義求解拋物線零點的函數(shù)先假設(shè)Z=0.01三，實驗程序及注釋%定義求解拋物線零點的函數(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為炮彈飛行的時間%計算17個點%定義采樣時間間隔%定義各點橫坐標:*t.A2; %定義各點縱坐標plot(x,y,x,y,'r*')Xmax=x(17)四，實驗數(shù)據(jù)結(jié)果及分析下圖為未經(jīng)過修正的原函數(shù)曲線

Xmax=1.0652e+004。已經(jīng)比較接近于真實值11800了。五，實驗結(jié)論實驗描述了通過對拋物曲線速度函數(shù)的修正，使得拋射曲線更加接近真實情況的方法，并通過程序化出了較理想的曲線圖，但為了逼近真實數(shù)據(jù)，還需要對阻力系數(shù)進行精確的求解。六，注記下面嘗試求出助力系數(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)*