算額（その494）

算額（その494）

北越長岡蔵王権現社（金峰神社） 寛政7年(1795)
涌田和芳，山田章，島津大，若山想思，市野梨保子，矢野敦大，山崎瑠聖: 長岡金峰神社の紛失算額，長岡工業高等専門学校研究紀要，第 55 巻，p.39-44，2019
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_55/55_39wakuta.pdf

長径，短径が a, b の楕円内に半径 r の円が 3 個入っている。左右の円は楕円と 2 点で接し中央の円と外接している。中央の円は楕円に内接している。

長径，短径が与えられたとき，円の半径を求めよ。

注: 和算では楕円の長径，短径は今でいう長径，短径の 2 倍である。同じく，円の径も直径を意味する。
ここでは他の記事と同じで，長径，短径は今の意味，円も半径で扱う。

eq2 は，山本賀前著『算法助術』の84番目の公式である。
山本賀前: 算法助術，天保12年(1841)，東北大学DB.
インターネットでのダウンロード・閲覧は
中村信弥編著:和算の図形公式-『算法助術』-，2008.
http://www.wasan.jp/kosiki/kosiki3.pdf

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r::positive, x::positive;

eq1 = x^2 + (b - r)^2 - 4r^2
eq2 = (b^2 - r^2)*(a^2 - b^2)/b^2 - x^2
res = solve([eq1, eq2], (r, x))

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

円の半径は以下のように簡約化される。これは，術で述べられている式と同じである。

res[1][1] |> simplify |> println

   a^2*b/(a^2 + 2*b^2)

右側の円の中心座標（x座標）も以下のように簡約化される。

res[1][2] |> simplify |> println

   2*b*sqrt(a^4 - b^4)/(a^2 + 2*b^2)

長径と短径が 18, 4.5 のとき，円の半径は 4 である。

(a, b) = (36, 9) .// 2
a^2*b/(a^2 + 2*b^2) |> println

   4//1

   r = 4;  x = 7.98436;  a = 18;  b = 4.5

算額の図は，長径が短径の 2 倍程度なので上の図のような感じであるが，長径と短径が 18, 4.5 のときは，かなり横長なものになる。

using Plots

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x) = (a^2*b/(a^2 + 2b^2), 2b*sqrt(a^4 - b^4)/(a^2 + 2b^2))
   @printf("r = %g;  x = %g;  a = %g;  b = %g\n", r, x, a, b)
   plot()
   circle(0, b - r, r, :red)
   circle(x, 0, r, :red)
   circle(-x, 0, r, :red)
   ellipse(0, 0, a, b, color=:green)
   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(x, 0, "x", :red, :center, :bottom, delta=delta)
       point(a, 0, " a", :green, :left, :bottom, delta=delta)
       point(0, b - r, " b-r", :red, :left, :vcenter)
       point(0, b, " b", :green, :left, :bottom, delta=delta)
   else
       plot!(showaxis=false)
   end
end;

アルベロス

永井信一:「アルベロスに関連した問題」
http://www2.ttcn.ne.jp/~nagai/waseda/wasan/arbe.pdf

図のように，半円内に 2 　つの半円と円が入っている。
外半円，大半円，小半円の半径 r0, r1, r2 がそれぞれ 5, 3, 2 であるとき，円の半径 r を求めよ。

円の中心座標を x, y とし，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive, r::positive,
     x::positive, y::positive;

eq1 = (x + r0 - r1)^2 + y^2 - (r1 + r)^2
eq2 = (r0 - r2 - x)^2 + y^2 - (r2 + r)^2
eq3 = x^2 + y^2 - (r0 - r)^2
res = solve([eq1, eq2, eq3], (r, x, y))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (r0*(r0 - r1)*(r0 - r2)/(r0^2 - r1*r2), r0^2*(r1 - r2)/(r0^2 - r1*r2), -2*r0*sqrt(r1)*sqrt(r2)*sqrt((r0 - r1)*(r0 - r2))/(r0^2 - r1*r2))
    (r0*(r0 - r1)*(r0 - r2)/(r0^2 - r1*r2), r0^2*(r1 - r2)/(r0^2 - r1*r2), 2*r0*sqrt(r1)*sqrt(r2)*sqrt((r0 - r1)*(r0 - r2))/(r0^2 - r1*r2))

2番目のものが適解である。

円の半径は以下の式で求まる。

r0*(r0 - r1)*(r0 - r2)/(r0^2 - r1*r2)

外半円，大半円，小半円の半径 r0, r1, r2 がそれぞれ 5, 3, 2 であるとき，円の半径は 30/19 である。

res[2][1](r0 => 5, r1 => 3, r2 => 2) |> println

   30/19

using Plots

function draw(r1, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = r1 + r2
   (r, x, y) = (r0*(r0 - r1)*(r0 - r2)/(r0^2 - r1*r2), r0^2*(r1 - r2)/(r0^2 - r1*r2), 2*r0*sqrt(r1)*sqrt(r2)*sqrt((r0 - r1)*(r0 - r2))/(r0^2 - r1*r2))
   @printf("r = %g;  x = %g;  y = %g\n", r, x, y)
   plot()
   circle(0, 0, r0, :red, beginangle=0, endangle=180)
   circle(r1 - r0, 0, r1, :blue, beginangle=0, endangle=180)
   circle(r0 - r2, 0, r2, :magenta, beginangle=0, endangle=180)
   circle(x, y, r, :orange)
   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(r0, 0, "r0", :red, :left, :bottom, delta=delta)
       point(r1 - r0, 0, "r1-r0", :blue, :center, :bottom, delta=delta)
       point(r0 - r2, 0, "r0-r2", :magenta, :center, :bottom, delta=delta)
       point(x, y, " r,(x,y)", :orange, :left, :vcenter)
   else
       plot!(showaxis=false)
   end
end;