算額（その290）

算額（その290）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(249)
長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

長方形内に右下の頂点で長辺に接する円弧の一部があり，他に大円，中円，小円がある。それぞれは図のように互いに接している。
中円の直径が 500 寸のとき，大円の直径を求めよ。

長方形の長辺の長さを a とする。
その一部が円弧である円の半径，中心座標を r0, (0, r0) とする。
長方形の短辺の長さは sqrt(r0^2 - a^2) である。
大円の半径，中心座標を r1, (a - r1, r1)
中円の半径，中心座標を r2, (x2, x2)
小円の半径，中心座標を r3, (x3, y3) および (a - r3, r3) とする。
以下の連立方程式を nlsolve() で解く。

include("julia-source.txt");

using SymPy

@syms a::positive, r0::positive, r1::positive, r2::positive,
     x2::positive, r3::positive, x3::positive, y3::positive;

r2 = 500 // 2
eq1 = (a - r1)^2 + (r0 - r1)^2 - (r0 + r1)^2
eq2 = x3^2 + (r0 - y3)^2 - (r0 + r3)^2
eq3 = x2^2 + (r0 - r2)^2 - (r0 + r2)^2
eq4 = 2(r1 - r3)^2 - (r1 + r3)^2
eq5 = (a - r1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq6 = (a - r1 - x3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq7 = (x3 - x2)^2 + (y3 - r2)^2 - (r2 + r3)^2;

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7])

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=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (a, r0, r1, x2, r3, x3, y3) = u
   return [
       (a - r1)^2 + (r0 - r1)^2 - (r0 + r1)^2,  # eq1
       x3^2 - (r0 + r3)^2 + (r0 - y3)^2,  # eq2
       x2^2 + (r0 - 250)^2 - (r0 + 250)^2,  # eq3
       2*(r1 - r3)^2 - (r1 + r3)^2,  # eq4
       (r1 - 250)^2 - (r1 + 250)^2 + (a - r1 - x2)^2,  # eq5
       (-r1 + y3)^2 - (r1 + r3)^2 + (a - r1 - x3)^2,  # eq6
       -(r3 + 250)^2 + (-x2 + x3)^2 + (y3 - 250)^2,  # eq7
      ]
end;

iniv = [big"3100.0", 3900, 450, 2000, 75, 2100, 550]
res = nls(H, ini=iniv);

names = [" a", "r0", "r1", "x2", "r3", "x3", "y3"]
for (i, name) in enumerate(names)
   @printf("%s = %.6f\n", name, res[1][i])
end

    a = 3086.670209
   r0 = 3867.992287
   r1 = 449.500799
   x2 = 1966.721202
   r3 = 77.122145
   x3 = 2118.355531
   y3 = 539.855012

中円の半径は 449.500799 で，直径はほぼ 899 寸である。

using Plots

function diamond(x0, y0, a)
   plot!([1, 0, -1, 0, 1] ./ √2 .* a .+ x0, [0, 1, 0, -1, 0] ./ √2 .* a .+ y0)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r0, r1, x2, r3, x3, y3) = res[1]
   r2 = 500 / 2
   b = r0 - sqrt(r0^2 - a^2)
   plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:black, lw=0.5)
   circle(0, r0, r0, beginangle=270, endangle=360)
   circle(a - r1, r1, r1, :blue)
   circle(x2, r2, r2, :green)
   circle(x3, y3, r3, :magenta)
   circle(a - r3, r3, r3, :magenta)
   if more
       point(a, b, "(a,b)")
       point(0, r0, " r0", :red)
       point(a - r1, r1, "大円:r1\n(a-r1,r1)", :blue, :center, :bottom)
       point(x2, r2, "中円:r2 \n(x2,r2) ", :green, :right, :top)
       point(x3, y3, "小円:r3  \n(x3,y3)  ", :magenta, :right, :bottom)
       point(a - r3, r3, "小円:r3  \n(a-r3,r3)  ", :magenta, :right, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false, xlims=(0, 1.1Float64(a)), ylims=(-20, 1.1Float64(b)))
   end
end;

算額（その289）

算額（その289）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(250)
長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

長方形内に斜線を描き，大円，小円が2個ずつ入っている。大円は長方形と斜線に接している。小円は長方形，大円，斜線に接している。
長方形の短辺の長さが 360 寸，小円の直径が 40 寸のとき，長方形の長辺の長さを求めよ。

長方形の長辺の長さを a とする。短辺の長さは大円の直径に等しい。

大円の半径と中心座標を r1, (r1, r1)
小円の半径と中心座標を r2, (x, 2r1 - r2) とする。
斜線と長方形の長辺の交点の座標を (b, 2r1), (a - b, 0) とする。
以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r1::positive, r2::positive, x::positive;

r1 = 360 // 2
r2 = 40 // 2
eq1 = (x - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = distance(b, 2r1, a - b, 0, r1, r1) - r1^2
eq3 = distance(b, 2r1, a - b, 0, x, 2r1 - r2) - r2^2

res = solve([eq1, eq2, eq3], (a, b, x))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (735, 315, 300)
    (960, 240, 300)

最初の組が適解である。すなわち，長方形の長辺は 735 寸である。

using Plots

function diamond(x0, y0, a)
   plot!([1, 0, -1, 0, 1] ./ √2 .* a .+ x0, [0, 1, 0, -1, 0] ./ √2 .* a .+ y0)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, x) = res[1]
   println("a = $a")
   plot([0, a, a, 0, 0], [0, 0, 2r1, 2r1, 0], color=:black, lw=0.5)
   circle(r1, r1, r1)
   circle(a - r1, r1, r1)
   circle(x, 2r1 - r2, r2, :blue)
   circle(a - x, r2, r2, :blue)
   segment(a - b, 0, b, 2r1, :green)
   if more
       point(r1, r1, "(r1,r1)", :red)
       point(x, 2r1 - r2, "(x,2r1-r2)  ", :blue, :right)
       point(b, 2r1, "(b,2r1)", :black, :left, :bottom)
       point(a, 0, " a", :black, :right, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その288）

算額（その288）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(250)
長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

円内に 2 本の平行な弦と 3 個の正方形が入っている。正方形の一辺の長さが 1 寸のとき，円の径を求めよ。

円の半径を R とする。正方形の一辺の長さを a とすれば，対角線の長さの半分は a/√2 = √2a/2 で，これを b とする。
⊿OAB において，OB = R, OA = 3b - R, AB = 2b

以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, a::positive, b::positive

b = a/√Sym(2)
eq = (3b-R)^2 + (2b)^2 - R^2
solve(eq, R)[1] |> println  # 半径 R

   13*sqrt(2)*a/12

直径は a*13√2/6 で，a = 1 寸のとき，円の直径は 3.0641293851417064 寸。

13√2/6

   3.0641293851417064

using Plots

function diamond(x0, y0, a)
   plot!([1, 0, -1, 0, 1] ./ √2 .* a .+ x0, [0, 1, 0, -1, 0] ./ √2 .* a .+ y0)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1
   b = a/√2
   R = a*13*sqrt(2)/12
   println("R = $R  2R = $(2R)")
   plot()
   circle(0, 0, R)
   diamond(0, R - b, a)
   diamond(b, R - 3b, a)
   diamond(-b, R - 3b, a)
   y = R - 2b
   x = sqrt(R^2 - y^2)
   segment(-x, y, x, y, :blue)
   y = R - 4b
   x = sqrt(R^2 - y^2)
   segment(-x, y, x, y, :blue)
   if more
       point(0, R - b, "  R-b")
       point(0, R - 3b, "  A")
       point(0, R - 3b, "R-3b  ", :green, :right)
       point(2b, 0, "2b", :green, :center)
       point(R, 0, " R", :green)
       point(0, 0, " O")
       point(2b, R - 3b, " B")
       plot!([0, 2b, 0, 0], [0, R - 3b, R - 3b, 0], color=:gray, linestyle=:dot)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;