ホールケーキを三等分する
ホールケーキを等分に切り分けるとき，普通は扇形に切り分ける。
今回は，ケーキを二本の平行な直線で分割して，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
