算額(その458)
長野県長野市 久保寺観音堂 享和3年(1803)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
円に正三角形が内接している。
正三角形の一辺と外積(円と正三角形の間の面積)の和を A,矢(弦の中点と弧の中点を結ぶ線分の長さ)と円周の和を B とするとき,円の面積を求めよ。
最初に,簡単な解を求める方法について書く。条件 A は不要で,かつ,条件 B はとてつもなく簡単なものである。
図において,「矢」は直径の 1/4 である。
したがって,円の半径を r とすれば,「矢と円周の和を B とする」は 「r/2 + 2r*π = B」 である。
半径は,この方程式を r について解けば得られる。
include("julia-source.txt")
using SymPy
@syms r, B
solve(r/2 + 2r*PI - B, r)[1] |> println
2*B/(1 + 4*pi)
円の半径が 2B/(1 + 4π) なので,面積は π*(2B/(1 + 4π))^2 である。
B = 8.139822368615503
r = 2B/(1 + 4π)
(r, π*(2B/(1 + 4π))^2)
(1.2, 4.523893421169302)
まともにやるととてつもなく面倒なことになる。
円の半径を r,正三角形の一辺の長さを 2a とする。
円,正三角形の面積をそれぞれ s1, s2 とおく。
include("julia-source.txt")
using SymPy
@syms r::positive, a::positive, s1::positive, s2::positive, 矢::positive, A::positive, B::positive
s1 = PI*r^2
a = r*sqrt(Sym(3))/2
s2 = a^2*sqrt(Sym(3))
矢 = r - r/2;
A, B を定義する。
eq1 = 2a + (s1 - s2) - A
eq2 = 矢 + PI*2r - B;
差を取って,r について解くのが常套手段ではあるが,eq2 を見るとこれを r について解けばよいことがわかる。
solve(eq2, r) |> println
Sym[2*B/(1 + 4*pi)]
ということで,冒頭の簡単な方法へ戻る。
なお,eq1 は不要であるし,差を取ったのでは条件 A が残り,面倒な式になる。無理にやると解は二通り得られるが,どちらが適解かは単純にはいえない。以下に示すように,どちらを取るかは r に依存するのである。これでは自分の尻尾を飲み込もうとしている蛇と同じである。
eq3 = eq1 - eq2
res = solve(eq3, r);
res[1] |> println
(-sqrt(-12*sqrt(3)*A + 16*pi*A - 16*pi*B + 12*sqrt(3)*B - 16*sqrt(3)*pi - 4*sqrt(3) + 13 + 8*pi + 16*pi^2) - 4*pi - 1 + 2*sqrt(3))/(-4*pi + 3*sqrt(3))
res[2] |> println
(sqrt(-12*sqrt(3)*A + 16*pi*A - 16*pi*B + 12*sqrt(3)*B - 16*sqrt(3)*pi - 4*sqrt(3) + 13 + 8*pi + 16*pi^2) - 4*pi - 1 + 2*sqrt(3))/(-4*pi + 3*sqrt(3))
円の半径を与えて,A, B を計算し,それを eq3 の解 res により円の半径を再計算する関数 rev を書いてみる。
r > (-4*pi - 1 + 2*sqrt(3))/(-4*pi + 3*sqrt(3) の場合には最初のもの,そうでない場合には2番めのものが適解である。
function rev(r) # 半径を与えて検算する
s1 = pi*r^2
a = r*sqrt(3)/2
s2 = a^2*sqrt(3)
矢 = r - r/2
A0 = 2a + (s1 - s2)
B0 = 矢 + pi*2r
println(A0, ", ", B0)
term1 = sqrt(-12*sqrt(3)*A0 + 16*pi*A0 - 16*pi*B0 + 12*sqrt(3)*B0 - 16*sqrt(3)*pi - 4*sqrt(3) + 13 + 8*pi + 16*pi^2)
term2 = - 4*pi - 1 + 2*sqrt(3)
ans1 = (-term1 + term2)/(-4*pi + 3*sqrt(3))
ans2 = ( term1 + term2)/(-4*pi + 3*sqrt(3))
return (r > (-4*pi - 1 + 2*sqrt(3))/(-4*pi + 3*sqrt(3)) ? ans1 : ans2)
end;
半径が 1.2 の場合
rev(1.2)
4.731739518077568, 8.139822368615503
1.1999999999999988
半径が 3.45 の場合
rev(3.45)
27.906580792648718, 23.401989309769576
3.4499999999999993
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r = 1.2
a = r*√3/2
b = -r/2
c = 3r*(2 - √3)/2
d = r*(11 - 6*√3)/2
plot([a, 0, -a, a], [-r/2, r, -r/2, -r/2], color=:blue, lw=0.5)
circle(0, 0, r)
segment(0, -r/2, 0, -r, :green, lw=2)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(r, 0, "r ", :red, :right, :bottom, delta=delta/2)
point(a, -r/2, "(a,-r/2) ", :blue, :right, :bottom, delta=delta/2)
point(0, -r/2, " -r/2", :blue, :left, :bottom, delta=delta/2)
point(0, -3r/4, " 矢", :green, mark=false)
hline!([0], color=:gray, lw=0.5)
vline!([0], color=:gray, lw=0.5)
else
plot!(showaxis=false)
end
end;