算額（その876）

算額（その876）

新潟県長岡市 蒼柴神社 亨和元年(1801)3月
http://www.wasan.jp/niigata/aoshi.html

涌田和芳，外川一仁: 長岡蒼柴神社の算額，長岡工業高等専門学校研究紀要，第42巻，第2号(2006)
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_41-45/vol_42_2/42_2_1wakuta.pdf

外円内に，弦と 2 本の斜線を引き，区分された領域に甲円と乙円を 2 個ずつ入れる。大円の直径が 521 寸で，弦と矢の長さの差を最大にするとき，乙円の直径を求めよ。
注：矢（し）とは，弓形の孤の中点から弦におろした垂線
外円の半径と中心座標を R, (0, 0)
弦，矢の長さおよびその差を（そのまま）弦，矢，弦矢差
甲円の半径と中心座標を r1, (0, R - r1), (0, r - 矢 + r1)
乙円の半径と中心座標を r2, (x2, y2)
とおき以下の連立方程式を解く。

まず，弦と矢の長さの差が最大になるときの矢と弦を決定する。

include("julia-source.txt");

using SymPy
@syms R::positive, 弦::positive, 矢::positive, x::positive,
     r1::positive, r2::positive, x2::positive, y2::positive, d
R = 521//2
弦 = 2sqrt(R^2 - (R - 矢)^2)
弦矢差 = 弦 - 矢;

pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(弦矢差, xlims=(9, R), xlabel="矢", ylabel="弦 - 矢")

矢が 150 前後のときに弦矢差が最大になることがわかる。

最大値を取るときの矢の値を求めるには，導関数を求め，導関数が 0 になるときの値を求めればよい。

g = diff(弦矢差, 矢);
g |> println

   2*(521/2 - 矢)/sqrt(271441/4 - (521/2 - 矢)^2) - 1

mx_a = solve(g, 矢)[1]
mx_a |> factor |> println
mx_a.evalf() |> println

   -521*(-5 + sqrt(5))/10
   144.000858372261

弦矢差(矢 => mx_a.evalf()) |> println

   321.995708138695

弦(矢 => mx_a.evalf()) |> println

   465.996566510956

矢が 521/2 - 521*sqrt(5)/10 = 144.000858372261 のとき，弦と矢の差が最大値 321.995708138695 になる。ちなみにそのときの弦は 465.996566510956 である。

465.996566510956 - 144.000858372261, 321.995708138695

   (321.99570813869497, 321.995708138695)

弦と矢の差が最大となるときの矢，弦が決まったので，その状況における甲円と乙円の半径を求める。

using SymPy
@syms R::positive, 弦::positive, 矢::positive, x::positive,
     r1::positive, r2::positive, x2::positive, y2::positive, d
矢 = (5 - sqrt(Sym(5)))R/5
弦 = 2sqrt(R^2 - (R - 矢)^2)
eq1 = x2^2 + y2^2 - (R - r2)^2
eq2 = dist(x, sqrt(R^2 - x^2), -弦/2, R - 矢, 0, R - r1) - r1^2
eq3 = dist(x, sqrt(R^2 - x^2), -弦/2, R - 矢, 0, R - 矢 + r1) - r1^2
eq4 = dist(x, sqrt(R^2 - x^2), -弦/2, R - 矢, x2, y2) - r2^2
eq5 = dist(-x, sqrt(R^2 - x^2), 弦/2, R - 矢, x2, y2) - r2^2;

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)
   (r1, r2, x2, y2, x) = u
   return [
x2^2 + y2^2 - (R - r2)^2,  # eq1
-r1^2 + (-x - (-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))*(-x*(-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (R - r1 - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2)))/((-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2 + (R - r1 - sqrt(R^2 - x^2) - (-x*(-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (R - r1 - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2)))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))/((-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2,  # eq2
-r1^2 + (-x - (-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))*(-x*(-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R + r1 - sqrt(R^2 - x^2)))/((-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2 + (-R*(5 - sqrt(5))/5 + R + r1 - sqrt(R^2 - x^2) - (-x*(-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R + r1 - sqrt(R^2 - x^2)))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))/((-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2,  # eq3
-r2^2 + (-x + x2 - (-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))*((-x + x2)*(-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (y2 - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2)))/((-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2 + (y2 - sqrt(R^2 - x^2) - ((-x + x2)*(-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (y2 - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2)))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))/((-x - sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2,  # eq4
-r2^2 + (x + x2 - (x + sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))*((x + x2)*(x + sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (y2 - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2)))/((x + sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2 + (y2 - sqrt(R^2 - x^2) - ((x + x2)*(x + sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2)) + (y2 - sqrt(R^2 - x^2))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2)))*(-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))/((x + sqrt(R^2 - (-R*(5 - sqrt(5))/5 + R)^2))^2 + (-R*(5 - sqrt(5))/5 + R - sqrt(R^2 - x^2))^2))^2,  # eq5
   ]
end;

R = 521//2
iniv = BigFloat[35, 37, 122, 188, 127]
res = nls(H, ini=iniv)

   ([35.179523802100135, 36.00021459306524, 121.93467698217249, 188.49957081386952, 126.65410684822777], true)

外円の直径が 521 寸のとき，乙円の直径は 72 寸有奇である。

その他のパラメータは以下のとおりである。
r1 = 35.1795;  r2 = 36.0002;  x2 = 121.935;  y2 = 188.5;  x = 126.654

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 521//2
   矢 = (5 - √5)R/5
   弦 = 2sqrt(R^2 - (R - 矢)^2)
   (r1, r2, x2, y2, x) = res[1]
   @printf("弦 = %g;  矢 = %g;  弦と矢の差 = %g\n", 弦, 矢, 弦 - 矢)
   @printf("乙円の直径 = %.15g\n", 2r2)
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g;  x = %g\n", r1, r2, x2, y2, x)
   plot()
   circle(0, 0, R)
   segment(-弦/2, R - 矢, 弦/2, R - 矢, :blue)
   circle2(x2, y2, r2, :green)
   circle(0, R - r1, r1, :magenta)
   circle(0, R - 矢 + r1, r1, :magenta)
   y = sqrt(R^2 - x^2)
   segment(x, y, -弦/2, R - 矢)
   segment(-x, y, 弦/2, R - 矢)
   if more        
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, R - 矢, "R-矢", :blue, :center, delta=-delta/2)
       point(x, y, "(x,y)", :green, :left, :bottom, delta=delta/2)
       point(弦/2, R - 矢, "弦/2 ", :blue, :right, delta=-delta/2)
       point(0, R - r1, "甲円:r1,(0,R-r1)", :black, :center, :bottom, delta=delta/2)
       point(0, R - 矢/2, "R-矢/2   ", :black, :right, :vcenter)
       point(0, R - 矢 + r1, "甲円:r1,(0,R-矢+r1)", :black, :center, delta=-delta/2)
       point(x2, y2, "乙円:r2,(x2,y2)", :black, :center, delta=-delta/2)
       point(0, R, "R", :red, :center, :bottom, delta=delta/2)
 end
end;