算額（その1180）
時岡郁夫：趣味の数学問題集
http://i4.gmobb.jp/tokioka/tokioka_mondai/index.html
キーワード：円8個，外円，デカルトの円定理
A問題 60. 互いに外接している3つの円の中心をO1，O2，O3，半径を r1=2，r2=3，r3=4とする。
いくつかの円を次のように定義する。
1. 円O1，O2，O3のすべてと外接する円の中心をO0，半径をr0，
1. 円O1，O2，O3のすべてと内接する円の中心をO4，半径をr4，
1. 円O2，O3 に外接し，円O4 に内接する円のうち，円O1 と異なる円の中心をO5，半径をr5，
1. 円O1，O3 に外接し，円O4 に内接する円のうち，円O2 と異なる円の中心をO6，半径をr6，
1. 円O1，O2 に外接し，円O4 に内接する円のうち，円O3 と異なる円の中心をO7，半径をr7，
1. 円O3，O5 に外接し，円O4 に内接する円のうち，円O2 と異なる円の中心をO8，半径をr8 とする。
このとき，各円の半径r0，r4～r8を求めよ。
円 O0 〜 O8 の曲率（r0 ~ r8 の逆数）を k0 〜 k8 とおき，以下のように順次決定して行く。
1. デカルトの円定理による場合
デカルトの円定理は，互いに接する4個の円の半径の間の関係式がある。これを利用すると3個の円の半径からもう一つの円の半径が決まる。
半径を求めるだけならこれで良いが，図を描くとなると，円の中心座標が必要なので，2 個の円の内接・外接関係を表す連立方程式を解く必要がある。
include("julia-source.txt");
using SymPy
(r1, r2, r3) = (2, 3, 4)
k0 = 1/r1 + 1/r2 + 1/r3 + 2sqrt(1/(r1*r2) + 1/(r2*r3) + 1/(r3*r1))
r0 = 1/k0
0.43326088256146433
他の3円を内包する円の曲率は負の値を取る。
k4 = 1/r1 + 1/r2 + 1/r3 - 2sqrt(1/(r1*r2) + 1/(r2*r3) + 1/(r3*r1))
r4 = 1/k4
-7.071558754901892
k5 = 1/r2 + 1/r3 + 1/r4 - 2sqrt(1/(r2*r3) + 1/(r3*r4) + 1/(r4*r2))
r5 = 1/k5
2.60522781835883
k6 = 1/r1 + 1/r3 + 1/r4 + 2sqrt(1/(r1*r3) + 1/(r3*r4) + 1/(r4*r1))
r6 = 1/k6
1.1314219061967086
k7 = 1/r1 + 1/r2 + 1/r4 + 2sqrt(1/(r1*r2) + 1/(r2*r4) + 1/(r4*r1))
r7 = 1/k7
0.881955860873886
k8 = 1/r3 + 1/r5 + 1/r4 + 2sqrt(1/(r3*r5) + 1/(r5*r4) + 1/(r4*r3))
r8 = 1/k8
1.534846922834955
r1 = 2; r2 = 3; r3 = 4; r4 = 7.07155875490189; r5 = 2.60522781835883; r6 = 1.13142190619671; r7 = 0.881955860873886; r8 = 1.53484692283495
2. 2 円の相互関係に関する連立方程式による場合
連立方程式は独立な部分集合ごとに順次解けばよいのであるが，どうせ SymPy の能力上の問題で一度に解けないので，一度に解ける数値解を求める。
include("julia-source.txt");
using SymPy
@syms r1, x1, y1, r2, x2, r3, x3, y3,
r5, x5, y5, r6, x6, y6,
r7, x7, y7, r8, x8, y8, R
y2 = 0
eq1 = x1^2 + y1^2 - (R - r1)^2
eq2 = x3^2 + y3^2 - (R - r3)^2
eq3 = (x1 - x2)^2 + y1^2 - (r1 + r2)^2
eq4 = (x3 - x1)^2 + (y1 - y3)^2 - (r1 + r3)^2
eq5 = (x3 - x2)^2 + y3^2 - (r2 + r3)^2
eq6 = x5^2 + y5^2 - (R - r5)^2
eq7 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2
eq8 = (x2 - x5)^2 + y5^2 - (r2 + r5)^2
eq9 = x6^2 + y6^2 - (R - r6)^2
eq10 = (x1 - x6)^2 + (y1 - y6)^2 - (r1 + r6)^2
eq11 = (x3 - x6)^2 + (y3 - y6)^2 - (r3 + r6)^2
eq12 = x7^2 + y7^2 - (R - r7)^2
eq13 = (x1 - x7)^2 + (y1 - y7)^2 - (r1 + r7)^2
eq14 = (x2 - x7)^2 + (y2 - y7)^2 - (r2 + r7)^2
eq15 = x8^2 + y8^2 - (R - r8)^2
eq16 = (x3 - x8)^2 + (y3 - y8)^2 - (r3 + r8)^2
eq17 = (x5 - x8)^2 + (y5 - y8)^2 - (r5 + r8)^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, x3, y3, r5, x5, y5, r6, x6, y6, r7, x7, y7, r8, x8, y8) = u
return [
x1^2 + y1^2 - (R - r1)^2, # eq1
x3^2 + y3^2 - (R - r3)^2, # eq2
y1^2 - (r1 + r2)^2 + (x1 - x2)^2, # eq3
-(r1 + r3)^2 + (-x1 + x3)^2 + (y1 - y3)^2, # eq4
y3^2 - (r2 + r3)^2 + (-x2 + x3)^2, # eq5
x5^2 + y5^2 - (R - r5)^2, # eq6
-(r3 + r5)^2 + (x3 - x5)^2 + (y3 - y5)^2, # eq7
y5^2 - (r2 + r5)^2 + (x2 - x5)^2, # eq8
x6^2 + y6^2 - (R - r6)^2, # eq9
-(r1 + r6)^2 + (x1 - x6)^2 + (y1 - y6)^2, # eq10
-(r3 + r6)^2 + (x3 - x6)^2 + (y3 - y6)^2, # eq11
x7^2 + y7^2 - (R - r7)^2, # eq12
-(r1 + r7)^2 + (x1 - x7)^2 + (y1 - y7)^2, # eq13
y7^2 - (r2 + r7)^2 + (x2 - x7)^2, # eq14
x8^2 + y8^2 - (R - r8)^2, # eq15
-(r3 + r8)^2 + (x3 - x8)^2 + (y3 - y8)^2, # eq16
-(r5 + r8)^2 + (x5 - x8)^2 + (y5 - y8)^2, # eq17
]
end;
(r1, r2, r3) = (2, 3, 4)
R = 7.071558754901892
iniv = BigFloat[-2.2, 4.6, -4.1, 2.8, 1.2, 2.6, -0.6, -4.4,
1.1, 0.8, 5.9, 0.88, -4.9, 3.8, 1.5, 3.5, -4.2]
res = nls(H, ini=iniv)
r1 = 2; r2 = 3; r3 = 4; r4 = 7.07155875490189; r5 = 2.60522781835883; r6 = 1.13142190619671; r7 = 0.881955860873886; r8 = 1.53484692283495
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, r3) = (2, 3, 4)
r4 = 7.071558754901892
(x1, y1, x2, x3, y3, r5, x5, y5, r6, x6, y6, r7, x7, y7, r8, x8, y8) = res[1]
@printf("r1 = %.15g; r2 = %.15g; r3 = %.15g; r4 = %.15g; r5 = %.15g; r6 = %.15g; r7 = %.15g; r8 = %.15g\n", r1, r2, r3, r4, r5, r6, r7, r8)
plot(showaxis=false)
circlef(0, 0, r4, :cyan3)
circlef(x1, y1, r1)
circlef(x2, 0, r2, :green)
circlef(x3, y3, r3, :blue)
circlef(x5, y5, r5, :tomato)
circlef(x6, y6, r6, :orange)
circlef(x7, y7, r7, :brown3)
circlef(x8, y8, r8, :magenta)
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, " O1", :white, :left, :vcenter)
point(x2, 0, " O2", :white, :left, :vcenter)
point(x3, y3, " O3", :white, :left, :vcenter)
point(0, 0, " O4", :white, :left, :vcenter)
point(x5, y5, " O5", :white, :left, :vcenter)
point(x6, y6, " O6", :white, :left, :vcenter)
point(x7, y7, " O7", :white, :left, :vcenter)
point(x8, y8, " O8", :white, :left, :vcenter)
end
end;
