算額（その1349）
長野県長野市戸隠 戸隠山中院権現堂 安永5年(1776)
『算法天元算題花 初遍上』
http://www.wasan.jp/terasima/terasima3.pdf
キーワード：球，分割，積分
#Julia, #SymPy, #算額, #和算
直径が 10 寸の球がある。図のように，道幅（隙間）が 1 寸になるように平行平面で切り，残りの体積（緑の部分）を 4 等分したい。このとき，「矢」を求める方程式を求めよ。
球の半径を r，矢を h，道幅を a とおき，以下の方程式を得る。
include("julia-source.txt");
using SymPy
@syms a, r, h, x
eq1 = integrate(PI*(r^2 - x^2), (x, r - h, r)) - integrate(PI*(r^2 - x^2), (x, a/2, r - a - h));
r = 10/2, a = 1 を代入し，π を消去して，方程式 -16*h^3 + 216*h^2 + 216*h - 1589 = 0 を得る。
eq11 = eq1(r => 10//2, a => 1) |> simplify |> x -> x*(24/PI)
eq11 |> println
-16*h^3 + 216*h^2 + 216*h - 1589
数値解は要求されていないが，solve() により数値解を求める。
res = solve(eq11); # 解は 3 個ある
(real(res[1].evalf()), real(res[2].evalf()), real(res[3].evalf())) |> println
(2.44853772716366, -2.90597254121110, 13.9574348140474)
pyplot(size=(500, 300), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(eq11, xlims=(−4, 15), xlabel="h", ylabel="体積の差")
hline!([0])
矢 = 2.44853772716366 が適解である。
# 検算
h = 2.44853772716366
r = 10/2
a = 1
integrate(PI*(r^2 - x^2), (x, r - h, r)) |> N |> println
integrate(PI*(r^2 - x^2), (x, a/2, r - a - h)) |> N |> println
78.80187353219384
78.80187353219385
