ホールケーキを三等分する

ホールケーキを三等分する

ホールケーキを等分に切り分けるとき，普通は扇形に切り分ける。
今回は，ケーキを二本の平行な直線で分割して，3分割されたケーキの体積がが等しくなるようにする方法を求める。
ケーキの大きさには無関係なので，単位円を二本の平行な直線で分割して，面積を等しくするように考える。

1.　結論

結論から言えば，円の中心から上下（または左右）へ，約 0.265 離れた平行線を引く。

include("julia-source.txt")

using Plots

function draw()
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plt = plot()
   θ = 0.26813348747179905/pi*180
   x, y = cosd(θ), sind(θ)
   @printf("(x, y) = (%g, %g)\n", x, y)
   circle(0, 0, 1, :red)
   segment(-x, y, x, y, :blue)
   segment(-x, -y, x, -y, :blue)
   point(x, y, " (x, y)")
   hline!([0], color=:gray, lw=0.5)
   vline!([0], color=:gray, lw=0.5)
end

   (x, y) = (0.964267, 0.264932)

2.　連立方程式を解く

平行線と円の交点座標を (x, y) として，以下の連立方程式を解く。

x^2 +  y^2 = 1
x*y/2 + atand(y/x)/360*pi = pi/12

SymPy で数式解を求めるには

@syms x, y
eq1 = x^2 +  y^2 - 1
eq2 = x*y/2 + atand(y/x)/360*PI - PI/12
solve([eq1, eq2], (x, y))

とすればよいが，以下のエラーが発生する。
NotImplementedError('could not solve -180*y*sqrt(-(y - 1)*(y + 1)) - pi*atan(pi*y/(180*sqrt(-(y - 1)*(y + 1)))) - 30*pi')

そこで，NLsolve により数値解を求める。

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (x, y) = u
   return [
       x^2 +  y^2 - 1
       x*y/2 + atand(y/x)/360*pi - pi/12
   ]
end;

iniv = BigFloat[0.95, 0.25]
res = nls(H, ini=iniv)
println("(x, y) = ($(Float64(res[1][1])), $(Float64(res[1][2])))")

   (x, y) = (0.9642670742838972, 0.26493208460277684)

(x, y) = (0.9642670742838972, 0.26493208460277684) となる。

3.　シミュレーションで検証

この解が適切かどうか，シミュレーションにより検証する。

第 1 象限の一辺の長さ 1 の正方形内に一様乱数点 (x, y) を n 組発生させ，(x, y) が四分円内にあり，かつ，y > y0 の点の割合が pi/3 になることを確認する。

(x0, y0) = (0.9642670742838972, 0.26493208460277684)
n = 100000000
x = rand(n)
y = rand(n)
res = x.^2 .+ y.^2 .< 1 .&& y .> y0
println(2*sum(res)/n, ",  ", pi/3)

   1.04723984,  1.0471975511965976

4.　定積分で検証

上部の面積を定積分により求める。

using SymPy
@syms x, y
println(integrate(sqrt(1 - x^2) - y0, (x, -x0, x0)))

   1.04719755119660

中央部の面積を定積分により求める。

println(4(x0*y0 + integrate(sqrt(1 - x^2), (x, x0, 1))).evalf())

   1.04719755119660