算額（その1444）
福島県福島市堂殿 黒沼神社 明治11年(1878)
街角の数学 Street Wasan ～落書き帳「○△□」～ 186．僕は十一歳
http://streetwasan.web.fc2.com/math16.9.22.html
キーワード：円6個，円弧
#Julia, #SymPy, #算額, #和算
円弧（弓形）の中に，甲円 2 個，乙円 2 個，丙円 3 個を容れる。丙円の直径が 1 寸のとき，甲円の直径はいかほどか。
円弧の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (x1, y + r1)
乙円の半径と中心座標を r2, (0, R - r2)
丙円の半径と中心座標を r3, (0, y + r3), (x3, y + r3)
y = R - 2r2 - 2r3
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms R::positive, r1::positive, x1::positive, r2::positive,
r3::positive, x3::positive
y = R - 2r2 - 2r3
x3 = 2x1
eq1 = x1^2 + (y + r1)^2 - (R - r1)^2
eq2 = x3^2 + (y + r3)^2 - (R - r3)^2
eq3 = x1^2 + (2r3 + r2 - r1)^2 - (r1 + r2)^2
eq4 = x1^2 + (r1 - r3)^2 - (r1 + r3)^2
res = solve([eq1, eq2, eq3, eq4], (r1, r2, R, x1))[1];
# r1 甲円の半径
res[1] |> println
res[1](r3 => 1/2).evalf() |> println
r3*(sqrt(3) + 2)/2
0.933012701892219
# r2 乙円の半径
res[2] |> println
res[2](r3 => 1/2).evalf() |> println
2*sqrt(3)*r3/3
0.577350269189626
# R 円弧の半径
res[3] |> println
res[3](r3 => 1/2).evalf() |> println
4*r3*(3 + 2*sqrt(3))/3
4.30940107675850
# x1 甲円の中心の x 座標
res[4] |> println
res[4](r3 => 1/2).evalf() |> println
r3*(1 + sqrt(3))
1.36602540378444
丙円の直径が 1 寸のとき，甲円の直径は 1.86602540378444 寸である。
「答曰、甲円径一寸八分有奇」と，切り捨てたであるが，「術曰、置七分五厘開平方加一個乗丙径、得甲径合問」で「甲円直径 = (sqrt(0.75) + 1)*丙円直径」と正確である。
ちなみに，乙円の直径は 1.15470053837925 寸である。
---
和算家の取り組みと違い，連立方程式を解くという方法は，連立方程式はそのままで求めるものを変えるだけで良いということがある。
・ 甲円の半径を与えて，乙円，丙円，円弧の半径を求めるときは以下のようになる。
res2 = solve([eq1, eq2, eq3, eq4], (r2, r3, R, x1))[1];
res2[1] |> println
res2[2] |> println
res2[3] |> println
4*r1*(-3 + 2*sqrt(3))/3
2*r1*(2 - sqrt(3))
8*sqrt(3)*r1/3
・ 乙円の半径を与えて，甲円，丙円，円弧の半径を求めるときは以下のようになる。
res3 = solve([eq1, eq2, eq3, eq4], (r1, r3, R, x1))[1];
res3[1] |> factor |> println
res3[2] |> println
res3[3] |> factor |> println
r2*(3 + 2*sqrt(3))/4
sqrt(3)*r2/2
2*r2*(sqrt(3) + 2)
・ 円弧の半径を与えて，甲円，乙円，丙円の半径を求めるときは以下のようになる。
res4 = solve([eq1, eq2, eq3, eq4], (r1, r2, r3, x1))[1];
res4[1] |> println
res4[2] |> factor |> println
res4[3] |> factor |> println
sqrt(3)*R/8
-R*(-2 + sqrt(3))/2
R*(-3 + 2*sqrt(3))/4
・ なんなら，甲円の中心の x 座標を与えて，甲円，乙円，丙円，円弧の半径を求めるときは以下のようになる。
res5 = solve([eq1, eq2, eq3, eq4], (r1, r2, r3, R))[1];
res5[1] |> println
res5[2] |> factor |> println
res5[3] |> factor |> println
res5[4] |> factor |> println
x1*(1 + sqrt(3))/4
-x1*(-3 + sqrt(3))/3
x1*(-1 + sqrt(3))/2
2*x1*(sqrt(3) + 3)/3
function draw(r3, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, r1, x1, r2) = (4r3*(3 + 2√3)/3, r3*(√3 + 2)/2, r3*sqrt(2√3 + 4), 2√3r3/3)
@printf("丙円の直径が %g のとき，甲円の直径は %g である。\n", 2r3, 2r1)
@printf("r3 = %g; R = %g; r1 = %g; x1 = %g; r2 = %g\n", r3, R, r1, x1, r2)
y = R - 2r2 - 2r3
x = sqrt(R^2 - y^2)
θ = atand(y, x)
x3 = 2x1
plot()
circle(0, 0, R, beginangle=θ, endangle=180-θ)
circle2(x1, y + r1, r1, :blue)
circle(0, R - r2, r2, :green)
circle(0, y + r3, r3, :magenta)
circle2(x3, y + r3, r3, :magenta)
segment(-x, y, x, y, :orange)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
point(x1, y + r1, "甲円:r1\n(x1,y+r1)", :blue, :center, delta=-3delta)
point(0, R - r2, "乙円:r2,(0,R-r2)", :green, :left, delta=-3delta, deltax=-15delta)
point(0, y + r3, "丙円:r3,(0,y+r3)", :magenta, :right, delta=-3delta, deltax=12delta)
point(x3, y + r3, "丙円:r3,(x3,y+r3)", :magenta, :right, delta=-3delta, deltax=12delta)
point(0, R, "R", :red, :center, :bottom, delta=3delta)
point(0, y, "y", :red, :center, delta=-3delta)
end
end;
draw(1/2, true)
