13 岩手県奥州市（旧水沢市佐倉河） 胆沢城跡鎮守府八幡宮 弘化2年(1845)
安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：円4個，外円，菱形
#Julia #SymPy #算額 #和算 #数学
大円内に互いに内接・外接する小円 3 個と，菱形を容れる。小円の直径が与えられたとき，黒積を求めよ。
外円の半径と中心座標を \(R,\ (0,\ 0)\)
小円の半径と中心座標を \(r,\ (0,\ R - r),\ (x,\ y)\)
∠AOB を \(\theta\)
とおき，以下の連立方程式を解く。
一度に解けないので，まず eq1, eq3 から \(x,\ y\) を求め，結果を eq2 に代入して \(b\) を求める。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms r, R, x, y, a, b
R = b + r
a = sqrt(R^2 - (b - R)^2)
eq1 = x^2 + y^2 - (R - r)^2
eq2 = dist2(0, 2b - R, a, b - R, x, y, r)
eq3 = x^2 + (R - r - y)^2 - 4r^2;
#res = solve([eq1, eq2, eq3], (x, y, b))
res = solve([eq1, eq3], (x, y))[2]
(2*r*sqrt(-(-b + r)*(b + r))/b, -(-b^2 + 2*r^2)/b)
eq12 = eq2(x => res[1], y => res[2]) |> simplify |> numerator |> x -> x/r^2 |> factor
eq12 |> println
(-b + 2*r)*(-3*b^2 - 2*b*r + 4*r^2 - 4*sqrt(b^2 - r^2)*sqrt(b^2 + 2*b*r))
eq12 は \( (2r - b)\) という因数を持つので，\(b = 2r\) はその解である。また，適解であることが確かめられる。
よって，\(x,\ y\) は，\(r = 1/2,\ b = 2r = 1\) を代入することにより，\(\sqrt{3}/2,\ 1/2\) である。
# x
res[1](r => 1//2, b => 1) |> println
sqrt(3)/2
# y
res[2](r => 1//2, b => 1) |> println
1/2
結果としてこの場合には，小円の中心は一辺の長さが \(2r\) の正六角形の頂点である。
黒積は「扇形AOBの面積 - 扇型ACDの面積 - ⊿OCDの面積」の 4 倍である。
∠ACD = 120°，∠AOB = 30° なので，面積 \(S\) は \(4(\pi R^2(30/360) - (R - r)\cos(30°)r/2 - \pi r^2(120/360))\) である。
簡約化すると \(-\sqrt{3} b r - 4 \pi r^2/3 + \pi(b + r)^2/3\) である。
S = 4*(PI*R^2*(30//360) - (R - r)*cosd(Sym(30))*r//2 - PI*r^2*(120//360)) |> simplify
S |> println
-sqrt(3)*b*r - 4*pi*r^2/3 + pi*(b + r)^2/3
術は，小円の直径が 1 のときで，そのとき \(b\) も 1 になるので，黒積の数値は \(5\pi/12 - \sqrt{3}/2 = 0.4429715352113086\)
S(b =>1, r => 1//2) |> simplify |> println
-sqrt(3)/2 + 5*pi/12
-sqrt(3)/2 + 5*pi/12
0.4429715352113086
術は円積率=0.76 を使っているので，黒積は「\( (10\cdot 0.76 - \sqrt{27})/6\)」で，計算結果は 0.40064126288222796 となる。
円周率を使うと \( (10\pi/4 - \sqrt{27})/6 = 0.44297153521130844\) である。
\( (10\pi/4 - \sqrt{27})/6\) は簡約化すると \(-\sqrt{3}/2 + 5\pi/12\) となり，上述の式に一致する。
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r = 1/2
b = 2r
x = 2r*sqrt( (b - r)*(b + r))/b
y = (b^2 - 2r^2)/b
R = b + r
a = sqrt(R^2 - (b - R)^2)
θ = asind(r/(R - r))
S = -sqrt(3)*b*r - 4*pi*r^2/3 + pi*(b + r)^2/3
@printf("小円の直径が %g のとき，黒積は %g である。\n", 2r, 4S)
@printf("r = %g; x = %g; y = %g; b = %g; a = %g; R = %g; θ = %g\n", r, x, y, b, a, R, θ)
plot([0, a, 0, -a, 0], [-R, b - R, 2b - R, b - R, -R], color=:green, lw=0.5)
circle(0, 0, R)
circle(0, R - r, r, :blue)
circle2(x, y, r, :blue)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
point(R, 0, " R", :red, :left, :vcenter)
point(x, y, "小円:r,(x,y)", :blue, :center, delta=-delta/2)
θ2 = 30:0.1:90
ax = R.*cosd.(θ2)
ay = R.*sind.(θ2)
θ3 = 90:-0.1:-30
append!(ax, r.*cosd.(θ3))
append!(ay, (R - r) .+ r.*sind.(θ3))
θ4 = 150:-0.1:30
append!(ax, x .+ r.*cosd.(θ4))
append!(ay, y .+ r.*sind.(θ4))
plot!(ax, ay, seriestype=:shape, color=:gray80, fillcolor=:gray80, lw=0.5)
plot!(-ax, ay, seriestype=:shape, color=:gray80, fillcolor=:gray80, lw=0.5)
segment(0, 0, R*cosd(60), R*sind(60), :gray60)
point(0, R, "A", :red, :center, :bottom, delta=delta/2)
point(R*cosd(60), R*sind(60), "B", :red, :left, :bottom, delta=delta/2)
point(0, R - r, "C ", :blue, :right, :vcenter)
point(x/2, (R - r + y)/2, " D", :blue, :left)
segment(0, R - r, x/2, (R - r + y)/2, :gray60)
segment(0,0, 0, R - r, :gray60)
point(0, 0, "O", :red, :center, delta=-delta/2)
point(0, 2b - R, "2b-r ", :green, :right, :vcenter)
point(0, b - R, "b-r ", :green, :right, :vcenter)
point(a, b - R, "(a,b-r) ", :green, :right, :vcenter)
end
end;
draw(true)
以下のアイコンをクリックして応援してください
