算額（その492）

算額（その492）

新潟県長岡市与板町本与板 与板八幡宮 文化5年(1808)
http://www.wasan.jp/niigata/yoitahatiman3.html
涌田和芳，外川一仁: 与板八幡宮の紛失算額（3），長岡工業高等専門学校研究紀要，第49巻，2013.
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_46-50/vol_49/49_01wakuta.pdf

鈎股弦の中に，大弧，小弧，円が入っている。大弧は股と弦の隅で股に接し，鈎と弦の隅で大弧に接している。円は 2 つの弧と股に接している。鈎と股が 28.83 寸，38.44 寸のとき円の直径はいかほどか。

鈎，股の長さを 鈎，股
大弧の半径と中心座標を r1, (0, r1)
小弧の半径と中心座標を r2, (x2, 鈎/2)
円の半径と中心座標を r3, (x3, r3)
とおき，以下の連立方程式を解く。
鈎，股は後に具体的な数値を与えるので，連立方程式の解の式には鈎，股は変数名として残っている。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x2::positive, r3::positive, x3::positive,
     鈎::positive, 股::positive;
x2 = 股 + sqrt(r2^2 - (鈎//2)^2)
eq1 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2
eq2 = (x2 - x3)^2 + (鈎//2 - r3)^2 - (r2 + r3)^2
eq3 = x2^2 + (r1 - 鈎//2)^2 - (r1 + r2)^2
eq4 = 股^2 + (r1 - 鈎)^2 - r1^2
res = solve([eq1, eq2, eq3, eq4], (r1, r2, r3, x3))

   2-element Vector{NTuple{4, Sym}}:
    ((股^2 + 鈎^2)/(2*鈎), 鈎*(股^2 + 鈎^2)/(2*(股 - 鈎)*(股 + 鈎)), 股^2*鈎*(鈎*sqrt((股^8*鈎^2 + 股^6*(2*鈎^2*sqrt(1/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + 1)*(股^4 - 2*股^2*鈎^2 + 鈎^4) - 2*股^4*鈎^6 - 股^2*鈎^4*(2*鈎^2*sqrt(1/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + 1)*(股^4 - 2*股^2*鈎^2 + 鈎^4) + 鈎^10)/(股^4 - 2*股^2*鈎^2 + 鈎^4)) - (股 - 鈎)*(股 + 鈎)*(股^2 + 鈎^2)*(鈎^2*sqrt(1/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + 1))^2/(2*(股^2 + 鈎^2)*(-股^4 + 股^2*鈎^2 + 鈎^4)^2), 股*(鈎*sqrt((股^8*鈎^2 + 股^6*(2*鈎^2*sqrt(1/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + 1)*(股^4 - 2*股^2*鈎^2 + 鈎^4) - 2*股^4*鈎^6 - 股^2*鈎^4*(2*鈎^2*sqrt(1/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + 1)*(股^4 - 2*股^2*鈎^2 + 鈎^4) + 鈎^10)/(股^4 - 2*股^2*鈎^2 + 鈎^4)) - (股 - 鈎)*(股 + 鈎)*(股^2 + 鈎^2)*(鈎^2*sqrt(1/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + 1))/(-股^4 + 股^2*鈎^2 + 鈎^4))
    ((股^2 + 鈎^2)/(2*鈎), 鈎*(股^2 + 鈎^2)/(2*(股 - 鈎)*(股 + 鈎)), 股^2*鈎*(鈎*sqrt((股^8*鈎^2 + 股^6*(2*鈎^2*sqrt(1/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + 1)*(股^4 - 2*股^2*鈎^2 + 鈎^4) - 2*股^4*鈎^6 - 股^2*鈎^4*(2*鈎^2*sqrt(1/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + 1)*(股^4 - 2*股^2*鈎^2 + 鈎^4) + 鈎^10)/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + (股 - 鈎)*(股 + 鈎)*(股^2 + 鈎^2)*(鈎^2*sqrt(1/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + 1))^2/(2*(股^2 + 鈎^2)*(-股^4 + 股^2*鈎^2 + 鈎^4)^2), -股*(鈎*sqrt((股^8*鈎^2 + 股^6*(2*鈎^2*sqrt(1/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + 1)*(股^4 - 2*股^2*鈎^2 + 鈎^4) - 2*股^4*鈎^6 - 股^2*鈎^4*(2*鈎^2*sqrt(1/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + 1)*(股^4 - 2*股^2*鈎^2 + 鈎^4) + 鈎^10)/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + (股 - 鈎)*(股 + 鈎)*(股^2 + 鈎^2)*(鈎^2*sqrt(1/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + 1))/(-股^4 + 股^2*鈎^2 + 鈎^4))

2 組の解が得られるが，最初のものが適解である。

術では円の直径は「股/(鈎/sqrt(鈎^2 + 股^2) + 股/鈎)^2」と簡約化されているが，SymPy ではとても長い式のままである。

半径 = 股^2*鈎*(鈎*sqrt((股^8*鈎^2 + 股^6*(2*鈎^2*sqrt(1/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + 1)*(股^4 - 2*股^2*鈎^2 + 鈎^4) - 2*股^4*鈎^6 - 股^2*鈎^4*(2*鈎^2*sqrt(1/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + 1)*(股^4 - 2*股^2*鈎^2 + 鈎^4) + 鈎^10)/(股^4 - 2*股^2*鈎^2 + 鈎^4)) - (股 - 鈎)*(股 + 鈎)*(股^2 + 鈎^2)*(鈎^2*sqrt(1/(股^4 - 2*股^2*鈎^2 + 鈎^4)) + 1))^2/(2*(股^2 + 鈎^2)*(-股^4 + 股^2*鈎^2 + 鈎^4)^2)

鈎, 股 が 28.83, 38.44のときには，数値としては以下のようになる。

names = ("r1", "r2", "r3", "x3")
for (i, name) in enumerate(names)
   println(name, " = ", res[1][i](鈎 => 28.83, 股 => 38.44).evalf())
end

   r1 = 40.0416666666667
   r2 = 51.4821428571429
   r3 = 6.00000000000000
   x3 = 31.0000000000000

鈎, 股 が 28.83, 38.44のときには，円の直径は 12 である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股) = (28.83, 38.44)
   (r1, r2, r3, x3) = (40.0416666666667, 51.4821428571429, 6.000000000000, 31.0000000000000)
   x2 = 股 + sqrt(r2^2 - (鈎/2)^2)
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", r1, r2, x2, r3, x3)
   @printf("円の直径は %g\n", 2r3)
   plot([0, 股, 股, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
   circle(0, r1, r1, :red)
   circle(x2, 鈎/2, r2, :magenta)
   circle(x3, r3, r3, :green)
   plot!(xlims=(0, x2), ylims=(0, r1))
   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=:black, lw=0.5)
       point(0, r1, " r1", :red)
       point(x2, 鈎/2, "(x2,鈎/2)", :magenta)
       point(x3, r3, "(x3,r3)", :green, :right)
       point(股, 鈎, " (股,鈎)", :blue, :left, :vcenter)
   else
       plot!(showaxis=false)
   end
end;