算額（その231）
福島県郡山市西田町丹伊田字宮作239番地 鹿島大神宮 令和3年 武田健司さんが奉納
https://kashimadaijingu.jp/archives/2785.html
半径 1 の 2 つの円が，互いの中心を通り重なり合っている。その中に 27 個の円が図のように接している。2 つの円が作る領域から，27個の円を除いた領域の面積を求めよ。
交差する 2 つの円の半径を 1 とする。
x軸上に並ぶ 3 個の円（円1）の半径を r1 とする。r1 = 1/2 である。中心座標は左から順に (-2r1, 0), (0, 0), (2r1, 0)
ついで半径の大きい順に 半径，中心座標を以下のように定める。
円2: r2, (x2, y2),
円3: r3, (0, r1 + r3),
円4: r4, (x4, y4),
円5: r5, (x5, y5),
円6: r6, (x6, y6),
円7: r7, (x7, y7),
円8: r8, (0, r1 + 2r3 + r8)
以下の方程式を解く。
using SymPy
@syms r1::positive,
r2::positive, x2::positive, y2::positive,
r3::positive,
r4::positive, x4::positive, y4::positive,
r5::positive, x5::positive, y5::positive,
r6::positive, x6::positive, y6::positive,
r7::positive, x7::positive, y7::positive,
r8::positive;
r1 = 1//2
eq01 = r1^2 + (r1 + 2r3 + r8)^2 - (2r1 - r8)^2 # 右外円に円8が内接
eq02 = r1^2 + (r1 + r3)^2 - (2r1 - r3)^2 # 右外円に円3が内接
eq03 = x6^2 + (r1 + r3 - y6)^2 - (r3 + r6)^2 # 円3と円6が外接
eq04 = x6^2 + y6^2 - (r1 + r6)^2 # 中円1と円6が外接
eq05 = (x6 + r1)^2 + y6^2 - (2r1 - r6)^2 # 左外円に円6が内接
eq06 = (x2 - r1)^2 + y2^2 - (2r1 - r2)^2 # 右外円に円2が内接
eq07 = (x2 + r1)^2 + y2^2 - (2r1 + r2)^2 # 左外円と円2が外接
eq08 = (2r1 - x2)^2 + y2^2 - (r1 + r2)^2 # 右円1と円2が外接
eq09 = (x4 - r1)^2 + y4^2 - (2r1 - r4)^2 # 右外円に円4が内接
eq10 = (x4 -2r1)^2 + y4^2 - (r1 + r4)^2 # 右円1と円4が外接
eq11 = (x4 - x2)^2 + (y2 - y4)^2 - (r2 + r4)^2 # 円2と円4が外接
eq12 = (x5 - r1)^2 + y5^2 - (2r1 - r5)^2 # 右外円に円5が内接
eq13 = (x2 - x5)^2 + (y5 - y2)^2 - (r2 + r5)^2 # 円2と円5が外接
eq14 = (x5 + r1)^2 + y5^2 - (2r1 + r5)^2 # 左外円と円5が外接
eq15 = (x7 - r1)^2 + y7^2 - (2r1 - r7)^2 # 右外円に円7が内接
eq16 = (x7 + r1)^2 + y7^2 - (2r1 + r7)^2 # 左外円と円7が外接
eq17 = (x5 - x7)^2 + (y7 - y5)^2 - (r5 + r7)^2; # 円5と円7が外接
res = solve([eq01, eq02, eq03, eq04, eq05, eq06, eq07, eq08, eq09, eq10, eq11, eq12, eq13, eq14, eq15, eq16, eq17],
(r2, x2, y2, r3, r4, x4, y4, r5, x5, y5, r6, x6, y6, r7, x7, y7, r8))
4-element Vector{NTuple{17, Sym}}:
(3/10, 3/5, 2*sqrt(3)/5, 1/6, 39/121 - 12*sqrt(3)/121, 36*sqrt(3)/121 + 129/242, 30/121 + 28*sqrt(3)/121, 9/82, 9/41, 20*sqrt(3)/41, 1/11, 5/22, 6/11, 27/730, 27/365, 182*sqrt(3)/365, 1/66)
(3/10, 3/5, 2*sqrt(3)/5, 1/6, 39/121 - 12*sqrt(3)/121, 36*sqrt(3)/121 + 129/242, 30/121 + 28*sqrt(3)/121, 9/82, 9/41, 20*sqrt(3)/41, 1/11, 5/22, 6/11, 3/10, 3/5, 2*sqrt(3)/5, 1/66)
(3/10, 3/5, 2*sqrt(3)/5, 1/6, 12*sqrt(3)/121 + 39/121, 129/242 - 36*sqrt(3)/121, -30/121 + 28*sqrt(3)/121, 9/82, 9/41, 20*sqrt(3)/41, 1/11, 5/22, 6/11, 27/730, 27/365, 182*sqrt(3)/365, 1/66)
(3/10, 3/5, 2*sqrt(3)/5, 1/6, 12*sqrt(3)/121 + 39/121, 129/242 - 36*sqrt(3)/121, -30/121 + 28*sqrt(3)/121, 9/82, 9/41, 20*sqrt(3)/41, 1/11, 5/22, 6/11, 3/10, 3/5, 2*sqrt(3)/5, 1/66)
1 番目の解が妥当なものである。
names = ("r2", "x2", "y2", "r3", "r4", "x4", "y4", "r5", "x5", "y5", "r6", "x6", "y6", "r7", "x7", "y7", "r8")
for i = 1:17
println("$(names[i]) = $(res[1][i]) = $(res[1][i].evalf())")
end
r2 = 3/10 = 0.300000000000000
x2 = 3/5 = 0.600000000000000
y2 = 2*sqrt(3)/5 = 0.692820323027551
r3 = 1/6 = 0.166666666666667
r4 = 39/121 - 12*sqrt(3)/121 = 0.150540415778293
x4 = 36*sqrt(3)/121 + 129/242 = 1.04837875266512
y4 = 30/121 + 28*sqrt(3)/121 = 0.648739029850649
r5 = 9/82 = 0.109756097560976
x5 = 9/41 = 0.219512195121951
y5 = 20*sqrt(3)/41 = 0.844902832960428
r6 = 1/11 = 0.0909090909090909
x6 = 5/22 = 0.227272727272727
y6 = 6/11 = 0.545454545454545
r7 = 27/730 = 0.0369863013698630
x7 = 27/365 = 0.0739726027397260
y7 = 182*sqrt(3)/365 = 0.863652731445303
r8 = 1/66 = 0.0151515151515152
円1〜円8 の円の面積（ひとつあたり）
s = PI .* [r1, res[1][1], res[1][4], res[1][5], res[1][8], res[1][11], res[1][14], res[1][17]].^2
s |> println
Sym[pi/4, 9*pi/100, pi/36, pi*(39/121 - 12*sqrt(3)/121)^2, 81*pi/6724, pi/121, 729*pi/532900, pi/4356]
図中にある個数
num = [3, 4, 2, 4, 4, 4, 4, 2]
8-element Vector{Int64}:
3
4
2
4
4
4
4
2
図中にある全ての円の面積の和
S1 = sum(s .* num).evalf()
S1 |> println
4.22035198660619
二円の構成する面積。第1象限の面積の 4 倍
@syms x
S2 = integrate(sqrt(4r1^2 - (x - r1)^2), (x, 0, 3r1)).evalf() * 4
S2 |> println
5.05481560857083
二円の構成する面積から全ての円の面積の和を引いたもの
S2 - S1 |> println
0.834463621964642
using Plots
using Printf
function circle4(x, y, r, color=:red; fill=false)
circle(x, y, r, color, fill=fill)
circle(x, -y, r, color, fill=fill)
circle(-x, y, r, color, fill=fill)
circle(-x, -y, r, color, fill=fill)
end;
function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
θ = beginangle:0.1:endangle
x = r.*cosd.(θ)
y = r.*sind.(θ)
if fill
plot!(ox .+ x, oy .+ y, linecolor=:black, linewidth=0.5, seriestype=:shape, fillcolor=:black)
else
plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end
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; fill=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 1/2
(r2, x2, y2, r3, r4, x4, y4, r5, x5, y5, r6, x6, y6, r7, x7, y7, r8) = res[1]
@printf("r1 = %.7f\n", r1)
@printf("r2 = %.7f; x2 = %.7f; y2 = %.7f\n", r2, x2, y2)
@printf("r3 = %.7f\n", r3)
@printf("r4 = %.7f; x4 = %.7f; y4 = %.7f\n", r4, x4, y4)
@printf("r5 = %.7f; x5 = %.7f; y5 = %.7f\n", r5, x5, y5)
@printf("r6 = %.7f; x6 = %.7f; y6 = %.7f\n", r6, x6, y6)
@printf("r7 = %.7f; x7 = %.7f; y7 = %.7f\n", r7, x7, y7)
@printf("r8 = %.7f\n", r8)
plot()
circle(r1, 0, 2r1, :red)
circle(-r1, 0, 2r1, :red)
circle(0, 0, r1, :green, fill=fill)
circle(2r1, 0, r1, :green, fill=fill)
circle(-2r1, 0, r1, :green, fill=fill)
circle(0, r1 + r3, r3, fill=fill)
circle(0, -r1 - r3, r3, fill=fill)
circle(0, r1 + 2r3 + r8, r8, :black, fill=fill)
circle(0, -r1 - 2r3 - r8, r8, :black, fill=fill)
circle4(x2, y2, r2, :gold, fill=fill)
circle4(x4, y4, r4, :magenta, fill=fill)
circle4(x5, y5, r5, :blue, fill=fill)
circle4(x6, y6, r6, :darkorange3, fill=fill)
circle4(x7, y7, r7, :cyan4, fill=fill)
if more == true
point(0, 0, "0 ", :green, :right)
point(r1, 0, "r1 ", :green, :right)
point(2r1, 0, "2r1 ", :green, :right)
point(3r1, 0, "3r1 ", :green, :right)
point(2r1, r1/3, "円1", :green, mark=false)
point(x2, y2, "円2", :gold, mark=false)
point(0, r1 + r3, "円3", :red, mark=false)
point(x4, y4, "円4", :magenta, mark=false)
point(x5, y5, "円5", :blue, mark=false)
point(x6, y6, "円6", :darkorange3, mark=false)
point(x7, y7, "円7", :cyan4, mark=false)
point(0, r1 + 2r3 + 7r8, "円8", :black, :center, :bottom, mark=false)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
