算額（その526）
和算図形問題あれこれ - 令和4年8月の問題-No.2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html
正方形内に大半円，中半円，小半円，甲円，乙円，丙円が入っている。正方形の一辺の長さが 10 寸のとき，丙円の直径はいかほどか。
大半円の半径と中心座標を r1, (x1, y1); r1 = a/2; x1 = 1; y1 = r1
中半円の半径と中心座標を r2, (x2, y2); x2 = 0; y2 = a - r2
小半円の半径と中心座標を r3, (x3, y3); x3 = r3; y3 = 0
甲円の半径と中心座標を r4, (x4, y4)
乙円の半径と中心座標を r5, (x5, y5)
丙円の半径と中心座標を r6, (x6, y6)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms a::positive,
r1::positive, x1::positive, y1::positive,
r2::positive, x2::positive, y2::positive,
r3::positive, x3::positive, y3::positive,
r4::positive, x4::positive, y4::positive,
r5::positive, x5::positive, y5::positive,
r6::positive, x6::positive, y6::positive;
r1 = a/2; x1 = a; y1 = r1
x2 = 0; y2 = a - r2
x3 = r3; y3 = 0
eq1 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2
eq2 = (x1 - x4)^2 + (y1 - y4)^2 - (r1 + r4)^2
eq3 = (x1 - x5)^2 + (y1 - y5)^2 - (r1 + r5)^2
eq4 = (x1 - x6)^2 + (y1 - y6)^2 - (r1 + r6)^2
eq5 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq6 = (x2 - x4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq7 = (x2 - x5)^2 + (y2 - y5)^2 - (r2 + r5)^2
eq8 = (x2 - x6)^2 + (y2 - y6)^2 - (r2 + r6)^2
eq9 = (x3 - x4)^2 + (y3 - y4)^2 - (r3 + r4)^2
eq10 = (x4 - x5)^2 + (y4 - y5)^2 - (r4 + r5)^2
eq11 = (x5 - x6)^2 + (y5 - y6)^2 - (r5 + r6)^2;
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11], (r2, r3, r4, x4, y4, r5, x5, y5, r6, x6, y6))
2-element Vector{NTuple{11, Sym}}:
(3*a/8, a/3, 2*a*(-10 + 7*sqrt(3))/47, a*(23 - 2*sqrt(3))/47, a*(23 - 2*sqrt(3))/47, 2*a*(18 - 7*sqrt(3))/177, a*(81 - 2*sqrt(3))/177, a*(51/59 - 10*sqrt(3)/59), a/11, 5*a/11, 8*a/11)
(3*a/8, a/3, -20*a/47 + 14*sqrt(3)*a/47, -2*sqrt(3)*a/47 + 23*a/47, -2*sqrt(3)*a/47 + 23*a/47, 2*a*(18 - 7*sqrt(3))/177, a*(81 - 2*sqrt(3))/177, a*(51 - 10*sqrt(3))/59, 2*a*(-10 + 7*sqrt(3))/47, a*(23 - 2*sqrt(3))/47, a*(23/47 - 2*sqrt(3)/47))
2 組の解が得られるが，1 番目のものが適解である。
丙円の半径は a/11 である。
a = 10 のとき，直径は 2*10/11 = 1.8181818181818181 である。
res[1][9] |> println
a/11
a = 10
2a/11
1.8181818181818181
その他のパラメータ
r2 = 3.75; r3 = 3.33333
r4 = 0.903981; x4 = 4.15657; y4 = 4.15657
r5 = 0.663915; x5 = 4.38056; y5 = 5.70839
r6 = 0.909091; x6 = 4.54545; y6 = 7.27273
丙円の直径 = 1.81818
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
a = 10
(r2, r3, r4, x4, y4, r5, x5, y5, r6, x6, y6) = (
3*a/8, a/3,
2*a*(-10 + 7*sqrt(3))/47, a*(23 - 2*sqrt(3))/47, a*(23 - 2*sqrt(3))/47,
2*a*(18 - 7*sqrt(3))/177, a*(81 - 2*sqrt(3))/177, a*(51/59 - 10*sqrt(3)/59),
a/11, 5*a/11, 8*a/11)
r1 = a/2; x1 = a; y1 = r1
x2 = 0; y2 = a - r2
x3 = r3; y3 = 0
@printf("r2 = %g; r3 = %g\nr4 = %g; x4 = %g; y4 = %g\nr5 = %g; x5 = %g; y5 = %g\nr6 = %g; x6 = %g; y6 = %g\n", r2, r3, r4, x4, y4, r5, x5, y5, r6, x6, y6)
@printf("丙円の直径 = %g\n", 2r6)
plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
circle(x1, y1, r1, beginangle=90, endangle=270)
circle(x2, y2, r2, :blue, beginangle=-90, endangle=90)
circle(x3, y3, r3, :magenta, beginangle=0, endangle=180)
circle(x4, y4, r4, :brown)
circle(x5, y5, r5, :green)
circle(x6, y6, r6, :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(x1, y1, "r1,(x1,y1) ", :black, :right, :top, delta=-delta/2)
point(x2, y2, " r2,(x2,y2)", :black, :left, :top, delta=-delta/2)
point(x3, y3, " r3,(x3,y3)", :black, :left, :bottom, delta=delta/2)
point(x4, y4, " r4,(x4,y4)", :black, :left, :top, delta=-delta/2)
point(x5, y5, " r5,(x5,y5)", :black, :left, :top, delta=-delta)
point(x6, y6, " r6,(x6,y6)", :black, :left, :bottom, delta=delta)
point(a, a, "(a,a) ", :black, :right, :bottom, delta=delta/2)
end
end;
