三五 埼玉県大宮市中釘 秋葉神社 天保11年(1840)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
埼玉県さいたま市西区中釘 秋葉神社 天保11年(1840)
山口正義：やまぶき，第20号
https://yamabukiwasan.sakura.ne.jp/ymbk20.pdf
キーワード：円5個，三角形
#Julia #SymPy #算額 #和算 #数学
三角形の中に甲円 3 個，乙円と丙円を 1 個ずつ容れる。甲円の直径が 20 寸のとき，乙円の直径はいかほどか。
解析解を求めるようにした
三角形の頂点を \( (0, 0), (x, y), (x_2, 0)\)
甲円の半径と中心座標を \(r_1, (x_1, r_1), (x_1 + 2r_1, r_1), (x_1 + 4r_1, r_1)\)
乙円の半径と中心座標を \(r_2, (x_1 + 3r_1, y_2)\)
丙円の半径と中心座標を \(r_3, (x_1 + r_1, y_3)\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a, b, c, d,
r1::positive, x1::positive,
r2::positive, y2::positive,
r3::positive, y3::positive,
x::positive, y::positive, x2::positive;
真ん中の甲円と乙円の中心の垂直距離を \(a\) とする。
甲円と乙円が外接することから，\(a^2 = (r_1 + r_2)^2 - r_1^2 = r_2(2r_1 + r_2)\) である。
eq1 = a^2 + r1^2 - (r1 + r2)^2
ans_a = solve(eq1, a)[2] # 2 of 2
ans_a |> println
sqrt(r2)*sqrt(2*r1 + r2)
乙円と丙円の中心の垂直距離を \(b\) とする。
乙円と丙円が外接することから \(b^2 = (r_2 + r_3)^2 - (2r_1)^2 = -4r_1^2 + r_2^2 + 2r_2 r_3 + r_3^2\) である。
eq2 = b^2 + (2r1)^2 - (r2 + r3)^2
ans_b = solve(eq2, b)[2] # 2 of 2
ans_b |> println
sqrt(-4*r1^2 + r2^2 + 2*r2*r3 + r3^2)
真ん中の甲円と丙円の中心の垂直距離を \(c\) とする。
甲円と丙円が外接することから，\(c^2 = (r_1 + r_3)^2 - r_1^2 = r_3(2r_1 + r_3)\) である。
eq3 = c^2 + r1^2 - (r1 + r3)^2
ans_c = solve(eq3, c)[2] # 2 of 2
ans_c |> println
sqrt(r3)*sqrt(2*r1 + r3)
\(a = b + c\) なのでこれを方程式としてもよいが，SymPy の能力的には両辺を二乗して \(a^2 = (b + c)^2 = b^2 + 2b c + c^2\) を方程式とするほうがよい。
eq123 = ans_a^2 - (ans_b + ans_c)^2 |> expand
eq123 |> println
4*r1^2 + 2*r1*r2 - 2*r1*r3 - 2*r2*r3 - 2*sqrt(r3)*sqrt(2*r1 + r3)*sqrt(-4*r1^2 + r2^2 + 2*r2*r3 + r3^2) - 2*r3^2
次に，左の甲円と乙円の中心間の距離を 2 通りの方法で記述する。
左の甲円と乙円の中心を斜辺 \(d_1\) とする直角三角形で，\(d_1^2 = a^2 + (3r)^2 = r_2(2r_1 + r_2) + 9r_1^2\) である。
@syms d1, d2
eq4 = d1^2 - (ans_a^2 + (3r1)^2)
ans_d1 = solve(eq4, d1)[2] # 2 of 2
ans_d1 |> println
sqrt(9*r1^2 + 2*r1*r2 + r2^2)
乙円と甲円の半径の差 \( (r_2 - r_1)\) と，左の甲円と乙円が三角形の斜辺の接点間の距離 \(2\sqrt{r_1 r_3} + 2\sqrt{r_2 r_3}\) を直角を挟む 2 辺とする直角三角形において，\(d_2^2 = (2\sqrt{r_1 r_3} + 2\sqrt{r_2 r_3})^2 + (r_2 - r_1)^2 = 4r_3(\sqrt{r_1} + \sqrt{r_2})^2 + (r_1 - r_2)^2\) である。
eq5 = d2 - ( (r2 - r1)^2 + (2sqrt(r2*r3) + 2sqrt(r1*r3))^2)
ans_d2 = solve(eq5, d2)[1]
ans_d2 |> println
4*r3*(sqrt(r1) + sqrt(r2))^2 + (r1 - r2)^2
\(d_1 = d_2\) を 2 つ目の方程式とする。
eq45 = (ans_d1)^2 - ans_d2 |> expand |> simplify
eq45 |> println
-8*sqrt(r1)*sqrt(r2)*r3 + 8*r1^2 + 4*r1*r2 - 4*r1*r3 - 4*r2*r3
eq123, eq45 の 2 つの方程式を解き \(r_2, r_3\) を求める。
res = solve([eq123, eq45], (r2, r3))[1]
(2*r1*(sqrt(7) + 4)/9, 2*r1*(sqrt(7) + 13)/(2*sqrt(7) + 17 + 6*sqrt(2)*sqrt(sqrt(7) + 4)))
# r2
ans_r2 = res[1]
ans_r2 |> println
ans_r2(r1 => 10).evalf() |> println
2*r1*(sqrt(7) + 4)/9
14.7683362468102
乙円の直径は，甲円の直径の \( (2\sqrt{7} + 8)/9\) 倍である。
甲円の直径が 20 寸のとき，乙円の直径は 29.5366724936204 寸である。
# r3
@syms d
ans_r3 = apart(res[2], d) |> sympy.sqrtdenest |> simplify |> factor
ans_r3 |> println
ans_r3(r1 => 10).evalf() |> println
-2*r1*(-3 + sqrt(7))
7.08497377870819
丙円の直径は，甲円の直径の \( (6 - 2\sqrt{7})\) 倍である。
甲円の直径が 20 寸のとき，乙円の直径は 14.169947557416371 寸である。
乙円，丙円の中心座標は以下のようになる。
# y2
eq3 = r1^2 + (y2 - r1)^2 - (r1 + r2)^2
ans_y2 = solve(eq3, y2)[2] # 2 of 2
ans_y2 |> println
ans_y2(r2 => ans_r2)(r1 => 10).evalf() |> println
r1 + sqrt(r2)*sqrt(2*r1 + r2)
32.6598870349137
左の甲円の中心座標が \( (x_1, r_1)\) のとき，乙円の中心座標は \( (x_1 + 3r_1, r_1 + \sqrt{r_2}\sqrt{2r_1 + r_2})\) である。
# y3
eq4 = r1^2 + (y3 - r1)^2 - (r1 + r3)^2
ans_y3 = solve(eq4, y3)[2] # 2 of 2
ans_y3 |> println
ans_y3(r3 => ans_r3)(r1 => 10).evalf() |> println
r1 + sqrt(r3)*sqrt(2*r1 + r3)
23.8526650511426
左の甲円の中心座標が \( (x_1, r_1)\) のとき，丙円の中心座標は \( (x_1 + r_1, r_1 + \sqrt{r_3}\sqrt{2r_1 + r_3})\) である。
三角形の頂点の座標を決定するのは複雑である。
別途数値計算を行い，甲円の直径が 20 寸のとき以下を得る。
x1 = 24.534049509916088
x = 63.120699623258076
y = 61.70734973660003
x2 = 77.04124269850618
function driver(r1, r2, r3, y2, y3)
function H(u)
(x1, x, y, x2) = u
return [
dist(0, 0, x, y, x1, r1) - r1^2,
dist(0, 0, x, y, x1 + r1, y3) - r3^2,
dist(x, y, x2, 0, x1 + 3r1, y2) - r2^2,
dist(x, y, x2, 0, x1 + 4r1, r1) - r1^2,
]
end;
iniv = BigFloat[24, 62.5, 62.5, 78] .* r1/10
res = nls(H, ini=iniv)
end;
描画関数プログラムのソースを見る
function draw(r1, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 = 2r1*(√7 + 4)/9
r3 = 2r1*(3 - √7)
y2 = r1 + sqrt(r2)*sqrt(2*r1 + r2)
y3 = r1 + sqrt(r3)*sqrt(2*r1 + r3)
res = driver(r1, r2, r3, y2, y3)
(x1, x, y, x2) = res[1]
@printf("甲円の直径 = %g; 乙円の直径 = %g; 丙円の直径 = %g\n", 2r1, 2r2, 2r3)
@printf("x1 = %g; r2 = %g; y2 = %g; r3 = %g; y3 = %g; x = %g; y = %g; x2 = %g\n", x1, r2, y2, r3, y3, x, y, x2)
plot([0, x2, x, 0], [0, 0, y, 0], color=:black, lw=0.5)
circle(x1, r1, r1, :green)
circle(x1 + 2r1, r1, r1, :green)
circle(x1 + 4r1, r1, r1, :green)
circle(x1 + r1, y3, r3, :blue)
circle(x1 + 3r1, y2, r2, :red)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
plot!(showaxis=false, xlims=(-5delta, x2+5delta))
point(x, y, "(x,y)", :black, :right, :bottom, delta=delta)
point(x1, r1, " 甲円:r1\n (x1,r1)", :green, :center, delta=-delta)
point(x1 + 2r1, r1, "甲円:r1\n(x1+2r1,r1)", :green, :center, delta=-delta)
point(x1 + 4r1, r1, "甲円:r1\n(x1+4r1,r1)", :green, :center, delta=-delta)
point(x1 + 3r1, y2, "乙円:r2,(x1+3r1,y2)", :red, :center, :bottom, delta=delta)
point(x1 + r1, y3, "丙円:r3,(x1+r1,y3) ", :blue, :right, :vcenter)
point(x2, 0, "(x2,0)", :black, :center, delta=-delta)
point(x1 + 3r1, r1)
point(0, 0, "(0,0)", :black, :center, delta=-delta)
dimension_line(x1, r1, x1 + 3r1, y2, "d", deltax=-2.5delta, length=0.2r1)
dimension_line(x1 + 3r1, r1, x1 + 3r1, y2, "a", deltax=1.5delta, length=0.2r1)
end
end;
draw(20/2, true)
以下のアイコンをクリックして応援してください
