算額（その735）

算額（その735）

福島県田村郡三春町 龍穏院 明治26年(1893)
～落書き帳「○△□」～ 64．手作りの封筒
http://streetwasan.web.fc2.com/math15.5.22.html

一辺の長さが 3 寸の正方形の中に斜線 4 本を描き，区画された領域に宇円4個，宙円 1 個を入れる。それぞれの円の直径を求めよ。

正方形の一辺の長さを 2a
宇円の半径と中心座標を r1, (a - r1, 0)
宙円の半径と中心座標を r2, (0, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive, d
eq1 = dist(0, 2r1 - a, a, a, a - r1, 0) - r1^2
eq2 = dist(0, 2r1 - a, a, a, 0, 0) - r2^2
eq1 = numerator(apart(eq1, d))  # a^2*(-a + r1)*(-a + 3*r1)
eq2 = numerator(apart(eq2, d))
res = solve([eq1, eq2], (r1, r2))

   2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (a/3, a/5)
    (a, a)

2 組の解が求まるが，最初のものが適解である。

正方形の一辺の長さを 2a としたとき，宇円，宙円の半径は a/3, a/5 である。
正方形の一辺の長さが 3 寸のとき，宇円，宙円の直径は 1 寸，6 分である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 3//2
   (r1, r2) = (1/2, 3/10)
   @printf("正方形の一辺の長さ = %g;  宇円の直径 = %g;  宙円の直径 = %g\n", 2a, 2r1, 2r2)
   plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:black, lw=0.5)
   circle(0, 0, r2, :blue)
   circle2(a - r1, 0, r1)
   circle22(0, a - r1, r1)
   plot!([-a, 0, a], [a, 2r1 - a, a], color=:green, lw=0.5)
   plot!([-a, 0, a], [-a, a - 2r1, -a], color=:green, 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(0, a, "a", :black, :center, :bottom, delta=delta)
       point(0, r2, "r2", :black, :center, :bottom, delta=delta)
       point(0, a - r1, " a-r1", :red, :center, :bottom, delta=delta)
       point(0, 0, "宙円:r2\n(0,0)", :blue, :center, :bottom, delta=delta)
       point(0, a - 2r1, "a-2r1", :red, :center, :bottom, delta=delta)
       point(a - r1, 0, "宇円:r1,(a-2r1,0)", :red, :center, :bottom, delta=delta)
       point(a, 0, " a", :black, :left, :vcenter)
   end
end;

算額（その734）

算額（その734）

福島県安達郡東和町木幡山（現二本松市） 隠津島神社 明治17年(1884)
～落書き帳「○△□」～ 16．また一寸（○△□のおにぎり弁当）
http://streetwasan.web.fc2.com/math15.5.23.html

直角三角形の中に，正方形，円，正三角形が互いに接して入っている。鈎の長さが 1 寸のとき，「子」の長さはいかほどか。

「鈎」，「股」，「弦」をそのまま変数名として使う。
正三角形の一辺の長さを a，正方形の一辺の長さを b，円の半径を r として以下の方程式を解く。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
     r::positive, a::positive, b::positive,
     xc::positive, x::positive, y::positive,
     x0::positive, y0::positive, d
弦 = sqrt(鈎^2 + 股^2)
eq1 = 股*y - xc*鈎
eq2 = sqrt(Sym(3))a/(股 - (x + a)) - 鈎/股
eq2 = sqrt(Sym(3))a*股 - 鈎*(股 - (x + a))
eq3 = b*弦 - 股*(鈎 - y)
eq4 = dist(x, 0, x + a, sqrt(Sym(3))a, r, r) - r^2
eq4 = numerator(apart(eq4, d))
eq5 = sqrt((鈎 - y)^2 - b^2) + b + sqrt(3a^2 + (股 - x - a)^2) - 弦
eq6 = y + xc - sqrt(y^2 + xc^2) - 2r
eq7 = dist(x, 0, x + a, sqrt(Sym(3))a, 0, y) - b^2
eq7 = numerator(apart(eq7, d))
eq8 = x0^2 + (y0 - y)^2 - b^2
eq9 = (鈎 - y0)/x0 - 鈎/股;

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)
   (股, r, a, b, x, y, xc, x0, y0) = u
   return [
       -xc + y*股,  # eq1
       sqrt(3)*a*股 + a + x - 股,  # eq2
       b*sqrt(股^2 + 1) - 股*(1 - y),  # eq3
       -sqrt(3)*r^2/2 - 3*r*x/2 + sqrt(3)*r*x/2 + 3*x^2/4,  # eq4
       b + sqrt(3*a^2 + (-a - x + 股)^2) + sqrt(-b^2 + (1 - y)^2) - sqrt(股^2 + 1),  # eq5
       -2*r + xc + y - sqrt(xc^2 + y^2),  # eq6
       -b^2 + 3*x^2/4 + sqrt(3)*x*y/2 + y^2/4,  # eq7
       -b^2 + x0^2 + (-y + y0)^2,  # eq8
       -1/股 + (1 - y0)/x0,  # eq9
   ]
end;

鈎 = 1
iniv = BigFloat[1.7, 0.17, 0.37, 0.5, 0.27, 0.46, 0.8, 0.23, 0.87]
res = nls(H, ini=iniv)

   ([1.7320508075688772, 0.16987298107780677, 0.36602540378443865, 0.4641016151377546, 0.2679491924311227, 0.4641016151377546, 0.8038475772933681, 0.2320508075688773, 0.8660254037844386], true)

直角三角形は辺の比が 2:1:√3 のよく見かけるものだ。
鈎が 1 寸のとき，「子」も 1 寸である。

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

   子 = 1.0, 弦 = 1.9999999999999998, b = 0.4641016151377546
   股 = 1.73205;  r = 0.169873;  a = 0.366025;  b = 0.464102;  x = 0.267949;  y = 0.464102;  xc = 0.803848;  x0 = 0.232051;  y0 = 0.866025

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   鈎 = 1
   (股, r, a, b, x, y, xc, x0, y0) = res[1]
   弦 = sqrt(鈎^2 + 股^2)
   println("子 = $(x + 2a), 弦 = $弦, b = $b")
   @printf("股 = %g;  r = %g;  a = %g;  b = %g;  x = %g;  y = %g;  xc = %g;  x0 = %g;  y0 = %g\n", 股, r, a, b, x, y, xc, x0, y0)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(r, r, r)
   segment(x, 0, x + a, √3a, :blue)
   segment(x + 2a, 0, x + a, √3a, :blue)
   segment(0, y, xc, 0)
   segment(0, y, x0, y0)
   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, r, "r,(r,r)", :red, :center, delta=-2delta)
       point(x0, y0, " (x0,y0)", :black, :left, :bottom)
       point(x + a, √3a, " (x+a,√3a)", :black, :left, :bottom)
       point(0, 鈎, "鈎 ", :black, :right, :vcenter)
       point(0, y, "y ", :black, :right, :vcenter)
       point(股, 0, "股", :black, :center, delta=-2delta)
       point(x, 0, "x", :black, :center, delta=-2delta)
       point(x + 2a, 0, "x+2a", :black, :center, delta=-2delta)
       point(xc, 0, "xc", :black, :center, delta=-2delta)
       plot!(xlims=(-0.15, 1.9), ylims=(-0.15, 1.1))
       arrow(0, 0, x + 2a, 0, :both)
       # plot!([0, x + 2a], [-0.1, -0.1], c=:red, arrow=:both)#legend=false, arrow=(:both, :open))
       dimension_line(0.0, -0.1, x + 2a, -0.1, "子", color=:red, horizontal=:center, vertical=:bottom, delta=delta/2)
   end
end;

算額（その733）

算額（その733）

福島県安達郡東和町木幡山（現二本松市） 隠津島神社 明治17年(1884)
～落書き帳「○△□」～ 15．一寸
http://streetwasan.web.fc2.com/math15.5.22.html

半径が 1 寸の扇の中に，甲半円と乙円が入っている。乙円の直径が最大になるとき，甲円の直径はいかほどか。

扇の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0. R - r1)
乙円の半径と中心座標を r2, (x2, R- r1 + r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive, x2::positive
eq1 = x2^2 + r2^2 - (r1 + r2)^2
eq2 = x2^2 + (R - r1 + r2)^2 - (R - r2)^2
res = solve([eq1, eq2], (r1, x2))

   2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    ((-R^(3/2)*sqrt(R - 8*r2)/2 - sqrt(R)*r2*sqrt(R - 8*r2) + R^2/2 + R*r2)/(R + 2*r2), sqrt(-R^(3/2)*sqrt(R - 8*r2)/2 - sqrt(R)*r2*sqrt(R - 8*r2) + R^2/2 - R*r2))
    ((R^(3/2)*sqrt(R - 8*r2)/2 + sqrt(R)*r2*sqrt(R - 8*r2) + R^2/2 + R*r2)/(R + 2*r2), sqrt(R^(3/2)*sqrt(R - 8*r2)/2 + sqrt(R)*r2*sqrt(R - 8*r2) + R^2/2 - R*r2))

2 組の解が求まるが，基本的にはどちらも適解である。
どちらの解も，R - 8r2 の平方根を求める操作が必要なので，乙円の半径 r2 は扇の半径 R の 1/8 よりも大きくはなれない。すなわち r2 の最大値は R/8 であり，このときの r1 が求める答えである。

扇の半径が 1 寸ならば，r2 の最大値は 1/8 = 0.125 であり，このとき r1 は 0.5 寸（直径は 1 寸）である。

res[1][1](R => 1, r2 => 0.125) |> println

   0.500000000000000

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 1
   r2 = R/8
   (r1, x2) = ((-R^(3/2)*sqrt(R - 8*r2)/2 - sqrt(R)*r2*sqrt(R - 8*r2) + R^2/2 + R*r2)/(R + 2*r2), sqrt(-R^(3/2)*sqrt(R - 8*r2)/2 - sqrt(R)*r2*sqrt(R - 8*r2) + R^2/2 - R*r2))
   @printf("扇の半径 = %g;  乙円の直径 = %g;  甲円の直径 = %g\n", R, 2r2, 2r1)
   r3 = R - r1
   y = R - r1
   x = sqrt(R^2 - y^2)
   θ = atand(y/ x)
   plot()
   circle(0, 0, R, beginangle=θ, endangle=180-θ, n=500)
   circle(0, 0, r3, :magenta, beginangle=θ, endangle=180-θ, n=500)
   circle(0, y, r1, :gray80)
   circle(0, y, r1, :blue, beginangle=0, endangle=180)
   circle2(x2, y + r2, r2, :green)
   segment(-x, y, x, y)
   segment(0, 0, x, y)
   segment(0, 0, -x, y)
   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(0, R, "R", :red, :center, :bottom, delta=delta)
       point(0, R - r1, "甲半円:r1,(0,R-r1)", :blue, :center, :bottom, delta=delta)
       point(x2, R - r1 + r2, "乙円:r2,(x2,R-r1+r2) ", :green, :right, :vcenter)
   end
end;