算額（その1139）

算額（その1139）

二五 武州妻沼 聖天宮 文政11年(1828)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円10個，正三角形

正三角形の中に大円 3 個，小円 7 個を容れる。正三角形の一辺の長さが 20 寸のとき，小円の直径はいかほどか。

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

include("julia-source.txt")

using SymPy

@syms a::positive, r1::positive, x1::positive, r2::positive, x2::positive
x2 = 2x1
eq1 = x2 - x1 - 2sqrt(r1*r2)
eq2 = r1/(a - x1) - √Sym(3)/3
eq3 = r2/(a - x2) - √Sym(3)/3
res = solve([eq1, eq2, eq3], (r1, x1, r2))[1]

   (sqrt(3)*a/5, 2*a/5, sqrt(3)*a/15)

小円の直径 2r2 は，正三角形の一辺の長さ 2a の √3/15 倍である
正三角形の一辺の長さが 20 寸のとき，小円の直径は 20*√3/15 = 2.309401076758503 寸である。
ちなみに，大円の直径は小円の直径の 3 倍で，6.928203230275509 寸である。

function draw(a, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, x1, r2) = (sqrt(3)*a/5, 2*a/5, sqrt(3)*a/15)
   @printf("正三角形の一辺の長さが %g のとき，小円の直径は %g である（大円の直径は %g）。\n", 2a, 2r2, 2r1)
   @printf("a = %g;  r1 = %g;  x1 = %g;  r2 = %g\n", a, r1, x1, r2)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:blue, lw=0.5)
   circle2(x1, r1, r1)
   circle(0, √3a/3 + r2 + r1, r1)
   circle2(2x1, r2, r2, :green)
   circle(0, r2, r2, :green)
   circle(0, √3a/3, r2, :green)
   circle2(x1, 2r1 + r2, r2, :green)
   circle(0, √3a/3 + 2r2 + 2r1, r2, :green)
   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(x1, r1, "大円:r1,(x1,r1)", :red, :center, delta=-delta/2)
       point(0, √3a/3 + r2 + r1, "大円:r1\n(0,√3a/3+r2+r1)", :red, :center, delta=-delta/2)
       point(0, √3a/3, "小円:r2,(0,√3a/3)", :green, :left, :vcenter)
       point(0, r2, "小円:r2,(0,r2) ", :green, :right, :vcenter)
       point(2x1, r2, "小円:r2,(2x1,r2) ", :green, :right, :vcenter)
       point(x1, 2r1+r2, "小円:r2,(x1,2r1+r2) ", :green, :right, :vcenter)
       point(0, √3a/3 + 2r2 + 2r1, " 小円:r2,(0,√3a/3+2r2+2r1)", :green, :left, :vcenter)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(0, √3a, " √3a", :blue, :left, :vcenter)
   end
end;

算額（その1138）

算額（その1138）

二〇 武州不動岡村 不動堂 文政元年(1818)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円6個，外円，正五角形

外円の中に等円 5 個と正五角形を容れる。外円の直径が 15 寸，等円の直径が 3 寸のとき，正五角形の一辺の長さはいかほどか。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x, y)
正五角形の外接円の半径を a
とおき，以下の連立方程式を解く。

正五角形の一辺の長さは 2a*sin(π/5) である。

include("julia-source.txt")

using SymPy

@syms R::positive, a::positive, r::positive, x::positive, y::positive,
     x1::positive, y1::positive,
     x2::positive, y2::positive,
     x0::positive, y0::positive
(x1, y1) = (a*cosd(Sym(18)), a*sind(Sym(18)))
(x2, y2) = (a*cosd(Sym(54)), -a*sind(Sym(54)))
eq1 = x^2 + y^2 - (R - r)^2
eq2 = dist(x1, y1, 0, a, x, y) - r^2
eq3 = dist(x1, y1, x2, y2, x, y) - r^2;

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (a, x, y) = u
   return [
x^2 + y^2 - (R - r)^2,  # eq1
sqrt(5)*a^2/8 + 3*a^2/8 - a*x*sqrt(2*sqrt(5) + 10)/4 - 3*a*y/4 - sqrt(5)*a*y/4 - r^2 - sqrt(5)*x^2/8 + 5*x^2/8 + x*y*sqrt(2*sqrt(5) + 10)/4 + sqrt(5)*y^2/8 + 3*y^2/8,  # eq2
sqrt(5)*a^2/8 + 3*a^2/8 - 9*a*x*sqrt(10*sqrt(5) + 50)/16 - 19*a*x*sqrt(2*sqrt(5) + 10)/16 + 13*a*x*sqrt(10 - 2*sqrt(5))/8 + 3*a*x*sqrt(50 - 10*sqrt(5))/4 + a*y/2 - r^2 + sqrt(5)*x^2/8 + 5*x^2/8 - 11*x*y*sqrt(10*sqrt(5) + 50)/8 - 23*x*y*sqrt(2*sqrt(5) + 10)/8 + 37*x*y*sqrt(10 - 2*sqrt(5))/8 + 17*x*y*sqrt(50 - 10*sqrt(5))/8 - sqrt(5)*y^2/8 + 3*y^2/8,  # eq3
   ]
end;

R = 15/2
r = 3/2
iniv = BigFloat[5, 4, 4]
res = nls(H, ini=iniv)

   ([5.430242979853682, 4.375871302877152, 4.105088347484888], true)

正五角形の外接円の半径は 5.430242979853682 なので，正五角形の一辺の長さは 5.430242979853682*sind(36)*2 = 6.3836334798454555 寸である。

res[1][1]*sind(36)*2

   6.3836334798454555

図を描くために必要な座標値を計算する。

@syms R
eq11 = x0^2 + y0^2 - R^2
eq12 = y0 - y2 - (y1 - y2)/(x1 - x2)*(x0 - x2)
res2 = solve([eq11, eq12], (x0, y0))[1]

   (-a*sqrt(50 - 10*sqrt(5))/40 + a*sqrt(10*sqrt(5) + 50)/40 + a*sqrt(10 - 2*sqrt(5))/8 + a*sqrt(2*sqrt(5) + 10)/8 + (a/4 - sqrt(10)*sqrt(-4*sqrt(5)*R^2 + 20*R^2 - 5*a^2 - sqrt(5)*a^2)/(2*(-sqrt(5 - sqrt(5))*sqrt(sqrt(5) + 5) + 10)))*(-sqrt(10*sqrt(5) + 50)/10 + sqrt(50 - 10*sqrt(5))/10), -a/4 + sqrt(-40*sqrt(5)*R^2 + 200*R^2 - 50*a^2 - 10*sqrt(5)*a^2)/16 + sqrt(-8*sqrt(5)*R^2 + 40*R^2 - 10*a^2 - 2*sqrt(5)*a^2)/16)

function draw(R, r, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, x, y) = [5.430242979853682, 4.375871302877152, 4.105088347484888]
   (x1, y1) = (a*cosd(18), a*sind(18))
   (x2, y2) = (a*cosd(54), -a*sind(54))
   (x0, y0) = (-a*sqrt(50 - 10*sqrt(5))/40 + a*sqrt(10*sqrt(5) + 50)/40 + a*sqrt(10 - 2*sqrt(5))/8 + a*sqrt(2*sqrt(5) + 10)/8 + (a/4 - sqrt(10)*sqrt(-4*sqrt(5)*R^2 + 20*R^2 - 5*a^2 - sqrt(5)*a^2)/(2*(-sqrt(5 - sqrt(5))*sqrt(sqrt(5) + 5) + 10)))*(-sqrt(10*sqrt(5) + 50)/10 + sqrt(50 - 10*sqrt(5))/10), -a/4 + sqrt(-40*sqrt(5)*R^2 + 200*R^2 - 50*a^2 - 10*sqrt(5)*a^2)/16 + sqrt(-8*sqrt(5)*R^2 + 40*R^2 - 10*a^2 - 2*sqrt(5)*a^2)/16)
   θ1 = -54:72:306 
   θ2 = atand(y0, x0):72:432
   plot()
   circle(0, 0, R)
   rotate(x, y, r, angle=72, :blue)
   for i = 1:5
       segment(a*cosd(θ1[i]), a*sind(θ1[i]), R*cosd(θ2[i]), R*sind(θ2[i]), :green)
   end
   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(x1, y1, " (x1,y1)", :green, :left, :vcenter)
       point(x2, y2, "(x2,y2) ", :green, :right, :bottom, delta=delta/2)
       point(0, a, " a", :green, :left, :bottom, delta=delta/2)
       point(0, R, " R", :red, :left, :bottom, delta=delta/2)
       point(x, y, "等円:r,(x,y)", :blue, :center, delta=-delta/2)
   end
end;

算額（その1137）

算額（その1137）

二〇 武州不動岡村 不動堂 文政元年(1818)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円2個，半円1個，直角三角形

直角三角形の中に大円（半円），中円，小円を容れる。大円と中円の直径がそれぞれ 2601 寸，1156 寸のとき，小円の直径を求めよ。

直角三角形の鈎，股を a, b
大円の半径と中心座標を r1, (r1, 0)
中円の半径と中心座標を r2, (r2, y2)
小円の半径と中心座標を r3, (x3, r3)
とおき，以下の連立方程式を解く。

SymPy の性能上，r1, r2 を未知数のままとして解を求めることができない。

include("julia-source.txt")

using SymPy

@syms a::positive, b::positive, r1::positive,
     r2::positive, y2::positive, r3::positive, x3::positive
(r1, r2) = (2601, 1156).//2
y2 = 2*sqrt(r1)*sqrt(r2)
eq1 = dist2(a, 0, 0, b, r1, 0, r1)/a
eq2 = b*r2 - r1*(b - y2)
eq3 = b - y2 + 2sqrt(r1*r2) + 2sqrt(r1*r3) + a - x3 - sqrt(a^2 + b^2)
eq4 = (x3 - r1)^2 + r3^2 - (r1 + r3)^2 |> expand;
res = solve([eq1, eq2, eq3, eq4], (a, b, r3, x3))[1]

   (3147.42857142857, 3121.20000000000, 162.000000000000, 2754.00000000000)

大円と中円の直径がそれぞれ 2601 寸，1156 寸のとき，小円の直径は 324 寸である。

その他のパラメータは以下のとおりである。

   r1 = 1300.5;  r2 = 578;  a = 3147.43;  b = 3121.2;  r3 = 162;  x3 = 2754

function draw(r1, r2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, r3, x3) = (3147.42857142857, 3121.20000000000, 162.000000000000, 2754.00000000000)
   y2 = 2*sqrt(r1)*sqrt(r2)
   @printf("大円と中円の直径が %g，%g のとき，小円の直径は %g である。\n", 2r1, 2r2, 2r3)
   @printf("r1 = %g;  r2 = %g;  a = %g;  b = %g;  r3 = %.15g;  x3 = %g\n", r1, r2, a, b, r3, x3)
   plot([0, a, 0, 0], [0, 0, b, 0], color=:green, lw=0.5)
   circle(r1, 0, r1, beginangle=0, endangle=180)
   circle(r2, y2, r2, :blue)
   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", :green, :left, :bottom, delta=delta/2)
       point(0, b, " b", :green, :left, :bottom, delta=delta/2)
       point(r1, 0, "大円:r1,(r1,0)", :red, :center, :bottom, delta=delta/2)
       point(r2, y2, "中円:r2,(r2,y2)", :blue, :center, delta=-delta/2)
       point(x3, r3, "小円:r3,(x3,r3)", :magenta, :right, :vcenter, deltax=-4delta)
   end
end;