算額（その503）

算額（その503）

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/8
https://dl.ndl.go.jp/pid/3508431/1/9

長方形内に2本の斜線を引き，分割された領域に甲円，乙円を入れる。
長方形の長辺と短辺の長さが与えられたとき，乙円の直径を求めよ。

長方形の長辺と短辺の長さを 2a, 2r1 とする。
甲円の半径と中心座標を r1, (a - r1, 0)
乙円の半径と中心座標を r2, (0, r1 - r2)
斜線と長方形の上辺の交点を (b, r1)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a, b, r1, r2
eq1 = b/sqrt(b^2 + r1^2) - r2/(r1 - r2)
eq2 = r1/(a - r1) - r1/sqrt(b^2 + r1^2);
eq2 = (a - r1) - sqrt(b^2 + r1^2);
res = solve([eq1, eq2], (r2, b))

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

2 組の解が得られるが，2 番目のものが適解である。
SymPy では式を簡約化できないが，容易に以下のようにできる。

b = sqrt(a^2 - 2a*r1)
r2 = (r1*b)/(b + a - r1)

長方形の長辺と短辺が 150, 72 のとき，乙円の半径は 10，直径は 20 である。

長 = 150;  平 = 72;  a = 75;  r1 = 36;  b = 15;  r2 = 10
乙円の直径 = 20

using Plots

function draw(長, 平, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r1) = (長, 平) ./ 2
   # (r2, b) = (r1*sqrt(a*(a - 2*r1))/(sqrt(a*(a - 2*r1)) + sqrt(a^2 - 2*a*r1 + r1^2)), sqrt(a*(a - 2*r1)))
   b = sqrt(a^2 - 2a*r1)
   r2 = (r1*b)/(b + a - r1)
   @printf("長 = %g;  平 = %g;  a = %g;  r1 = %g;  b = %g;  r2 = %g\n", 2a, 2r1, a, r1, b, r2)
   @printf("乙円の直径 = %g\n", 2r2)
   plot([a, a, -a, -a, a], [-r1, r1, r1, -r1, -r1], color=:black, lw=0.5)
   circle(a - r1, 0, r1)
   circle(r1 - a, 0, r1)
   circle(0, r1 - r2, r2, :blue)
   circle(0, r2 - r1, r2, :blue)
   segment(b, r1, -b, -r1, :green)
   segment(-b, r1, b, -r1, :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(a, 0, "a ", :red, :right, :bottom, delta=delta)
       point(0, r1 - r2, " 乙円:r2,(0,r1-r2)", :blue, :left, :vcenter)
       point(a - r1, 0, "甲円:r1,(a-r1,0)", :red, :center, :bottom, delta=2delta)
       point(b, r1, "(b,r1)", :green, :center, :bottom, delta=2delta)
   else
       plot!(showaxis=false)
   end
end;