算額(その1293)
十七 群馬県高崎市八幡町 八幡宮 文化7年(1810)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:円5個,菱形
菱形の中に,大円 2 個,中円 2 個,小円 1 個を容れる。菱形の対角線の長い方が 12 寸,短い方が 5 寸のとき,小円の直径はいかほどか。
菱形の対角線を 2a, 2b; a > b
大円の半径と中心座標を r1, (r3 + r1, 0)
中円の半径と中心座標を r2, (0, r3 + r2)
小円の半径と中心座標を r3, (0, 0)
とおき,以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms a::positive, b::positive,
r1::positive, r2::positive, r3::positive;
sinθ = b/sqrt(a^2 + b^2)
cosθ = a/sqrt(a^2 + b^2)
eq1 = (r3 + r1)^2 + (r3 + r2)^2 - (r1 + r2)^2 |> expand
eq2 = (a - r3 - r1)*sinθ - r1
eq3 = (b - r3 - r2)*cosθ - r2;
SymPy では能力的に一度に解けないので,逐次解いていく。
まず,eq1 を解いて r1 を求める。
ans_r1 = solve(eq1, r1)[1]
ans_r1 |> println
r3*(r2 + r3)/(r2 - r3)
eq2 の r1 に,上で求めた r1 を代入し,r2 を求める。
eq12 = eq2(r1 => ans_r1)
ans_r2 = solve(eq12, r2)[1]
ans_r2 |> println
r3*(-a*b - r3*sqrt(a^2 + b^2))/(-a*b + 2*b*r3 + r3*sqrt(a^2 + b^2))
eq3 に,上で求めた r1, r2 を代入し,r3 を求める。
eq13 = eq3(r1 => ans_r1, r2 => ans_r2) |> simplify |> numerator
ans_r3 = solve(eq13, r3)[2] |> factor
ans_r3 |> println
-a*b*(a + b + sqrt(a^2 + b^2) - sqrt(3*a^2 + 2*a*sqrt(a^2 + b^2) + 3*b^2 + 2*b*sqrt(a^2 + b^2)))/(a - b)^2
解として得られた r3 を 代入して r2,r1 を確定する。
a = 12/2
b = 5/2
r3 = -a*b*(a + b + sqrt(a^2 + b^2) - sqrt(3*a^2 + 2*a*sqrt(a^2 + b^2) + 3*b^2 + 2*b*sqrt(a^2 + b^2)))/(a - b)^2
r2 = r3*(-a*b - r3*sqrt(a^2 + b^2))/(-a*b + 2*b*r3 + r3*sqrt(a^2 + b^2))
r1 = r3*(r2 + r3)/(r2 - r3)
(r1, r2, r3) |> println
(1.5296184351192643, 0.9631806558860881, 0.49337363357064845)
菱形の対角線の長い方が 12 寸,短い方が 5 寸のとき,小円の直径は 2*0.49337363357064845 = 0.9867472671412969 寸である。「答」では「小円の径一寸三分二厘半微弱」としているが,正しいとは思えない。
すべてのパラメータは以下のとおりである。
a = 6; b = 2.5; r1 = 1.52962; r2 = 0.963181; r3 = 0.493374
function draw(a, b, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r3 = -a*b*(a + b + sqrt(a^2 + b^2) - sqrt(3*a^2 + 2*a*sqrt(a^2 + b^2) + 3*b^2 + 2*b*sqrt(a^2 + b^2)))/(a - b)^2
r2 = r3*(-a*b - r3*sqrt(a^2 + b^2))/(-a*b + 2*b*r3 + r3*sqrt(a^2 + b^2))
r1 = r3*(r2 + r3)/(r2 - r3)
@printf("菱形の対角線の長い方が %g,短い方が %g のとき,小円の直径は %g である。\n", 2a, 2b, 2r3)
@printf("a = %g; b = %g; r1 = %g; r2 = %g; r3 = %g\n", a, b, r1, r2, r3)
plot([a , 0, -a, 0, a], [0, b, 0, -b, 0], color=:green, lw=0.5)
circle2(r3 + r1, 0, r1)
circle22(0, r3 + r2, r2, :blue)
circle(0, 0, r3, :magenta)
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(a, 0, " a", :green, :left, :vcenter)
point(0, b, "b", :green, :center, :bottom, delta=delta)
point(r3 + r1, 0, "大円:r1\n(r3+r1,0)", :red, :center, delta=-delta)
point(0, r3 + r2, "中円:r2\n(0,r3+r2)", :blue, :center, :bottom, delta=delta)
point(0, 0, "小円:r3,(0,0)", :magenta, :right, :vcenter, deltax=-7delta)
end
end;
draw(12/2, 5/2, true)