算額（その969）
一七 大里郡岡部村岡 稲荷社 文化13年(1816)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
横倒しされた台形の中に 2 本の斜線を入れ，区分された領域に 甲円，乙円，丙円を入れる。大頭（台形の下底），甲円の直径，乙円の直径がそれぞれ 16 寸，12 寸，8 寸のとき，丙円の直径はいかほどか。
台形の上底，下底，高さを h2, h1, a
斜線と台形の辺の交点座標を (b, 0), (c, 0), (d, (h1 - h2)*d/a)
甲円の半径と中心座標を r1, (r1, r1)
乙円の半径と中心座標を r2, (a -r2, r2)
丙円の半径と中心座標を r3, (x3, r3)
とおき，以下の連立方程式の数値解を求める。
include("julia-source.txt");
using SymPy
@syms a::positive, b::positive,
c::positive, d::positive, h1::positive,
h2::positive, r1::positive, r2::positive,
r3::positive, x3::positive
eq1 = dist(0, h1, a, h2, a - r2, r2) - r2^2
eq2 = dist(b, 0, d, h1 - (h1 - h2)*d/a, a - r2, r2) - r2^2
eq3 = dist(b, 0, d, h1 - (h1 - h2)*d/a, x3, r3) - r3^2
eq4 = dist(b, 0, d, h1 - (h1 - h2)*d/a, r1, r1) - r1^2
eq5 = dist(0, h1, c, 0, a - r2, r2) - r2^2
eq6 = dist(0, h1, c, 0, x3, r3) - r3^2
eq7 = dist(0, h1, c, 0, r1, r1) - r1^2
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r3, x3, h2, a, b, c, d))
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)
(r3, x3, h2, a, b, c, d) = u
return [
-r2^2 + (a - a*(a*(a - r2) + (-h1 + h2)*(-h1 + r2))/(a^2 + (-h1 + h2)^2) - r2)^2 + (-h1 + r2 - (-h1 + h2)*(a*(a - r2) + (-h1 + h2)*(-h1 + r2))/(a^2 + (-h1 + h2)^2))^2, # eq1
-r2^2 + (r2 - (h1 - d*(h1 - h2)/a)*(r2*(h1 - d*(h1 - h2)/a) + (-b + d)*(a - b - r2))/((-b + d)^2 + (h1 - d*(h1 - h2)/a)^2))^2 + (a - b - r2 - (-b + d)*(r2*(h1 - d*(h1 - h2)/a) + (-b + d)*(a - b - r2))/((-b + d)^2 + (h1 - d*(h1 - h2)/a)^2))^2, # eq2
-r3^2 + (r3 - (h1 - d*(h1 - h2)/a)*(r3*(h1 - d*(h1 - h2)/a) + (-b + d)*(-b + x3))/((-b + d)^2 + (h1 - d*(h1 - h2)/a)^2))^2 + (-b + x3 - (-b + d)*(r3*(h1 - d*(h1 - h2)/a) + (-b + d)*(-b + x3))/((-b + d)^2 + (h1 - d*(h1 - h2)/a)^2))^2, # eq3
-r1^2 + (r1 - (h1 - d*(h1 - h2)/a)*(r1*(h1 - d*(h1 - h2)/a) + (-b + d)*(-b + r1))/((-b + d)^2 + (h1 - d*(h1 - h2)/a)^2))^2 + (-b + r1 - (-b + d)*(r1*(h1 - d*(h1 - h2)/a) + (-b + d)*(-b + r1))/((-b + d)^2 + (h1 - d*(h1 - h2)/a)^2))^2, # eq4
-r2^2 + (a - c*(c*(a - r2) - h1*(-h1 + r2))/(c^2 + h1^2) - r2)^2 + (-h1 + h1*(c*(a - r2) - h1*(-h1 + r2))/(c^2 + h1^2) + r2)^2, # eq5
-r3^2 + (-c*(c*x3 - h1*(-h1 + r3))/(c^2 + h1^2) + x3)^2 + (-h1 + h1*(c*x3 - h1*(-h1 + r3))/(c^2 + h1^2) + r3)^2, # eq6
-r1^2 + (-c*(c*r1 - h1*(-h1 + r1))/(c^2 + h1^2) + r1)^2 + (-h1 + h1*(c*r1 - h1*(-h1 + r1))/(c^2 + h1^2) + r1)^2, # eq7
]
end;
(h1, r1, r2) = (16, 12/2, 8/2)
iniv = BigFloat[2.3, 20.8, 7.11, 35.0, 7.0, 30.0, 30.83]
res = nls(H, ini=iniv)
([2.3, 20.8, 7.111111111111111, 35.0, 7.0, 30.0, 30.829787234042552], true)
大頭 = 16 寸, 甲円径 = 12 寸, 乙円径 = 8 寸のとき，丙円径 = 4.6 寸である。
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(h1, r1, r2) = (16, 12/2, 8/2)
(r3, x3, h2, a, b, c, d) = [2.3, 20.8, 7.111111111111111, 35.0, 7.0, 30.0, 30.829787234042552]
@printf("大頭 = %g, 甲円径 = %g, 乙円径 = %g のとき，丙円径 = %g\n", h1, 2r1, 2r2, 2r3)
plot([0, a, a, 0, 0], [0, 0, h2, h1, 0], color=:blue, lw=0.5)
circle(r1, r1, r1)
circle(a - r2, r2, r2, :green)
circle(x3, r3, r3, :magenta)
segment(0, h1, c, 0, :black)
segment(b, 0, d, h1 -(h1 - h2)*d/a, :black)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:black, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
point(0, h1, " h1", :blue, :left, :bottom, delta=delta)
point(a, h2, " (a,h2)", :blue, :center, :bottom, delta=2delta)
point(b, 0, "b", :black, :center, delta=-2delta)
point(c, 0, "c", :black, :center, delta=-2delta)
point(d, h1 - (h1 - h2)*d/a, "(d,h1-(h1-h2)*d/a)", :black, :center, :bottom, delta=5delta)
point(r1, r1, "甲円:r1,(r1,r1)", :red, :center, delta=-2delta)
point(a - r2, r2, "乙円:r2\n(a-r2,r2)", :green, :center, delta=-2delta)
point(x3, r3, "丙円:r3,(x3,r3)", :black, :center, delta=-2delta)
plot!(ylims=(-8delta, h1 + 3delta))
end
end;
