46 岩手県一関市舞川字龍ヶ沢 観福寺観音堂 明治34年(1901)
安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：円11個，楕円
#Julia #SymPy #算額 #和算 #数学
小円を中心として楕円を 5 個配置する。楕円の内外に大円を 10個置く。小円の直径が与えられたときに大円の直径を求める術をのべよ。
楕円の長半径と短半径と中心座標を \(a,\ b,\ (0,\ r_2 + b)\)
大円の半径と中心座標を \(r_1,\ (0,\ r_2 + b)\)
小円の半径と中心座標を \(r_2,\ (0,\ 0)\)
楕円と右隣の楕円の接点座標を \( (x_0,\ y_0)\)
とおき，以下の連立方程式を解く。
なお，方程式中に \(\tan(54°) = \sqrt{2\sqrt{5}/5 + 1}\) の定数が出てくるが，これをそのまま用いると式が複雑になりすぎて SymPy では解がもとまらない。\(t = \tan(54°)\) として，\(t\) を使って方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
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\sqrt{5} + 5)\cdot(\sqrt{5} - 2)\cdot(\sqrt{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)
以下のアイコンをクリックして応援してください
