算額(その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;