以棱長為1的正四面體的各個頂點為球心作半徑為1的球,那麼四個球的公共部分的體積是多少?
是否能用高中知識來做,一定要定積分才能算嗎?
這個立體叫萊洛四面體。
我想了好久,發現想清楚它的形狀都很有挑戰性……這兒有一篇聲稱用純幾何方法求出萊洛四面體的表面積和體積的文章:
http://www.math.unl.edu/~bharbourne1/ST/sphericaltetrahedron.html用 PolyhedronData 獲得正四面體的頂點,作四個球,然後取其交集,
r=RegionIntersection@@(Ball/@N@PolyhedronData["Tetrahedron","VertexCoordinates"]);
可以數值地計算出區域的測度,
RegionMeasure[r]
(* 0.422158 *)
RegionPlot3D[r,PlotPoints-&>150,PlotStyle-&>Opacity[0.5]]
瀉藥,這題求解的正統方法可能是多重積分,但是關於積分邊界的確定貌似就很麻煩,我就嘗試一下數值方法了,權當拋磚引玉。
第一次看到這題的直覺是通過幾何概型來實現(Monte-Carlo Method),不妨假設這個正四面體的各頂點坐標分別為,不難發現前三個點恰好在平面內。
那麼我們以這四個頂點分別為球心,半徑為1做球,不難得出球內所有點的取值範圍:
.接下來做的事情很簡單,就是在這樣的一個長方體(由上面三個區間唯一確定)內隨機地生成大量的點,如個,要求這些點在各個軸上滿足均勻分布(Uniformly Distributed),我們統計其中恰好同時處於四個球內部的點的個數,那麼按照幾何概型的思想,我們所要求的目標區域的體積就是
代碼實現你們就湊合著看,MATLAB我實在不太熟悉:%%
p1 = [0 0 0]; p2 = [1 0 0]; p3 = [.5 sqrt(3)/6 sqrt(6)/3]; p4 = [.5 sqrt(3)/2 0];
n = 1000000;
x = -1 + 3*rand(n, 1);
y = -1 + (2+sqrt(3)/2)*rand(n, 1);
z = -1 + (2+sqrt(6)/3)*rand(n, 1);
points = [x y z];
d1 = bsxfun(@minus, points, p1);
d2 = bsxfun(@minus, points, p2);
d3 = bsxfun(@minus, points, p3);
d4 = bsxfun(@minus, points, p4);
is_in_sphere1 = sum(d1.^2, 2) &<= 1; is_in_sphere2 = sum(d2.^2, 2) &<= 1; is_in_sphere3 = sum(d3.^2, 2) &<= 1; is_in_sphere4 = sum(d4.^2, 2) &<= 1; total = sum(prod([is_in_sphere1, is_in_sphere2, is_in_sphere3, is_in_sphere4], 2)); volume = total/n*3*(2+sqrt(3)/2)*(2+sqrt(6)/3) volume = 0.4219
而根據 @王贇 Maigo給出的文獻中的公式:
令,計算結果為 0.4222,相對來說誤差不算太大。謝邀,但是我不會.感覺題主的圖形是個粽子.我之前生成過一個粽子的形狀.
vertices = D1:100 D2:100
u = from 0 to (PI*2) D1
v = from 0 to (PI*2) D2
x = cos(u)
y = cos(v)
z = cos(u + v)
給題主做個參考吧.
在包含這個粽子的長方體里隨機取點然後判斷是否屬於這個粽子23333用概率求。。
k=Integrate[ Boole[x^2 + y^2 + z^2 &< 2 (x - 1)^2 + (y - 1)^2 + z^2 &< 2 (x - 1)^2 + y^2 + (z - 1)^2 &< 2 x^2 + (y - 1)^2 + (z - 1)^2 &< 2], {x, -3, 3}, {y, -3, 3}, {z, -3, 3}] 答案就是k/(2sqrt(2))然而我電腦算不出,前面加個N倒是能求k=1.19404,本題答案0.422158
matlab代碼:
f = @(x,y,z)(x.^2 + y.^2 + (z-sqrt(3/8)).^2 &< 1) .* ((x+1/sqrt(12)).^2 + (y+1/2).^2 + (z+1/sqrt(24)).^2 &< 1) .* ((x+1/sqrt(12)).^2 + (y-1/2).^2 + (z+1/sqrt(24)).^2 &< 1) .* ((x-1/sqrt(3)).^2+y.^2+(z+1/sqrt(24)).^2&<1);
integral3(f, -1, 1, -1, 1, -1, 1,method,iterated)
但是為什麼算這麼慢?
@Falccm推薦閱讀:
TAG:數學 | 微積分 | MATLAB | WolframMathematica | 解析幾何 |