算額（その681）

算額（その681）

山口正義: やまぶき4，第58号，2018/12/06.
千葉県君津市鹿野山 鹿野山神野寺 万延二年(1861)
https://yamabukiwasan.sakura.ne.jp/ymbk58.pdf
山口正義：やまぶき，和算と歴史随想 からリンク
https://yamabukiwasan.sakura.ne.jp/page3.html

短径が 3 寸の楕円の中に，直径 2寸4分の等円が 2 個入っている。楕円の周長はいかほどか。ただし，円積率として 7分8厘5毛4糸を使うものとする。

楕円の長半径，短半径，中心座標を a, b, (0, 0)
等円の半径と，中心座標を r, (r, 0)
楕円と等円の接点の座標を (x1, y1)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms a::positive, b::positive, r::positive, x1::positive, y1::positive

eq1 = (x1 - r)^2 + y1^2 - r^2
eq2 = -b^2*x1/(a^2*y1) - (r - x1)/y1
eq3 = x1^2/a^2 + y1^2/b^2 - 1
res = solve([eq1, eq2, eq3], (x1, y1, a))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (b^2/r, b*sqrt(-b^2 + 2*r^2)/r, -b^2*sqrt(1/(b^2 - r^2)))
    (b^2/r, b*sqrt(-b^2 + 2*r^2)/r, b^2*sqrt(1/(b^2 - r^2)))

2 組の解が得られるが，2 番目のものが適解である（最初のものは x 軸に関して線対称な解である）。

x1, y2, a はそれぞれ短半径 b の b/r 倍，sqrt(2r^2 - b^2) 倍，b/sqrt(b^2 - r^2) 倍である。

短径（差渡し径）が 3 寸（短半径が 1.5 寸），等円の直径が 2.4 寸のとき，長半径は 2.5 寸である（差渡し径 5 寸）。

   b = 1.5;  r = 1.2;  x1 = 1.875;  y1 = 0.992157;  a = 2.5

楕円の周長は第二種完全楕円積分による正確な値として 12.763499431699064 が得られる。

問において，「円周率は 0.7854 を使う」とあるのは，円積率のことで円周率になおすと 3.1416 である。当時，円積率としては 0.79（円周率は 3.16）が使われることが多かったであろうにしては，正確すぎる。
しかし，この円積率を使ったとしても，「答」にある 12.7634994317113 有奇を導いた近似式などが不明（術を読み解く力がない）。有効数字 11 桁が一致している。驚異的。

長径と短径の比がきれいな値なので，事前に計算された数表などがあったのかもしれない。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (b, r) = (3//2, 24/20)
   (x1, y1, a) = b .* (b/r, sqrt(2r^2 - b^2)/r, b/sqrt(b^2 - r^2))
   @printf("b = %g;  r = %g;  x1 = %g;  y1 = %g;  a = %g\n", b, r, x1, y1, a)
   plot()
   ellipse(0, 0, a, b, color=:blue)
   circle(r, 0, r)
   circle(-r, 0, r)
   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(x1, y1, " (x1,y1)", :black, :left, :bottom)
       point(r, 0, "r", :red, :center, :bottom, delta=delta)
       point(0, b, " b", :blue, :left, :bottom, delta=delta)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta)
   end
end;