算額（その185）

算額（その185）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（85）
長野県中野市田上 田上観音堂 文化6年(1809)

第2問 台形の隔斜（対角線）で区切られた領域に甲，乙，丙の円を入れる。甲円の直径が61寸6分，乙円の直径が20寸，丙円の直径が31寸5分のとき，小頭（台形の上底）はいかほどか。

以下のように記号を定め方程式を解く。
台形の右下隅 A，右上隅 b， の座標をそれぞれ (a, 0)，(b, y) と置く。

甲円，乙円，丙円の中心から隔斜（対角線）までの距離が円の半径に等しいという連立方程式を解く。

小数付きを避けるために 10 倍する。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
   p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
   l = sympy.Line(p1, p2)
   l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

@syms a::positive, b::positive, x1::positive, y1::positive, y::positive;

eq1 = distance(b, y, a, 0, x1, y1) - 315^2
eq2 = distance(-b, y, a, 0, x1, y1) - 315^2
eq3 = distance(-a, 0, b, y, x1, y1) - 315^2
eq4 = distance(-a, 0, b, y, 0, y - 200) - 200^2
eq5 = distance(-b, y, a, 0, 0, 616) - 616^2;

# res = solve([eq1, eq2, eq3, eq4, eq5], (a, b, x1, y1, y))

nlsolve() で解く。

println(eq1, ",")
println(eq2, ",")
println(eq3, ",")
println(eq4, ",")
println(eq5, ",")

   (x1 - (a^2*x1 - 2*a*b*x1 + a*y^2 - a*y*y1 + b^2*x1 + b*y*y1)/(a^2 - 2*a*b + b^2 + y^2))^2 + (-y*(a^2 - a*b - a*x1 + b*x1 + y*y1)/(a^2 - 2*a*b + b^2 + y^2) + y1)^2 - 99225,
   (x1 - (x1*(a^2 + 2*a*b + b^2 + y^2) - y*(-a*y + a*y1 + b*y1 + x1*y))/(a^2 + 2*a*b + b^2 + y^2))^2 + (-y*(a^2 + a*b - a*x1 - b*x1 + y*y1)/(a^2 + 2*a*b + b^2 + y^2) + y1)^2 - 99225,
   (x1 - (a^2*x1 + 2*a*b*x1 - a*y^2 + a*y*y1 + b^2*x1 + b*y*y1)/(a^2 + 2*a*b + b^2 + y^2))^2 + (-y*(a^2 + a*b + a*x1 + b*x1 + y*y1)/(a^2 + 2*a*b + b^2 + y^2) + y1)^2 - 99225,
   y^2*(-200*a + b*y - 200*b)^2/(a^2 + 2*a*b + b^2 + y^2)^2 + (-y*(a^2 + a*b + y^2 - 200*y)/(a^2 + 2*a*b + b^2 + y^2) + y - 200)^2 - 40000,
   y^2*(a*y - 616*a - 616*b)^2/(a^2 + 2*a*b + b^2 + y^2)^2 + (-y*(a^2 + a*b + 616*y)/(a^2 + 2*a*b + b^2 + y^2) + 616)^2 - 379456,

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=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (a, b, x1, y1, y) = u
   return [
       (x1 - (x1*(a^2 - 2*a*b + b^2 + y^2) + y*(a*y - a*y1 + b*y1 - x1*y))/(a^2 - 2*a*b + b^2 + y^2))^2 + (-y*(a^2 - a*b - a*x1 + b*x1 + y*y1)/(a^2 - 2*a*b + b^2 + y^2) + y1)^2 - 99225,
       (x1 - (x1*(a^2 + 2*a*b + b^2 + y^2) - y*(-a*y + a*y1 + b*y1 + x1*y))/(a^2 + 2*a*b + b^2 + y^2))^2 + (-y*(a^2 + a*b - a*x1 - b*x1 + y*y1)/(a^2 + 2*a*b + b^2 + y^2) + y1)^2 - 99225,
       (x1 - (a^2*x1 + 2*a*b*x1 - a*y^2 + a*y*y1 + b^2*x1 + b*y*y1)/(a^2 + 2*a*b + b^2 + y^2))^2 + (-y*(a^2 + a*b + a*x1 + b*x1 + y*y1)/(a^2 + 2*a*b + b^2 + y^2) + y1)^2 - 99225,
       y^2*(-200*a + b*y - 200*b)^2/(a^2 + 2*a*b + b^2 + y^2)^2 + (-y*(a^2 + a*b + y^2 - 200*y)/(a^2 + 2*a*b + b^2 + y^2) + y - 200)^2 - 40000,
       y^2*(a*y - 616*a - 616*b)^2/(a^2 + 2*a*b + b^2 + y^2)^2 + (-y*(a^2 + a*b + 616*y)/(a^2 + 2*a*b + b^2 + y^2) + 616)^2 - 379456,
   ]
end;

iniv = [big"3172", 2187.0, 770, 1586, 2649]
res = nls(H, ini=iniv)
println(res)
   (BigFloat[1848.000000000000000000070344227075033039909771395411031613018606702915177638875, 599.9999999999999999996728744427154150093629200433208589556725712328973843043104, 524.9999999999999999996719021835952608434605926603387970940837613822641515845599, 1385.999999999999999999661107977232676125151761705765721400241291345685648773989, 1835.999999999999999999496031313197167629938625864150348051177861910483984084011], true)

b = 600 であるから，元の単位でいえば小頭は 60 寸である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=:black, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
   end
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, x1, y1, y) = res[1]  # [3172, 2187.0, 770, 1586, 2649]
   @printf("a = %6.3f; b = %6.3f;  x1 = %6.3f; y1 = %6.3f; y = %6.3f\n", a, b, x1, y1, y)
   plot([a, b, -b, -a, a], [0, y, y, 0, 0], color=:black, lw=0.5)
   circle(0, 616, 616, :green)
   circle(x1, y1, 315, :blue)
   circle(-x1, y1, 315, :blue)
   circle(0, y - 200, 200, :red)
   segment(-b, y, a, 0, :black, linewidth=0.25)
   segment(-a, 0, b, y, :black, linewidth=0.25)
   if more == true
       point(0, 616, " 甲: (0,616)", :green)
       point(0, y - 200, "乙: (0,y-200)", :red, :center)
       point(x1, y1, "丙: (x1,y1,315)", :blue, :center)
       point(b, y, "  B(b,y)", :black)
       point(a, 0, " A(a,0)", :black, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

   a = 1848.000; b = 600.000;  x1 = 525.000; y1 = 1386.000; y = 1836.000