算額（その1404）

算額（その1404）

三十二 一関市舞川 観福寺内地蔵堂 明治43年(1901)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円2個，直角三角形4個，等脚台形
#Julia, #SymPy, #算額, #和算

問の原文は以下のとおりである。
「今有如図半円内設等円二個従其親所洩二斜容等円一個其等円径若干問得半円径術如何」

最初にいっておく。本問は条件不足であり，解けない。

原文には曖昧な点もあるので，補足的説明も加えると以下のようになろう。
「円弧（注1）」の中に等円 2 個を容れ，等円と円弧の接点と円弧の中心を結ぶ半径を 2 本引く（注2）。2 本の半径と円弧の弦が作る三角形の中に，円弧内の等円と同じ大きさの等円が内接する。等円の直径が与えられたとき，円弧の半径を求める方法を述べよ。
注1：原文は「半円」であるが，それでは図が成り立たない。
注2：この半径は，（当然ではあるが）円弧内の等円の中心を通る。

「条件不足」という意味は，以下に示すように，「等円径若干」で与えられた等円の直径が同じでも，円弧の直径が一意に定まらないということである。

円弧の半径と中心座標を R, (0, 0)
弦と y 軸の交点座標を (0, y)
等円の半径と中心座標を r, (x, y + r), (0, y - r)
とおくと以下の 2 方程式が立つ。
特定の r が与えられたとき，図を描くために必要なのは R, x, y の 3 変数であるが条件式が 2 本なので，「条件不足」ということになる。

術は「置等円径三之得半円径」ということで，「等円の直径を 3 倍すれば半円（円弧）の径が得られる」ということであるが，算額の図を見てもそんな比率でないことは明らかである。

r = 1/2 の場合について，連立方程式を解くが，「等円の直径が同じでも，円弧の直径が一意に定まらない」ことを示すために，R を色々変えて，x, y を求める。

R = 3 のときの x, y は以下のようになる。

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

@syms x, y, r, R
r = 1/2
R = 3
eq1 = x^2 + (y + r)^2 - (R - r)^2
eq2 = dist2(0, 0, x, y + r, 0, y - r, r);
res = solve([eq1, eq2], (x, y));

8 組の虚数解が得られるが，虚数部は実質 0 である。また，x < y という条件をつけて，以下の手順で実質的に適切な解を選択・表示する。

for (x, y) in res
   if abs(imag(x)) < 1e-15 && abs(imag(y)) < 1e-15 &&
      real(x) > 0 && real(y) > 0 && real(x) < real(y)
       println("(R, x, y) = ", (R, real(x), real(y)))
   end
end

   (R, x, y) = (3, 0.953431800017409, 1.81105339676857)

等円の直径が 1 のとき，円弧の直径は 6 である。ただし，x = 0.953432，y = 1.81105

しかし，円弧の直径は一意に決まらない。たとえば，他の例として，
等円の直径が 1 のとき，円弧の直径は 5.7 である。ただし，x = 1.08171，y = 1.58624
等円の直径が 1 のとき，円弧の直径は 8 である。ただし，x = 0.721715，y = 2.92478
　　：
当然ではあるが，円弧の直径はどのような値でもよいわけではなく，一定の範囲内に収まっていなければならない。

function draw(r, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, x, y) = (3, 0.953431800017409, 1.81105339676857)
   @printf("等円の直径が %g のとき，円弧の直径は %g である。ただし，x = %g，y = %g である。\n", 2r, 2R, x, y)
   @printf("r = %g;  R = %g;  x = %g;  y = %g\n", r, R, x, y)
   plot()
   x0 = sqrt(R^2 - y^2)
   segment(-x0, y, x0, y)
   θ = atand(y, x0)
   circle(0, 0, R, beginangle=θ, endangle=180 - θ)
   circle2(x, y + r, r, :blue)
   circle(0, y - r, r, :blue)
   segment(0, 0, x*R/(R - r), (y + r)*R/(R - r))
   segment(0, 0, -x*R/(R - r), (y + r)*R/(R - 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(x, y + r, "等円:r,(x,y+r) ", :blue, :right, :vcenter)
       point(0, y - r, "等円:r\n(0,y-r)", :blue, :center, :bottom, delta=2delta)
       point(0, y, "y", :black, :center, :bottom, delta=2delta)
       point(x0, y, "(sqrt(R^2-y^2),y)", :red, :right, delta=-delta)
       point(x - r, y - 1.5r, @sprintf("r = %g\nR = %.15g\nx = %.15g\ny = %.15g", r, R, x, y),
           :black, :left, deltax=5delta, mark=false)
   end
end;

draw(1/2, true)