算額（その516）

算額（その516）

『算法新書』 文政13年(1830)
一関市博物館>>和算に挑戦>>平成26年度出題問題>>平成26年度出題問題(2)[中級問題]（中学・高校生向き）
https://www.city.ichinoseki.iwate.jp/museum/wasan/h26/normal.html

外円の中に弦 2 本と小円が 3 個入っている。左右の小円は弦にその中点で接している。外円の直径が 1 寸のとき，小円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
小円の半径と中心座標を r, (x, y), (0, r - R)
弦と円の右下の交点座標を (a, -sqrt(R^2 - a^2))
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R, r, a, x, y

eq1 = r/(2R - r) - (R - 2r)/R;

eq1 だけで，R が与えられたときに r を求めることができる。

solve(eq1, r) |> println

   Sym[R*(3 - sqrt(5))/2, R*(sqrt(5) + 3)/2]

r = R*(3 - sqrt(5))/2 が適解である。
すなわち，外円の直径が 1 のとき，小円の直径は (3 - sqrt(5))/2 ≒ 0.3819660112501051 である。

(3 - sqrt(5))/2

   0.3819660112501051

図を描くために必要な a, x, y を求めるには以下の eq2, eq3, eq4 を加えて，4元連立方程式を解けばよい。

eq2 = 2sqrt(R^2 - (R - 2r)^2) - sqrt((R + sqrt(R^2 - a^2))^2 + a^2)
eq3 = x^2 + y^2 - (R - r)^2
eq4 = (y/x)*((R + sqrt(R^2 - a^2))/a) + 1
res = solve([eq1, eq2, eq3, eq4], (r, a, x, y));

16 組の解が得られるが，6 番目のものが適解である。

res[6][1] |> println  # r

   R*(3 - sqrt(5))/2

res[6][2] |> println  # a

   4*R*sqrt(-38 + 17*sqrt(5))

res[6][3] |> println  # x

   sqrt(2)*sqrt(R^2*(-8*R*(-199 + 89*sqrt(5))/(R + sqrt(609 - 272*sqrt(5))*sqrt(R^2)) - sqrt(5) + 3))/2

res[6][4] |> println  # y

   2*sqrt(R^3/(R + sqrt(609 - 272*sqrt(5))*sqrt(R^2)))*sqrt(-199 + 89*sqrt(5))

res[6]

   (R*(3 - sqrt(5))/2, 4*R*sqrt(-38 + 17*sqrt(5)), sqrt(2)*sqrt(R^2*(-8*R*(-199 + 89*sqrt(5))/(R + sqrt(609 - 272*sqrt(5))*sqrt(R^2)) - sqrt(5) + 3))/2, 2*sqrt(R^3/(R + sqrt(609 - 272*sqrt(5))*sqrt(R^2)))*sqrt(-199 + 89*sqrt(5)))

   R = 1;  r = 0.381966;  a = 0.458792;  x = 0.600566;  y = 0.145898

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 1
   (r, a, x, y) = (
       R*(3 - sqrt(5))/2,
       4*R*sqrt(-38 + 17*sqrt(5)),
       sqrt(2)*sqrt(R^2*(-8*R*(-199 + 89*sqrt(5))/(R + sqrt(609 - 272*sqrt(5))*sqrt(R^2)) - sqrt(5) + 3))/2,
       2*sqrt(R^3/(R + sqrt(609 - 272*sqrt(5))*sqrt(R^2)))*sqrt(-199 + 89*sqrt(5)))
   @printf("R = %g;  r = %g;  a = %g;  x = %g;  y = %g\n", R, r, a, x, y)
   plot()
   circle(0, 0, R)
   circle(x, y, r, :blue)
   circle(-x, y, r, :blue)
   circle(0, r - R, r, :blue)
   segment(0, R, a, -sqrt(R^2 - a^2), :green)
   segment(0, R, -a, -sqrt(R^2 - a^2), :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, y, " (x,y)", :blue, :left, :vcenter)
       point(a, -√(R^2 - a^2), "(a, -√(R^2-a^2))", :blue)
       point(0, r - R, " r-R", :blue, :left, :vcenter)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta)
   end
end;