算額（その171） - 算額あれこれ

岐阜県大垣市外野釜笛釜笛八幡神社慶應元年(1865)

http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF
キーワード：円7個，半円，直角二等辺三角形
#Julia #SymPy #算額 #和算 #数学

第三問　半円内に直角二等辺三角形を作り，その隙間に青，赤の 6 個の円を容れる，半円の直径を知って，赤円の径を求めよ。

半円の径を \(R\) とする。後々 \(R\) は 8 の倍数にしておくほうが方程式の係数が簡単になことがわかる。

白円，青円，赤円の半径を \(r_0, r_1, r_2\) とし，赤円の中心座標を \( (x_2, y_2)\) とする。

\(r_0, r_1\) はすぐに計算できる。\(r_2, x_2, y_2\) が未知数であるが，solve() では以下の 4 個の方程式のどの 3 つをとっても，有限の時間内には解が見つからないようである。
そこでいつものように nlsolve で解を求める。

include("julia-source.txt"); # julia-source.txt ソース

using SymPy

@syms R::positive, r0::positive, r1::positive, r2::positive, x1::positive, x2::positive, y2::positive;

R = 16
r0 = R/(1 + sqrt(Sym(2)))
r1 = R*(2 - sqrt(Sym(2)))/4
x1 = (R - r1) / sqrt(Sym(2))
r2 = R // 8
eq1 = distance(R, 0, 0, R, x2, y2) - r2^2
eq2 = x2^2 + y2^2 - (R - r2)^2
eq3 = (x2 - x1)^2 + (y2 - x1)^2 - (r1 + r2)^2
eq4 = x2^2 + (y2 - r0) - (r0 + r2)^2;

println(eq1, ",")
println(eq2, ",")
println(eq3, ",")
println(eq4, ",")

2*(x2/2 + y2/2 - 8)^2 - 4,
x2^2 + y2^2 - 196,
(x2 - sqrt(2)*(4*sqrt(2) + 8)/2)^2 + (y2 - sqrt(2)*(4*sqrt(2) + 8)/2)^2 - (10 - 4*sqrt(2))^2,
x2^2 + y2 - (2 + 16/(1 + sqrt(2)))^2 - 16/(1 + sqrt(2)),

function H(u)
(r2, x2, y2) = u
return [
2*(x2/2 + y2/2 - 8)^2 - 4,
x2^2 + y2^2 - 196,
(x2 - sqrt(2)*(4*sqrt(2) + 8)/2)^2 + (y2 - sqrt(2)*(4*sqrt(2) + 8)/2)^2 - (10 - 4*sqrt(2))^2,
# x2^2 + y2 - (2 + 16/(1 + sqrt(2)))^2 - 16/(1 + sqrt(2)),
]
end;

iniv = [2.0, 6, 12]
res = nls(H, ini=iniv)
println(res)

([2.0, 6.352746103452378, 12.475681021293815], true)

経験則で \(r_2\) は \(R\) の 1/8 となる。
そこで \(r_2 = R // 8\) を条件に加えて eq1, eq2 の二元連立方程式を解く。

res = solve([eq2, eq1], (x2, y2))

2-element Vector{Tuple{Sym, Sym}}:
(-(sqrt(2) + sqrt(32 - 16*sqrt(2)) + 8)^3/34 - 30*sqrt(32 - 16*sqrt(2))/17 - 30*sqrt(2)/17 - 32/17 + 8*(sqrt(2) + sqrt(32 - 16*sqrt(2)) + 8)^2/17, sqrt(2) + sqrt(32 - 16*sqrt(2)) + 8)
(-(-sqrt(32 - 16*sqrt(2)) + sqrt(2) + 8)^3/34 - 30*sqrt(2)/17 - 32/17 + 30*sqrt(32 - 16*sqrt(2))/17 + 8*(-sqrt(32 - 16*sqrt(2)) + sqrt(2) + 8)^2/17, -sqrt(32 - 16*sqrt(2)) + sqrt(2) + 8)

\(x_2\) が特に複雑な式になり，これが solve() で解が求まらない原因と思われる。

赤円の半径は半円の径の 1/8 である。

描画関数プログラムのソースを見る

function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
plot()
R = 16
r0 = R / (1 + √2)
r1 = (R - R/√2)/2
x1 = (R -r1)/√2
r2 = R // 8
(x2, y2) = res[1]
(x22, y22) = res[2]
circle(0, 0, R, :yellow, fill=true, beginangle=0, endangle=180)
plot!([R, 0, -R, R], [0, R, 0, 0], linecolor=:palegreen2, linewidth=0.5, seriestype=:shape, fillcolor=:palegreen2)
circle(0, r0, r0, :white, fill=true)
circle(x2, y2, r2, :indianred1, fill=true)
circle(x22, y22, r2, :indianred1, fill=true)
circle(x1, x1, r1, :steelblue1, fill=true)
circle(-x2, y2, r2, :indianred1, fill=true)
circle(-x22, y22, r2, :indianred1, fill=true)
circle(-x1, x1, r1, :steelblue1, fill=true)
if more
point(x2, y2, "(x2,y2)", :black)
point(x1, x1, "( (R-r1)/√2,(R-r1)/√2)", :black, :center)
point(0, r0, " r0", :black)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
end
end;