算額（その1181）
(7) 京都府 妙見堂
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円7個，外円，デカルトの円定理
大円の中に天円，地円，人円，甲円，乙円，丙円が入っている。天円，地円，人円の直径がそれぞれ 14 寸，7 寸，3 寸のとき，乙円の直径はいかほどか。
大円，天円，地円，人円，甲円，乙円，丙円の半径をそれぞれ r0, r1, r2, r3, r4, r5, r6 とする。また，それぞれの円の曲率 1/ri を ki とする。ki = 1/ri である。
1. デカルトの円定理を用いる場合
デカルトの円定理を用いて，互いに内接・外接する 4 個の円のうちの 3 個の円の半径から 4 番目の円の半径を求める。
まず最初に，天円，地円，人円が内接する大円の半径を求める。
include("julia-source.txt");
using SymPy
@syms r0, r1, r2, r3, r4, r5, r6
k0 = 1/r1 + 1/r2 + 1/r3 - 2sqrt(1/(r1*r2) + 1/(r2*r3) + 1/(r3*r1))
r0 = 1/k0
r0 |> println
1/(-2*sqrt(1/(r2*r3) + 1/(r1*r3) + 1/(r1*r2)) + 1/r3 + 1/r2 + 1/r1)
天円，地円，人円の直径がそれぞれ 14 寸，7 寸，3 寸のとき，天円，地円，人円が内接する大円の半径は 21 寸である。
r0(r1 => 14//2, r2 => 7//2, r3 => 3//2) |> println
-21
地円，人円，乙円と，それらが内接する外円の半径（曲率）の関係式から，乙円の半径を求める。
式は複雑になるが，式に天円，地円，人円の数値を代入すると，半径は 0.6 であることがわかる。
あるいは最初から数値演算をしてもよい。
k5 = 1/r2 + 1/r3 + 1/r0 + 2sqrt(1/(r2*r3) + 1/(r3*r0) + 1/(r0*r2))
r5 = 1/k5
r5(r1 => 14//2, r2 => 7//2, r3 => 3//2) |> println
3/5
以下同様に，甲円，丙円の半径も求めることができる。
k4 = 1/r1 + 1/r3 + 1/r0 + 2sqrt(1/(r1*r3) + 1/(r3*r0) + 1/(r0*r1))
r4 = 1/k4
r4(r1 => 14//2, r2 => 7//2, r3 => 3//2) |> println
21/26
k6 = 1/r1 + 1/r2 + 1/r3 + 2sqrt(1/(r1*r2) + 1/(r2*r3) + 1/(r3*r1))
r6 = 1/k6
r6(r1 => 14//2, r2 => 7//2, r3 => 3//2) |> println
21/47
r0 = 21; r4 = 21/26 ≒ 0.807692; r5 = 3/5 = 0.6; r6 = 21/47 ≒ 0.446809
2. 二円の内接・外接関係式による場合
前節の円定理では，半径のみが求まる。図を描く場合には，それぞれの円の中心座標を求めるために二円の内接・外接関係式による連立方程式を解く必要がある。
大円，天円，地円，人円，甲円，乙円，丙円の半径と中心座標を ri, (xi, yi) とする。
パラメータを削減するために，外円の中心座標は (x0, y0) = (0, 0)，人円の中心座標は (x2, y2) = (0, r0 - r3) とする。また，r0 は デカルトの円定理から決定する。
using SymPy
@syms r0, r1, x1, y1, r2, x2, y2, r3,
r4, x4, y4, r5, x5, y5, r6, x6, y6
x3 = 0
y3 = r0 - r3
eq1 = x1^2 + y1^2 - (r0 - r1)^2
eq2 = x2^2 + y2^2 - (r0 - r2)^2
eq3 = (x3 - x1)^2 + (y3 - y1)^2 - (r3 + r1)^2
eq4 = x2^2 + (y3 - y2)^2 - (r2 + r3)^2
eq5 = x4^2 + y4^2 - (r0 - r4)^2
eq6 = (x4 - x1)^2 + (y4 - y1)^2 - (r4 + r1)^2
eq7 = (x4 - x3)^2 + (y4 - y3)^2 - (r4 + r3)^2
eq8 = x5^2 + y5^2 - (r0 - r5)^2
eq9 = (x5 - x2)^2 + (y5 - y2)^2 - (r5 + r2)^2
eq10 = (x5 - x3)^2 + (y5 - y3)^2 - (r5 + r3)^2
eq11 = (x6 - x1)^2 + (y6 - y1)^2 - (r6 + r1)^2
eq12 = (x6 - x2)^2 + (y6 - y2)^2 - (r6 + r2)^2
eq13 = (x6 - x3)^2 + (y6 - y3)^2 - (r6 + r3)^2;
using NLsolve
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
v = r.zero
end
return Float64.(v), r.f_converged
end;
function H(u)
(x1, y1, x2, y2, r4, x4, y4, r5, x5, y5, r6, x6, y6) = u
return [
x1^2 + y1^2 - (r0 - r1)^2, # eq1
x2^2 + y2^2 - (r0 - r2)^2, # eq2
x1^2 - (r1 + r3)^2 + (r0 - r3 - y1)^2, # eq3
x2^2 - (r2 + r3)^2 + (r0 - r3 - y2)^2, # eq4
x4^2 + y4^2 - (r0 - r4)^2, # eq5
-(r1 + r4)^2 + (-x1 + x4)^2 + (-y1 + y4)^2, # eq6
x4^2 - (r3 + r4)^2 + (-r0 + r3 + y4)^2, # eq7
x5^2 + y5^2 - (r0 - r5)^2, # eq8
-(r2 + r5)^2 + (-x2 + x5)^2 + (-y2 + y5)^2, # eq9
x5^2 - (r3 + r5)^2 + (-r0 + r3 + y5)^2, # eq10
-(r1 + r6)^2 + (-x1 + x6)^2 + (-y1 + y6)^2, # eq11
-(r2 + r6)^2 + (-x2 + x6)^2 + (-y2 + y6)^2, # eq12
x6^2 - (r3 + r6)^2 + (-r0 + r3 + y6)^2, # eq13
]
end;
(r0, r1, r2, r3) = (21, 14/2, 7/2, 3/2)
iniv = BigFloat[-5.395, 12.879, 4.33, 16.935, 0.8, -2.24, 20.06, 0.6, 2.2, 19.7, 0.54, 0.34, 17.5]
res = nls(H, ini=iniv)
([-5.384615384615385, 12.923076923076923, 4.3076923076923075, 16.96153846153846, 0.8076923076923077, -2.2366863905325443, 20.068047337278106, 0.6, 1.9384615384615385, 20.307692307692307, 0.44680851063829785, 0.41243862520458263, 17.597381342062192], true)
function draw(more)
pyplot(size=(1000, 1000), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r0, r1, r2, r3) = (21, 14/2, 7/2, 3/2)
(x1, y1, x2, y2, r4, x4, y4, r5, x5, y5, r6, x6, y6) = res[1]
@printf("r0 = %g; r1 = %g; r2 = %g; r3 = %g; r4 = %g; r5 = %g; r6 = %g\n", r0, r1, r2, r3, r4, r5, r6)
x3 = 0
y3 = r0 - r3
plot()
circlef(0, 0, r0, :cyan3)
circlef(x1, y1, r1)
circlef(x2, y2, r2, :green)
circlef(x3, y3, r3, :blue)
circlef(x4, y4, r4, :orange)
circlef(x5, y5, r5, :tomato)
circlef(x6, y6, r6, :darkslateblue)
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(x1, y1, "天", :white, :center, :vcenter, mark=false)
point(x2, y2, "地", :white, :center, :vcenter, mark=false)
point(x3, y3, "人", :white, :center, :vcenter, mark=false)
point(x4, y4, "甲", :white, :center, :vcenter, mark=false)
point(x5, y5, "乙", :white, :center, :vcenter, mark=false)
point(x6, y6, "丙", :white, :center, :vcenter, mark=false)
end
end;
draw(true)
