算額（その1267）

算額（その1267）

百三十 現存するが奉納場所不明 明治14年(1881)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円3個，正三角形，正方形

正三角形と正方形が交わっており，その区分領域に大円，中円，小円を容れる。正方形の一辺の長さが与えられたとき小円の直径を求めよ。

正三角形の一辺の長さを a
正方形の一辺の長さを b
大円の半径と中心座標を r1, (r1, r1)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (x3, r3)
とおき，以下の連立方程式を解く。

1.　a, r2, x2 を求める

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     r1::positive, r2::positive, x2::positive;
# r1 = b/2
b = 2r1
eq1 = b/(a - b) - √Sym(3)
eq2 = (a - x2)/r2 - √Sym(3)
eq3 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2 |> expand;

ans_a = solve(eq1, a)[1]
ans_a |> println

   2*r1*(sqrt(3) + 3)/3

eq12 = eq2(a => ans_a) |> simplify;
eq13 = eq3(a => ans_a);
res1 = solve([eq12, eq13], (r2, x2))[1]

   (r1*(-2*sqrt(sqrt(3) + 3) + sqrt(3) + 4)/3, r1*(-2*sqrt(3) + 3 + 2*sqrt(3*sqrt(3) + 9))/3)

2.　r3, x3 を求める

前節で求めた a, r2, x2 を既知として r3, x3 を求める。

@syms a::positive, b::positive,
     r1::positive, r2::positive, x2::positive,
     r3::positive, x3::positive;
b = 2r1
# r1 = b/2
a = 2r1*(√Sym(3) + 3)/3
r2 = r1*(-2sqrt(√Sym(3) + 3) + √Sym(3) + 4)/3
x2 = r1*(-2√Sym(3) + 3 + 2sqrt(3√Sym(3) + 9))/3
eq4 = (x3 - r1)^2 + (r1 - r3)^2 - (r1 + r3)^2 |> expand
eq5 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2 |> expand;
res2 = solve([eq4, eq5], r3, x3)[2];  # 2 of 2

@syms d
ans_r3 = apart(res2[1]) |> simplify
ans_r3 |> println

   r1*(-204*sqrt(sqrt(3) + 3) - 92*sqrt(3*sqrt(3) + 9) + 229*sqrt(3) + 449)/338

r1 = b/2 として，r3 は r1 の (-204*sqrt(sqrt(3) + 3) - 92*sqrt(3*sqrt(3) + 9) + 229*sqrt(3) + 449)/338 倍である。

たとえば，b = 10 のとき，r3 = 0.8171180745896389，直径は 1.6342361491792778 である。

ans_x3 = apart(res2[2]) |> simplify
ans_x3 |> println

   r1*(-9*sqrt(sqrt(3) + 3) - sqrt(3*sqrt(3) + 9) + 8*sqrt(3) + 33)/13

b = 10 のとき，各パラメータは以下のとおりである。

   b = 10;  a = 15.7735;  r1 = 5;  r2 = 2.30233;  x2 = 11.7858;  r3 = 0.817118;  x3 = 9.04257

function draw(b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = b/2
   a  = 2r1*(√3 + 3)/3
   r2 = r1*(√3 + 4 - 2sqrt(√3 + 3))/3
   x2 = r1*(2sqrt(3√3 + 9) -2√3 + 3)/3
   r3 = r1*(229√3 + 449 - 204sqrt(√3 + 3) - 92sqrt(3√3 + 9))/338
   x3 = r1*(8√3 + 33 - 9sqrt(√3 + 3) - sqrt(3√3 + 9))/13
   @printf("b = %g;  a = %g;  r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", b, a, r1, r2, x2, r3, x3)
   plot([0, a, a/2, 0], [0, 0, √3a/2, 0], color=:blue, lw=0.5)
   plot!([0, b, b, 0, 0], [0, 0, b, b, 0], color=:green, lw=0.5) 
   circle(r1, r1, r1)
   circle(x2, r2, r2, :orange)
   circle(x3, r3, 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", :blue, :center, delta=-delta)
       point(0, b, "b ", :green, :right, :vcenter)
       point(b, 0, "b", :green, :left, delta=-delta)
       point(a/2, √3a/2, "(a/2,√3a/2)", :blue, :left, :bottom, delta=delta/2)
       point(r1, r1, "大円:r1,(r1,r1)", :red, :center, delta=-delta)
       point(x2, r2, "中円:r2\n(x2,r2)", :orange, :center, delta=-delta)
       point(x3, r3, "小円:r3,(x3,r3)", :magenta, :right, :bottom, delta=2delta, deltax=-4delta)
       plot!(xlims=(-5delta, a + 4delta), ylims=(-5delta, √3a/2 + 5delta))
   end
end;

draw(10, true)