算額（その425）
長野県諏訪市中洲 諏訪大社上社 明治12年(1897)
中村信弥(1999)：算額への招待
http://www.wasan.jp/syotai/syotai.html
長方形内に，甲，乙，丙，丁，戊，己の 6 個の円がある。それぞれ長方形に内接したり，互いに外接しあっている。
己円の直径を丁円と戊円の直径で表わせ。
長方形の長辺の長さを a とする（短辺の長さは 甲円の直径に等しい）。
甲円の半径と中心座標 r1, (a - r1, r1)
乙円の半径と中心座標 r2, (r2, r2)
丙円の半径と中心座標 r3, (x3, 2r1 - r3)
丁円の半径と中心座標 r4, (a - r4, r4)
戊円の半径と中心座標 r5, (r5, r5)
己円の半径と中心座標 r6, (x6, y6)
include("julia-source.txt");
using SymPy
@syms a::positive, r1::positive, r2::positive, r3::positive,
x3::positive, r4::positive, r5::positive, r6::positive,
x6::positive, y6::positive;
eq1 = (a - r1 - r2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (a - r1 - r3)^2 + (2r1 - r3 - r1)^2 - (r1 + r3)^2
eq3 = 2(r1 - r4)^2 - (r1 + r4)^2
eq4 = (a - r1 - x6)^2 + (r1 - y6)^2 - (r1 + r6)^2
eq5 = (x3 - r2)^2 + (2r1 - r3 - r2)^2 - (r2 + r3)^2
eq6 = 2(r2 - r5)^2 - (r2 + r5)^2
eq7 = (r2 - x6)^2 + (r2 - y6)^2 - (r2 + r6)^2
eq8 = (x3 - x6)^2 + (2r1 - r3 - y6)^2 - (r3 + r6)^2
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (a, r1, r2, r3, x3, r6, x6, y6))
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 v, r.f_converged
end;
function H(u)
(a, r1, r2, r3, x3, r6, x6, y6) = u
return [
(r1 - r2)^2 - (r1 + r2)^2 + (a - r1 - r2)^2, # eq1
(r1 - r3)^2 - (r1 + r3)^2 + (a - r1 - r3)^2, # eq2
2*(r1 - r4)^2 - (r1 + r4)^2, # eq3
-(r1 + r6)^2 + (r1 - y6)^2 + (a - r1 - x6)^2, # eq4
(-r2 + x3)^2 - (r2 + r3)^2 + (2*r1 - r2 - r3)^2, # eq5
2*(r2 - r5)^2 - (r2 + r5)^2, # eq6
-(r2 + r6)^2 + (r2 - x6)^2 + (r2 - y6)^2, # eq7
-(r3 + r6)^2 + (x3 - x6)^2 + (2*r1 - r3 - y6)^2, # eq8
]
end;
(r4, r5) = (2, 1)
iniv = [big"180.0", 59, 35, 21, 55, 5, 65, 62] ./2
res = nls(H, ini=iniv);
println(res);
(BigFloat[33.97056274847714058562026469051637694283606250452337687812015685588878974154529, 11.65685424949238019520675489683879231427868750150779229270671895196292991384848, 5.828427124746190097603377448419396157139343750753896146353359475981464956924239, 5.828427124746190097603377448419396157139343750753896146353359475981464956924308, 5.82842712474619009760948182491957344489198207867602126024818280529793350621324, 1.093836321356054313599895375610329707159527028223527974012807033844387485155129, 9.563017928136325881606859521228462607119160473284264318693911918118542430249862, 11.65685424949238019520479919353753625575412910164858126097863477985370829598065], true)
r4 = 2; r5 = 1
a = 33.9706; r1 = 11.6569; r2 = 5.82843; r3 = 5.82843; x3 = 5.82843; r6 = 1.09384; x6 = 9.56302; y6 = 11.6569
己円の直径は 2.18767 であるが，術では以下のように 7.31370849898476 であり，答えが違う。中村が術を読み違えているのか，式に起こすときに誤記があるのか。
(丁, 戊) = (2, 1) .* 2
天 = 丁 * 戊
地 = 丁 + 戊
法 = 2(sqrt((4地 * 戊 + 丁^2)天) + 2戊^2) + 地*丁
法*丁/(3 + √8)天
7.31370849898476
中村は下記の x1 式を簡約化して x2 を提示し術と同じとしている。
しかし，x1 は正しいが x2 は式の記述ミスがある。正しくは x3（分母の式の一箇所の - が + でなければならない）。
@syms d, e
x1 = (3 + 2sqrt(Sym(2)))*d^2*e*((d^2 + d*e + 4e^2) - 2(d + 2e)*sqrt(d*e))/(d^4 - 2d^3*e - 7d^2*e^2 - 8d*e^3 + 16e^4)
x2 = (3 + 2sqrt(Sym(2)))*d^2*e / (d^2 + d*e + 4e^2 - 2(d + 2e)*sqrt(d*e))
x3 = (3 + 2sqrt(Sym(2)))*d^2*e / (d^2 + d*e + 4e^2 + 2(d + 2e)*sqrt(d*e))
x1(d => 4, e => 2).evalf() |> println
x2(d => 4, e => 2).evalf() |> println
x3(d => 4, e => 2).evalf() |> println
2.18767264271211
-35.4929704984046
2.18767264271211
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r4, r5) = (2, 1)
(a, r1, r2, r3, x3, r6, x6, y6) = res[1]
@printf("r4 = %g; r5 = %g\n", r4, r5)
@printf("a = %g; r1 = %g; r2 = %g; r3 = %g; x3 = %g; r6 = %g; x6 = %g; y6 = %g\n", a, r1, r2, r3, x3, r6, x6, y6)
@printf("己円の直径 = %g\n", 2r6)
plot([0, a, a, 0, 0], [0, 0, 2r1, 2r1, 0], color=:black, lw=0.5)
circle(a - r1, r1, r1)
circle(r2, r2, r2, :blue)
circle(x3, 2r1 - r3, r3, :green)
circle(a - r4, r4, r4, :magenta)
circle(r5, r5, r5, :orange)
circle(x6, y6, r6, :purple)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(a - r1, r1, " 甲円:r1,(a-r1,r1)", :red, :left, :vcenter)
point(r2, r2, " 乙円:r2,(r2,r2)", :blue, :center, :top, delta=-delta)
point(x3, 2r1 - r3, " 丙円:r3,(x3,2r1-r3)", :green, :center, :top, delta=-delta)
point(a - r4, r4, "丁円:r4,(a-r4,r4) ", :magenta, :right, :vcenter)
point(r5, r5, " 戊円:r5,(r5,r5) ", :orange, :left, :vcenter)
point(x6, y6, " 己円:r6,(x5,y5) ", :purple, :left, :vcenter)
point(a, 2r1, "(a,2r1) ", :black, :right, :top, delta=-delta)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
