算額（その41）

算額（その41）

大阪府茨木市 井於神社 弘化3年
http://www.wasan.jp/osaka/iyo.html

3辺が 3,4,5の直角三角形に，図のように，大円，中円，小円が入っている。それぞれの円の径を求めよ。

図のように，それぞれの円の中心座標と半径を (0, r1, r1), (x2, r2, r2), (x3, r3, r3), 大円と斜辺の接点の座標を (x1, y1) とする。

using SymPy

@syms r1::positive, r2::positive, r3::positive;
@syms x2::positive, x3::positive, x0::positive, x1::positive, y1::positive, l::positive;

まず，AB = AE = 3 - r1, BC = CD = 4 - r1, AB + BC = 5  から 7 - 2r1 = 5。よって r1 = 1 である。

eq1 = (3 - r1) + (4 - r1) - 5
solve(eq1, r1)[1] |> println

   1

次に，以下の方程式を解く。

r1 = 1;
eq2 = (x3 - x2)^2 + (r2 - r3)^2 - (r2 + r3)^2;
eq3 = (4 - x3)/(4 - x2) - r3/r2;
eq4 = (4 - r1)/(4 - x2) - r1/r2;
eq6 = (x3 - r1)^2 + (r1 - r3)^2 - (r1 + 2r2 + r3)^2;
eq8 = 4y1 - 3(4 - x1);
eq9 = (x1 - r1)^2 + (y1 - r1)^2 - r1^2;

res = solve([eq2, eq3, eq4, eq6, eq8, eq9], (r2, r3, x2, x3, x1, y1))

   2-element Vector{NTuple{6, Sym}}:
    (-179/79 + 62*sqrt(10)/79, -81/79 + 36*sqrt(10)/79, 853/79 - 186*sqrt(10)/79, 559/79 - 108*sqrt(10)/79, 8/5, 9/5)
    (11/9 - 2*sqrt(10)/9, 161/81 - 44*sqrt(10)/81, 1/3 + 2*sqrt(10)/3, -53/27 + 44*sqrt(10)/27, 8/5, 9/5)

2 組の解が得られるが，最初の解は r2 < r3 となり不適切解なので，2番目の解が適切。

name = ("r2", "r3", "x2", "x3", "x1", "y1")
j = 2
for i in 1:6
   println("$(name[i]) = $(res[j][i]) = $(res[j][i].evalf())")
end

   r2 = 11/9 - 2*sqrt(10)/9 = 0.519493853295916
   r3 = 161/81 - 44*sqrt(10)/81 = 0.269873863612238
   x2 = 1/3 + 2*sqrt(10)/3 = 2.44151844011225
   x3 = -53/27 + 44*sqrt(10)/27 = 3.19037840916328
   x1 = 8/5 = 1.60000000000000
   y1 = 9/5 = 1.80000000000000

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1
   (r2, r3, x2, x3, x1, y1) = (11/9 - 2*sqrt(10)/9, 161/81 - 44*sqrt(10)/81, 1/3 + 2*sqrt(10)/3, -53/27 + 44*sqrt(10)/27, 8/5, 9/5)
   println("r1 = $r1;  r2 = $r2;  r3 = $r3;  x2 = $x2;  x3 = $x3")
   plot()
   circle(r1, r1, r1, :red)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :brown)
   plot!([0, 4, 0, 0], [0, 0, 3, 0], color=:black, linewidth=0.25)
   if more
       point(r1, r1, "大円 r1", :red)
       point(x2, r2, "中円 r2", :blue)
       point(x3, r3, "小円 r3", :brown)
       point(0, 3, "  A", :red)
       point(x1, y1, "  B", :red)
       point(4, 0, "  C", :red)
       point(r1, 0, " D=r1", :red, :left, :bottom)
       point(0, r1, " E=r1", :red, :left, :bottom)
       point(x2, 0, "x2", :blue, :left, :bottom)
       point(x3, 0, "x3", :brown, :left, :bottom)
   end
end;