算額（その1498）

算額（その1498）

二十五 一関市萩荘 萩荘春日神社 弘化4年(1847)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

今有如図 3032
https://w.atwiki.jp/sangaku/pages/192.html

意訳：長方形の中に，等円 2 個を容れる。長方形の対角線を 2 本引き，その延長線と長方形の上辺に 3 箇所で接する楕円を置く。楕円の短径と等円の直径が与えられたとき，楕円の長径を求める術を述べよ。

山村の図では何がなんだかわからないが，「今有如図」で解決策が見えてくる。

算法助術の公式97 を用いる。
公式の記述中では変数は大文字で書き（eq97），変数が取る実際の値（定数とは限らない）は小文字で代入する（eq)。
その後 eq を方程式として解く。eq で大文字のものが求めるものである（今の場合は P すなわち楕円の長軸）。
公式を用いるに当たり，図に補助線として長方形内の対角線の延長と楕円に接する水平線を描く。対角線の延長と水平線は二等辺三角形を作り，楕円は二等辺三角形に内接する。これが公式97が対象とする図形である。
この二等辺三角形は長方形内の二等辺三角形と相似で，相似比は r:(2b + r) である。

include("julia-source.txt");

using SymPy

@syms A::positive, B::positive, C::positive, H::positive,
      P::positive, Q::positive
eq97 = -(B^2 - C^2)^2*H*Q^2 + (B^2 - C^2)^2*Q^3 + A^4*H*(2H - Q)^2 - A^4*Q*(2H - Q)^2 - A^2*H*P^2*(2H -  Q)^2
eq97 |> println

    A^4*H*(2*H - Q)^2 - A^4*Q*(2*H - Q)^2 - A^2*H*P^2*(2*H - Q)^2 - H*Q^2*(B^2 - C^2)^2 + Q^3*(B^2 - C^2)^2

@syms b, r
eq = eq97(A => 4(2b + r), B => sqrt((2b + r)^2 + 4(2b + r)^2), C => sqrt((2b + r)^2 + 4(2b + r)^2), H => 2b + r, Q => 2b)
eq = eq |> simplify
eq |> println

    64*(b + r)^2*(2*b + r)^3*(-P^2 - 32*b*(2*b + r) + 16*(2*b + r)^2)

# P:楕円の長軸
ans_P = solve(eq, P)[2]
ans_P |> println

    4*sqrt(r*(2*b + r))

等円の直径が 2，楕円の短径が 3 のとき，楕円の長径は

ans_P(r => 2/2, b => 3/2).evalf() |> println

    8.00000000000000

術は「長径 = sqrt(等円径*(2短径 + 等円径))」となっているが，術のとおりだと 4 になってしまう(山村はさらに間違えている）。
これを 2 倍しないと正しくない。

短径 = 3
等円径 = 2
sqrt(等円径*(2短径 + 等円径))

    4.0

function draw(b, r, more=false)  # b:楕円の短半径と，r:円の半径
    pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    a = 4*sqrt(r*(2*b + r))/2  # 長半径
    @printf("短径 = %g，等円径 = %g のとき，長径は %g である。\n", 2b, 2r, 2a)
    plot(r.*[2, 2, -2, -2, 2], r.*[-1, 1, 1, -1, -1], color=:blue, lw=0.5)
    plot!([-2r, 2(2b+r), -2(2b + r), 2r], [-r, 2b + r, 2b + r, -r], color=:blue, lw=0.5)
    circle2(r, 0, r)
    ellipse(0, r + b, a, b, color=:green)
    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)
    end
end;

draw(3/2, 2/2, true)