算額（その632）

算額（その632）

和算図形問題あれこれ
令和3年11月の問題2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html

直角三角形の中に，正三角形と等円 2 個を入れる。股の長さが 487 寸のとき，正三角形の一辺の長さはいかほどか。

正三角形の一辺の長さを a とする。
等円の半径と中心座標を r, (x, r), (r, y)
とおき，以下の連立方程式を解く

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms r::positive, x::positive, y::positive, a::positive,
     鈎::positive, 股::positive
eq1 = distance(0, 0, a/2, sqrt(Sym(3))a/2, r, y) - r^2
eq2 = distance(a, 0, a/2, sqrt(Sym(3))a/2, x, r) - r^2
eq3 = distance(0, 鈎, 股, 0, r, y) - r^2
eq4 = distance(0, 鈎, 股, 0, x, r) - r^2
eq5 = (sqrt(Sym(3))a/2)/(股 - a/2) - 鈎/股;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r, x, y, a, 鈎))

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 v, r.f_converged
end;

function H(u)
   (r, x, y, a, 鈎) = u
   return [
       -r^2 + (3*r/4 - sqrt(3)*y/4)^2 + (-sqrt(3)*r/4 + y/4)^2,  # eq1
       -r^2 + (r - sqrt(3)*(a + sqrt(3)*r - x)/4)^2 + (-3*a/4 + sqrt(3)*r/4 + 3*x/4)^2,  # eq2
       -r^2 + (r - 股*(r*股 - y*鈎 + 鈎^2)/(股^2 + 鈎^2))^2 + (y - 鈎*(-r*股 + y*鈎 + 股^2)/(股^2 + 鈎^2))^2,  # eq3
       -r^2 + (r - 鈎*(r*鈎 - x*股 + 股^2)/(股^2 + 鈎^2))^2 + (x - 股*(-r*鈎 + x*股 + 鈎^2)/(股^2 + 鈎^2))^2,  # eq4
       sqrt(3)*a/(2*(-a/2 + 股)) - 鈎/股,  # eq5
   ]
end;

股 = 487
iniv = BigFloat[56, 295, 209, 263, 312]
res = nls(H, ini=iniv

   (BigFloat[56.11416043991601035847741706383402473599898003302003798949966598270036579889394, 295.3980059252043265908683636932107368735334613723332636723610589520722219482431, 209.4208977858381008824304265789438848352208122172291274041824122497419576667868, 263.0004802898689700383037446882855378574200007464121970616614846371699418629093, 312.0159697201720972932632403222578260473761673054778157588165017227055663499631], true)

股の長さが 487 寸のとき，正三角形の一辺の長さはほぼ 263 寸である。

その他のパラメータは以下の通り。
r = 56.1142;  x = 295.398;  y = 209.421;  a = 263;  鈎 = 312.016

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   股 = 487
   (r, x, y, a, 鈎) = res[1]
   @printf("正三角形の一辺の長さ = %.10g;  r = %g;  x = %g;  y = %g;  a = %g;  鈎 = %g\n", a, r, x, y, a, 鈎)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:gray80, lw=0.5)
   plot!([0, a, a/2, 0], [0, 0, √3a/2, 0], color=:gray80, lw=0.5)
   circle(r, y, r, :blue)
   circle(x, r, r, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:gray80, lw=0.5)
       hline!([0], color=:gray80, lw=0.5)
       point(a/2, √3a/2, "(a/2,√3a/2)")
       point(x, r, "(x,r)", :blue, :center, delta=-delta/2)
       point(r, y, "(r,y)", :blue, :center, delta=-delta/2)
       point(股, 0, "股", :black, :left, :bottom, delta=delta/2)
       point(a, 0, "a ", :black, :right, :bottom, delta=delta/2)
       point(0, 鈎, "鈎", :black, :left, :bottom, delta=delta/2)
   end
end;