算額（その1646）
三十二 一関市舞川 観福寺内地蔵堂後額 明治34年(1901)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
今有如図 03086
https://w.atwiki.jp/sangaku/pages/158.html
小円を中心として楕円を 5 個配置する。楕円の内外に大円を 10個置く。小円の直径が与えられたときに大円の直径を求める術をのべよ。
図の解釈としては「今有如図」を採った。
楕円の長半径と短半径と中心座標を a, b, (0, r2 + b)
大円の半径と中心座標を r1, (0, r2 + b)
小円の半径と中心座標を r2, (0, 0)
楕円と右隣の楕円の接点座標を (x0, y0)
とおき，以下の連立方程式を解く。
なお，方程式中に tan(54°) = sqrt(2*sqrt(5)/5 + 1) の定数が出てくるが，これをそのまま用いると式が複雑になりすぎて SymPy では解がもとまらない。t = tan(54°) として，t をつかって方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive,
x0::positive, y0::positive
@syms a, b, r1, r2, x0, y0, t
# t = tand(Sym(54))
a = r1
eq1 = x0^2/a^2 + (y0 - (r2 + b))^2/b^2 - 1
eq2 = -b^2*x0/(a^2*(y0 - (r2 + b))) - t
eq3 = (r2 + b)/2r1 - t;
eq4 = t - y0/x0;
res = solve([eq1, eq2, eq3, eq4], (r1, b, x0, y0))[2] # 2 of 2
(r2*(4*sqrt(3) + 7)/(t*(sqrt(3) + 2)), r2*(3 + 2*sqrt(3)), r2*(sqrt(3) + 2)/(2*t), r2*(sqrt(3) + 2)/2)
大円の半径は以下のようになる。
# r1
res[1] |> println
r2*(4*sqrt(3) + 7)/(t*(sqrt(3) + 2))
t に tan(54°) を代入すると以下のようになり，SymPy では直接的には簡約化できない。
res[1](t => tand(Sym(54)))
SymPy の力を借りながら手作業で簡約化して最終的に以下を得る。
大円の半径は，小円の半径の sqrt(2√5 + 5)*(√5 - 2)*(√3 + 2) = 2.71149362837554 倍である。
たとえば，小円の直径が 1 のとき，大円の直径は 2.71149362837554 である。
# r1
r1 = r2*sqrt(2√5 + 5)*(√5 - 2)*(√3 + 2)
r1 |> println
2.71149362837554*r2
r1(r2 => 1)
2.71149362837554
術は以下のようであるが，小円径 = 1 のとき，大円径 = 3.804226065180614 になってしまう。
術がおかしいのか，算額の図形の解釈が間違えているのかわからない。
一つの可能性としては山村の図のように 10 個の楕円の中心が同一円の円周上にあるような図形であるが，ちょっと立式が難しいか？
法 = sqrt(sqrt(0.8) + 1)
小円径 = 1
大円径 = (sqrt(5) + 3)*小円径/法
3.804226065180614
function draw(r2, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
t = tand(54)
(r1, b, x0, y0) = (r2*(4*sqrt(3) + 7)/(t*(sqrt(3) + 2)), r2*(3 + 2*sqrt(3)), r2*(sqrt(3) + 2)/(2*t), r2*(sqrt(3) + 2)/2)
r1 = r2*sqrt(2√5 + 5)*(√5 - 2)*(√3 + 2)
a = r1
@printf("r2 = %g; r1 = %g; a = %g; b = %g; x0 = %g; y0 = %g\n", r2, r1, a, b, x0, y0)
plot()
circle(0, 0, r2, :green)
l = (r2 + b)
l2 = l/cosd(36)
for i = 0:4
θ = 90 + 72i
ellipse(l*cosd(θ), l*sind(θ), a, b, color=:blue, φ=θ+90)
circle(l*cosd(θ), l*sind(θ), r1, :red)
circle(l2*cosd(θ + 36), l2*sind(θ + 36), r1, :red)
i == 4 && segment(0, 0, l2*cosd(θ + 36), l2*sind(θ + 36), :black)
end
x = l2.*cosd.(54:72:414)
y = l2.*sind.(54:72:414)
plot!(x, y, color=:gray80, lw=0.5)
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(x0, y0, " (x0,y0)", :blue, :left, :vcenter)
point(0, r2 + b, "大円:r1\n(0,r2+b)", :red, :center, delta=-delta/2)
end
end;
draw(1/2, true)
