算額（その621）

算額（その621）

和算図形問題あれこれ
令和 2 年 10 月の問題 - No.2?
『絵本工夫之錦』より
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html

垂線（中鈎）で区切られた三角形内に直径が 3 寸と 2 寸の大円と小円を入れる。2 円の直径は変えずに中鈎の長さを変えると三角形の面積も変化する。面積が最小となる中鈎の長さはいかほどか。

大円の半径と中心座標を r1, (r1, r1)
小円の半径と中心座標を r2, (-r2, r2)
中鈎の長さを h，左右の頂点の座標を (-x2, 0), (x1, 0)
とおき，以下の連立方程式を解く。それぞれのパラメータは h を含む式になる。

include("julia-source.txt");
# julia-source.txt ソース

using SymPy

@syms x1, x2, h, 弦1, 弦2, s, r1, r2

(r1, r2) = (3, 2) .// 2

eq1 = h^2 + x1^2 - 弦1^2
eq2 = h^2 + x2^2 - 弦2^2
eq3 = h + x1 - 弦1 - 2r1
eq4 = h + x2 - 弦2 - 2r2
eq5 = h*(x1 + x2)/2 - s;

res = solve([eq1, eq2, eq3, eq4, eq5], (s, x1, x2, 弦1, 弦2))

   1-element Vector{NTuple{5, Sym}}:
    (h*(2*h - 5)*(5*h - 6)/(4*(h - 3)*(h - 2)), 3*(2*h - 3)/(2*(h - 3)), 2*(h - 1)/(h - 2), (2*h^2 - 6*h + 9)/(2*(h - 3)), (h^2 - 2*h + 2)/(h - 2))

h と s の関係を図に示すと，4.7 < h < 4.8 において，面積が最小になるのがわかる。

using Plots
s = h*(2*h - 5)*(5*h - 6)/(4*(h - 3)*(h - 2))
pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(s, xlims=(4.5, 5), ylims=(19.7, 19.8), xlabel="h", ylabel="s")
hline!([19.7092342844582])
vline!([4.73850950038331])

面積 s が最小になる h を求めるには，s(h) の式を h で微分して，s(h)′ = 0 を解く（接線の傾きが 0 になるときの h を求める）。

s′ = diff(s, h)
plot(s′, xlims=(4.5, 5), xlabel="h", ylabel="s′")

res2 = solve(s′);

4 組の解が得られるが，2番目のものが適解である。得られる数値解は「虚数解」になっている。
ただ，虚数部は 「0.e-21*I」となっており，実質的に 0 であり，実数部 4.73850950038331 が求める解である。

res2[2].evalf(), real(res2[2].evalf())

   (4.73850950038331 + 0.e-21*I, 4.73850950038331)

そのときの x1, x2, 面積は 19.7092342844582 である。

h = real(res2[2].evalf())
(x1, x2) = (3*(2*h - 3)/(2*(h - 3)), 2*(h - 1)/(h - 2))
s = h*(x1 + x2)/2
s |> println

   19.7092342844582

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   h  = 4.73850950038331
   (x1, x2) = (3*(2*h - 3)/(2*(h - 3)), 2*(h - 1)/(h - 2))
   plot([-x2, x1, 0, -x2], [0, 0, h, 0], color=:black, lw=0.5)
   circle(r1, r1, r1)
   circle(-r2, r2, r2, :blue)
   segment(0, 0, 0, h, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(r1, r1, "大円:r1,(r1,r1)", :red, :center, delta=-delta)
       point(-r2, r2, "小円:r2,(-r2,r2)", :blue, :center, delta=-delta)
       point(x1, 0, "x1", :red, :center, delta=-delta)
       point(-x2, 0, "x2", :red, :center, delta=-delta)
       point(0, h, " h", :red, :left, :bottom, delta=delta/2)
       plot!(ylims=(-0.5, 5))
   end
end;