算額（その751）
三五 大宮市中釘 秋葉神社 天保11年(1840)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
埼玉県さいたま市西区中釘 秋葉神社 天保11年(1840)
山口正義：やまぶき，第20号
https://yamabukiwasan.sakura.ne.jp/ymbk20.pdf
キーワード：円5個，不等辺三角形
三角形の中に甲円 3 個，乙円と丙円を 1 個ずついれる。甲円の直径が 20 寸のとき，乙円の直径はいかほどか。
三角形の頂点を (0, 0, (x, y), (x2, 0)
甲円の半径と中心座標を r1, (x1, r1), (x1 + 2r1, r1), (x1 + 4r1, r1)
乙円の半径と中心座標を r2, (x1 + 3r1, y2)
丙円の半径と中心座標を r3, (x1 + r1, y3)
と置き，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms d, r1::positive, x1::negative,
r2::positive, y2::negative,
r3::positive, y3::positive,
x::positive, y::positive, x2::positive
eq1 = r1^2 + (y3 - r1)^2 - (r1 + r3)^2
eq2 = 4r1^2 + (y2 - y3)^2 - (r2 + r3)^2
eq3 = r1^2 + (y2 - r1)^2 - (r1 + r2)^2
eq4 = dist(0, 0, x, y, x1, r1) - r1^2
eq5 = dist(0, 0, x, y, x1 + r1, y3) - r3^2
eq6 = dist(0, 0, x, y, x1 + 3r1, y2) - r2^2
eq7 = dist(x, y, x2, 0, x1 + 3r1, y2) - r2^2
eq8 = dist(x, y, x2, 0, x1 + 4r1, r1) - r1^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, r2, y2, r3, y3, x, y, x2) = u
return [
r1^2 + (-r1 + y3)^2 - (r1 + r3)^2, # eq1
4*r1^2 - (r2 + r3)^2 + (y2 - y3)^2, # eq2
r1^2 + (-r1 + y2)^2 - (r1 + r2)^2, # eq3
-r1^2 + (r1 - y*(r1*y + x*x1)/(x^2 + y^2))^2 + (-x*(r1*y + x*x1)/(x^2 + y^2) + x1)^2, # eq4
-r3^2 + (-y*(x*(r1 + x1) + y*y3)/(x^2 + y^2) + y3)^2 + (r1 - x*(x*(r1 + x1) + y*y3)/(x^2 + y^2) + x1)^2, # eq5
-r2^2 + (-y*(x*(3*r1 + x1) + y*y2)/(x^2 + y^2) + y2)^2 + (3*r1 - x*(x*(3*r1 + x1) + y*y2)/(x^2 + y^2) + x1)^2, # eq6
-r2^2 + (-y + y*(-y*(-y + y2) + (-x + x2)*(3*r1 - x + x1))/(y^2 + (-x + x2)^2) + y2)^2 + (3*r1 - x + x1 - (-x + x2)*(-y*(-y + y2) + (-x + x2)*(3*r1 - x + x1))/(y^2 + (-x + x2)^2))^2, # eq7
-r1^2 + (r1 - y + y*(-y*(r1 - y) + (-x + x2)*(4*r1 - x + x1))/(y^2 + (-x + x2)^2))^2 + (4*r1 - x + x1 - (-x + x2)*(-y*(r1 - y) + (-x + x2)*(4*r1 - x + x1))/(y^2 + (-x + x2)^2))^2, # eq8
]
end;
r1 = 10
iniv = BigFloat[24, 15, 33, 6.4, 24, 62.5, 62.5, 78]
res = nls(H, ini=iniv)
([24.534049509916077, 14.768336246810202, 32.659887034913744, 7.084973778708188, 23.852665051142555, 63.12069962325806, 61.70734973660005, 77.04124269850617], true)
甲円の直径が 20 寸のとき，乙円の直径は 29.5367 寸である。
その他のパラメータは以下のとおりである。
甲円の直径 = 20; 乙円の直径 = 29.5367; 丙円の直径 = 14.1699
x1 = 24.534; r2 = 14.7683; y2 = 32.6599; r3 = 7.08497; y3 = 23.8527; x = 63.1207; y = 61.7073; x2 = 77.0412
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 10
(x1, r2, y2, r3, y3, 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, :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(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)", :magenta, :center, delta=-delta)
point(x1 + r1, y3, " 丙円:r3,(x1+r1,y3)", :blue, :left, :vcenter)
point(x2, 0, " x2", :black, :left, :bottom, delta=delta)
end
end;
