八六 埼玉県加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治13年（部分拡大図）（埼玉県加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg
キーワード：円7個，菱形，弦
#Julia #SymPy #算額 #和算 #数学
外円内に弦を境界として，上部に甲円 4 個，乙円 2 個を容れる。下部には菱形と，菱形に内接する楕円を容れる。
乙円の径を▢寸，楕円の短径を一寸五分としたとき，楕円の長径はいかほどか。
注：欠損した一文字は「答」から推測すると「一」であろう。なお，「答」にも「長径▢寸七分五厘有奇」と欠損文字があるが，こちらは「二」であろう。
外円の半径と中心座標を \(R,\ (0,\ 0)\)
甲円の半径と中心座標を \(r_1,\ (0,\ R - r_1),\ (x_1,\ R - 2r_1);\ x_1 = \sqrt{3}r_1\)
乙円の半径と中心座標を \(r_2,\ (x_2,\ y_a + r_2);\ y_a = R - 4r_1 = 2r_1 - R\)
楕円の長半径と短半径を \(a,\ b\)
楕円と菱形の接点の座標を \( (x_0,\ y_0)\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R::positive, r1::positive, x1::positive,
r2::positive, x2::positive,
a::positive, b::positive,
x0::positive, y0::positive
x1 = sqrt(Sym(3))r1
yb = r1 - R
xb = sqrt(R^2 - yb^2)
eq1 = x1^2 + (R - 2r1)^2 - (R - r1)^2
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq3 = (x1 - x2)^2 + (2r1 - r2)^2 - (r1 + r2)^2
eq4 = x0^2/a^2 + (y0 - r1 + R)^2/b^2 - 1
eq5 = -b^2*x0/(a^2*(y0 - r1 + R)) + r1/xb
eq6 = (y0 - yb)/(xb - x0) - r1/xb;
図の上半分は \(R,\ r_1,\ x_2\) がわかれば描ける。
res = solve([eq1, eq2, eq3], (R, r1, x2))
2-element Vector{Tuple{Sym, Sym, Sym}}:
(-2*sqrt(2)*r2 + 3*r2, r2*(3 - 2*sqrt(2))/3, -2*sqrt(6)*r2/3 + 2*sqrt(3)*r2/3)
(2*sqrt(2)*r2 + 3*r2, r2*(2*sqrt(2) + 3)/3, 2*sqrt(3)*r2/3 + 2*sqrt(6)*r2/3)
2 組の解が得られるが，2 番目のものが適解である。
res[2][1] |> factor |> println
res[2][2] |> factor |> println
res[2][3] |> factor |> println
r2*(2*sqrt(2) + 3)
r2*(2*sqrt(2) + 3)/3
2*r2*(sqrt(3) + sqrt(6))/3
r2 = 1/2
(2*sqrt(2)*r2 + 3*r2, r2*(2*sqrt(2) + 3)/3, 2*sqrt(3)*r2/3 + 2*sqrt(6)*r2/3)
(2.914213562373095, 0.9714045207910317, 1.3938468501173515)
下半分も描くには SymPy は力不足なので，数値解を求める。
function H(u)
(R, r1, x2, x0, y0, a) = u
return [
3*r1^2 + (R - 2*r1)^2 - (R - r1)^2, # eq1
x2^2 + (r1 - r2)^2 - (r1 + r2)^2, # eq2
-(r1 + r2)^2 + (2*r1 - r2)^2 + (sqrt(3)*r1 - x2)^2, # eq3
-1 + (R - r1 + y0)^2/b^2 + x0^2/a^2, # eq4
r1/sqrt(R^2 - (-R + r1)^2) - b^2*x0/(a^2*(R - r1 + y0)), # eq5
-r1/sqrt(R^2 - (-R + r1)^2) + (R - r1 + y0)/(-x0 + sqrt(R^2 - (-R + r1)^2)), # eq6
]
end;
r2 = 1//2
b = 15//20
iniv = BigFloat[2.91, 0.97, 1.39, 0.88, -1.36, 1.38]
res = nls(H, ini=iniv)
([2.914213562373095, 0.9714045207910317, 1.3938468501173518, 0.8773124760905445, -1.363750587600455, 1.3804469258418701], true)
楕円の長径は 2.7608938516837402 （長半径は 1.3804469258418701）である。
「答」では「長径は▢寸七分5厘有奇」となっている。
その他のパラメータは以下のとおりである。
\(R = 2.91421;\ r_1 = 0.971405;\ x_2 = 1.39385;\ y_a = -0.971405\)
\(x_0 = 0.877312;\ y_0 = -1.36375;\ a = 1.38045\)
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r2, b) = (1/2, 1.5/2)
(R, r1, x2, x0, y0, a) = res[1]
x1 = √3r1
ya = R - 4r1 # = 2r1 - R
xa = sqrt(R^2 - ya^2)
yb = r1 - R
xb = sqrt(R^2 - yb^2)
@printf("長径 = %g（長半径:a =%g）\n", 2a, a)
@printf("r2 = %g; b = %g\n", r2, b)
@printf("R = %g; r1 = %g; x2 = %g; ya = %g; x0 = %g; y0 = %g; a = %g\n",
R, r1, x2, ya, x0, y0, a)
plot()
circle(0, 0 , R, :blue)
circle(0, R - r1, r1)
circle(0, R - 3r1, r1)
circle(x1, R - 2r1, r1)
circle(-x1, R - 2r1, r1)
circle(0, r1 - R, r1, :gray90)
circle(x2, ya + r2, r2, :orange)
circle(-x2, ya + r2, r2, :orange)
segment(-xa, ya, xa, ya)
plot!([xb, 0, -xb, 0, xb], [yb, ya, yb, -R, yb], color=:green, lw=0.5)
ellipse(0, r1 - R, a, b)
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(0, R - r1, " 甲円:r1\n (0,R-r1)", :red, :left, :vcenter)
point(0, R - 3r1, "", :red)
point(x1, R - 2r1, "(x1,R-2r1)", :red, :center, delta=-delta/2)
point(x2, ya + r2, " 乙円:r2\n (x2,ya+r2)", :black, :left, :vcenter)
point(0, ya, " ya=R-4r1", :black, :center, :bottom, delta=delta)
segment(-xb, yb, xb, yb)
segment(0, ya, 0, -R)
point(0, r1 - R, "r1-R ", :black, :right, :bottom, delta=delta)
point(a, r1 - R, "(a,r1-R) ", :red, :right, :bottom, delta=delta)
point(0, r1 - R + b, " r1-R+b", :red, :left, delta=-delta)
end
end;
以下のアイコンをクリックして応援してください
