https://sangaku0418.hatenablog.com/entry/2023/08/10/210002

静岡県伊豆市江川邸享和2年(1802)

和算の館
http://www.wasan.jp/sizuoka/egawa.html

長野県長野市若穂清水寺観音堂寛政6年(1794)

中村信弥「改訂増補長野県の算額」県内の算額(P.45)
http://www.wasan.jp/zoho/zoho.html

キーワード：円4個，正三角形，矢
#Julia #SymPy #算額 #和算 #数学

正三角形内に 3 本の斜線を引き，分割された領域に 4 個の等円を内接させる。正三角形の一辺の長さがが与えられるとき，等円の直径はいかほどか。

正三角形の一辺の長さを \(a\)，等円の半径を \(r\) とする。

以下の連立方程式を解けばよいが，solve() では無理のようなので，nlsove() で数値解を求める。解析解は後半に記載。

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

using SymPy

@syms a::positive, r::positive, x::positive, y::positive,
x1::positive,y1::positive;

eq1 = distance(0, √Sym(3)a/2, a/2, 0, x, y) - r^2
eq2 = distance(0, √Sym(3)a/2, x1, y1, x, y) - r^2
eq3 = distance(0, √Sym(3)a/2, x1, y1, 0, a/2√Sym(3)) - r^2
eq4 = distance(x1, y1, a/2, 0, x, y) - r^2
eq5 = distance(x1, y1, a/2, 0, 0, a/2√Sym(3)) - r^2;

# res = solve([eq1, eq2, eq3, eq4, eq5])

function H(u)
(r, x, y, x1, y1) = u
return [
-r^2 + (-3*a/8 + 3*x/4 + sqrt(3)*y/4)^2 + (-sqrt(3)*a/8 + sqrt(3)*x/4 + y/4)^2, # eq1
-r^2 + (x - x1*(3*a^2 - 2*sqrt(3)*a*y - 2*sqrt(3)*a*y1 + 4*x*x1 + 4*y*y1)/(3*a^2 - 4*sqrt(3)*a*y1 + 4*x1^2 + 4*y1^2))^2 + (y - (3*a^2*y - 2*sqrt(3)*a*x*x1 + 2*sqrt(3)*a*x1^2 - 4*sqrt(3)*a*y*y1 + 4*x*x1*y1 + 4*y*y1^2)/(3*a^2 - 4*sqrt(3)*a*y1 + 4*x1^2 + 4*y1^2))^2, # eq2
4*a^2*x1^2*(3*a - 2*sqrt(3)*y1)^2/(9*(3*a^2 - 4*sqrt(3)*a*y1 + 4*x1^2 + 4*y1^2)^2) - r^2 + (sqrt(3)*a/6 - sqrt(3)*a*(3*a^2 - 4*sqrt(3)*a*y1 + 12*x1^2 + 4*y1^2)/(6*(3*a^2 - 4*sqrt(3)*a*y1 + 4*x1^2 + 4*y1^2)))^2, # eq3
-r^2 + (x - (a^2*x - 4*a*x*x1 - 2*a*y*y1 + 2*a*y1^2 + 4*x*x1^2 + 4*x1*y*y1)/(a^2 - 4*a*x1 + 4*x1^2 + 4*y1^2))^2 + (y - y1*(a^2 - 2*a*x - 2*a*x1 + 4*x*x1 + 4*y*y1)/(a^2 - 4*a*x1 + 4*x1^2 + 4*y1^2))^2, # eq4
a^2*y1^2*(-sqrt(3)*a + 2*sqrt(3)*x1 + 6*y1)^2/(9*(a^2 - 4*a*x1 + 4*x1^2 + 4*y1^2)^2) - r^2 + (-a*y1*(3*a - 6*x1 + 2*sqrt(3)*y1)/(3*(a^2 - 4*a*x1 + 4*x1^2 + 4*y1^2)) + sqrt(3)*a/6)^2, # eq5
]
end;

a = 10
iniv = [big"0.09", 0.27, 0.19, 0.17, 0.09] .* a
res = nls(H, ini=iniv);
println(res);
([1.1421256293696416, 2.5, 2.04587576018291, 1.5108901906526, 1.1735629018936662], true)

\(r = 1.14213;\ a = 10;\ x = 2.5;\ y = 2.04588;\ x_1 = 1.51089;\ y_1 = 1.17356\)
\(小円の直径 = 2r = 2.28425\)

解析解

正三角形の面積を 3 個の相似な三角形と中央の小さな正三角形の面積に分割する。
中央の正三角形の面積は \(3\sqrt{3}r^2\) である（注）。
下側の三角形に注目する。鈍角は 120°であるので，3 辺の和は \(2a + 2r/\sqrt{3}\) で，面積は \( \left (2a + 2r/\sqrt{3}\right)r/2\)である。
よって，
\(\sqrt{3}a^2/4 = 3\left (2a + 2r/\sqrt{3}\right )r/2 + 3\sqrt{3}r^2\)
を \(r\) について解けばよい。

注: 正三角形の内接円の半径が \(r\) のとき，正三角形の面積は \(3\sqrt{3}r^2\) である。

@syms a::positive, r::positive
eq = 3(2a + 2*r/√Sym(3))r/2 + 3√Sym(3)r^2 - √Sym(3)a^2/4 |> simplify
eq |> println

-sqrt(3)*a^2/4 + 3*a*r + 4*sqrt(3)*r^2

solve(eq, r)[1] |> println

a*(-sqrt(3) + sqrt(7))/8

直径は \(a \left (\sqrt{7} - \sqrt{3}\right )/4\) である。

\(a = 10\) のとき \(2.2842512587392836\) である。

(√7 - √3)/4 * 10

2.2842512587392836

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

function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
a = 10
(r, x, y, x1, y1) = res[1]
@printf("r = %g; a = %g; x = %g; y = %g; x1 = %g; y1 = %g\n", r, a, x, y, x1, y1)
@printf("小円の直径 = 2r = %g\n", 2r)
plot([a/2, 0, -a/2, a/2], [0, √3a/2, 0, 0], color=:black, linewidth=0.25)
circle(0, a/2√3, r)
for deg = 0:120:240
A = [cosd(deg) -sind(deg); sind(deg) cosd(deg)]
(x2, y2, p1, p2, p3, p4) = A * [x 0 x1; y-a/2√3 √3a/2-a/2√3 y1-a/2√3]
segment(p1, p2 + a/2√3, p3, p4 + a/2√3, :blue)
circle(x2, y2 + a/2√3, r)
end
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(0, a/2√3, " a/2√3", :red)
point(x, y, "(x,y)", :red)
point(x1, y1, "(x1,y1) ", :black, :right)
point(0, √3a/2, " √3a/2", :black)
point(a/2, 0, " a/2", :black, :left, :bottom)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
end
end;