算額（その2047）

算額（その2047）

四十九 群馬県安中市板鼻 鷹巣神社 文政11年(1828)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

キーワード：円5個，長方形

長方形の中に，甲，乙，丙，丁，戊の 5 円を容れる。長方形の長辺が 38 寸のとき， 短辺はいかほどか。

長方形の長編，短辺を a, b
甲円の半径と中心座標を r1, (a - r1, r1)
乙円の半径と中心座標を r2, (r2, b - r2)
丙円の半径と中心座標を r3, (r3, r3)
丁円の半径と中心座標を r4, (x4, b - r4)
戊円の半径と中心座標を r5, (a - r5, b - r5)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r1::positive, r2::positive,
     r3::positive, r4::positive, x4::positive, r5::positive
@syms a, b, r1, r2, r3, r4, x4, r5
a = Sym(38)
eq1 = (a - r1 - r2)^2 + (b - r2 -r1)^2 - (r1 + r2)^2
eq2 = (a - r1 - r3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (a - r1 - x4)^2 + (r1 - b + r4)^2 - (r1 + r4)^2
eq4 = (r1 - r5)^2 + (b - r5 - r1)^2 - (r1 + r5)^2
eq5 = (r2 - r3)^2 + (b - r2 - r3)^2 - (r2 + r3)^2
eq6 = (x4 - r2)^2 + (r2 - r4)^2 - (r2 + r4)^2
eq7 = (a - r5 - x4)^2 + (r5 - r4)^2 - (r4 + r5)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (b, r1, r2, r3, r4, x4, r5))

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)
   (b, r1, r2, r3, r4, x4, r5) = u
   return [
       -(r1 + r2)^2 + (a - r1 - r2)^2 + (b - r1 - r2)^2,  # eq1
       (r1 - r3)^2 - (r1 + r3)^2 + (a - r1 - r3)^2,  # eq2
       -(r1 + r4)^2 + (a - r1 - x4)^2 + (-b + r1 + r4)^2,  # eq3
       (r1 - r5)^2 - (r1 + r5)^2 + (b - r1 - r5)^2,  # eq4
       (r2 - r3)^2 - (r2 + r3)^2 + (b - r2 - r3)^2,  # eq5
       (-r2 + x4)^2 + (r2 - r4)^2 - (r2 + r4)^2,  # eq6
       (-r4 + r5)^2 - (r4 + r5)^2 + (a - r5 - x4)^2,  # eq7
   ]
end;

a = 38
iniv = BigFloat[33, 12, 9, 7.5, 4.5, 23, 5.5]
res = nls(H, ini=iniv)

   ([33.00311463837703, 11.814579503859783, 9.106235741252021, 7.437509234825292, 4.897027068931824, 22.461906135655383, 5.325015237401601], true)

長方形の長辺が 38 寸のとき， 短辺は 33.00311463837703 寸である。

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

  a = 38;  b = 33.0031;  r1 = 11.8146;  r2 = 9.10624;  r3 = 7.43751;  r4 = 4.89703;  x4 = 22.4619;  r5 = 5.32502

function draw(a, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (b, r1, r2, r3, r4, x4, r5) = res[1]
   @printf("a = %g;  b = %g;  r1 = %g;  r2 = %g;  r3 = %g;  r4 = %g;  x4 = %g;  r5 = %g\n",
       a, b, r1, r2, r3, r4, x4, r5)
   plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:black, lw=0.5)
   circle(a - r1, r1, r1)
   circle(r2, b - r2, r2, :green)
   circle(r3, r3, r3, :blue)
   circle(x4, b - r4, r4, :magenta)
   circle(a - r5, b - r5, r5, :orange)
   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 - r1, r1, "甲円:r1,(a-r1,r1)", :red, :center, delta=-delta)
       point(r2, b - r2, "乙円:r2,(r2,b-r2)", :green, :center, delta=-delta)
       point(r3, r3, "丙円:r3,(r3,r3)", :blue, :center, delta=-delta)
       point(x4, b - r4, "丁円:r4\n(x4,b-r4)", :magenta, :center, delta=-delta)
       point(a - r5, b - r5, "戊円:r5\n(a-r5,b-r5)", :orange, :center, delta=-delta)
       point(a, b, "(a,b)", :black, :right, :bottom, delta=delta)
   end
end;

draw(38, true)

 

算額（その2046）

算額（その2046）

三十八 群馬県前橋市下大屋町 産泰神社 文政5年(1822)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

半円の円周内接する n 個の等円を描く。図の灰色部分の面積はいかほどか。

注：n は奇数とする。
半円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x, y); x = (R - r)*cos(pi/2n), y = (R - r)*sin(pi/2n)
とおき，以下の方程式を解いて等円の半径を求める。

include("julia-source.txt");

using SymPy

@syms n::positive, θ::positive, R::positive, r::positive, S::positive
θ = PI/2n
eq1 = sin(θ)*(R - r) - r

res = solve(eq1, r)[1]
res |> println

   R*sin(pi/(2*n))/(sin(pi/(2*n)) + 1)

求める面積は，図の右下の直角三角形の面積から白の扇型の面積を引き，2n 倍したものである。

S = (R - r)*cos(θ)*r/2 - PI*r^2*(PI/2 - θ)/2PI |> simplify
S |> println

   r*(2*n*(R - r)*cos(pi/(2*n)) + pi*r*(1 - n))/(4*n)

外円の半径が 1/2，等円の個数が 7 のとき，求める面積は 0.175958422104866 である。

14*S(R => 1/2, n=> 7, r => res(R => 1/2, n=> 7).evalf()).evalf() |> println

   0.175958422104866

function rotate2(ox, oy, r, color=:red; angle=120, beginangle=0, endangle=360, by=0.5, n=0)
   for deg in 0:angle:180-1
       (ox2, oy2) = [cosd(deg) -sind(deg); sind(deg) cosd(deg)] * [ox; oy]
       circlef(ox2, oy2, r, color; beginangle, endangle, by, n)
       circle(ox2, oy2, r, "red"; beginangle, endangle, by, n)
   end
end;

function draw(R, n, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   θ = 90/n
   r = R*sind(θ)/(sind(θ) + 1)
   S = r*(2*n*(R - r)*cosd(θ) + pi*r*(1 - n))/(4*n)
   @printf("n = %d;  R = %g;  r = %g;  S = %g\n", n, R, r, 2n*S)
   plot()
   circle(0, 0, R, :blue, beginangle=0, endangle=180)
   circlef(0, 0, (R - r)*cosd(θ), :gray70, beginangle=0, endangle=180)
   x = (R - r)*cosd(θ)
   y = (R - r)*sind(θ)
   rotate2(x, y, r, :white, angle=180/n)
   plot!([0, x, x], [0, y, 0], color=:black, lw=0.5)
   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, :bottom, delta=delta/2)
       point(x, y, "(x,y)", :red, :center, :bottom, delta=delta)
   end
   segment(-R, 0, R, 0, :blue)
end;

draw(1/2, 7, true)