算額（その969）

算額（その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;