算額（その248）

算額（その248）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（229)
長野県諏訪市下諏訪 諏訪大社下社 慶応4年(1868)

長方形内に大円 1 個，中円，小円がそれぞれ 4 個入っている。大円の径が 1 寸のとき，長方形の長辺の長さはいかほどか。

長方形の長辺，短辺の長さをそれぞれ x，y とする。
第 1 象限の図形を対象とする。
大円の中心を原点とする。
大円の半径，中心座標を r1, (0, 0) とする。y = 2r1 である。
中円の半径，中心座標を r2, (x/2 - r2, y/2 - r2) とする。
小円の半径，中心座標を r3, (x/2 - r3, r3) とする。

問では「大円の径が 1 寸のとき」としているが，術では大円と長方形の長辺の関係を述べているので，以下では r1 も変数として連立方程式を解く。解は，r1 を含む式になる。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive,
   x::positive, y::positive;

r1 = 1//2
y = 2r1
eq1 = (x/2 - r2)^2 + (y/2 - r2)^2 - (r1 + r2)^2
eq2 = (x/2 - r3)^2 + r3^2 - (r1 + r3)^2
eq3 = (r2 - r3)^2 + (y/2 - r2 - r3)^2 - (r2 + r3)^2;

res = solve([eq1, eq2, eq3], [r2, r3, x])

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (4*sqrt(2)*r1/49 + 9*r1/49, 4*r1*(11 - 6*sqrt(2))/49, 2*r1*(23/49 + 32*sqrt(2)/49))

r1 = 0.5 とすると

   r1 = 0.50000;  r2 = 0.14956;  r3 = 0.10264;  x = 1.39296;  y = 1.00000

長辺の長さは 1.39296 である。

x は 2\*r1\*(23/49 + 32\*sqrt(2)/49) である。2r1 は大円の直径（径）であるので， x は 大円の径の (23//49 + 32\*sqrt(Sym(2))/49) 倍である。

2*r1*(23//49 + 32*sqrt(Sym(2))/49) |> simplify |> println

   2*r1*(23 + 32*sqrt(2))/49

術は「2048 の平方根に，23 を加え，大円の径を掛けて，49 で割る」

(sqrt(Sym(2048)) + 23) * 2r1 / 49 |> println

   2*r1*(23 + 32*sqrt(2))/49

同じである。

using Plots
using Printf

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1/2
   y = 2r1
   (r2, r3, x) = (4*sqrt(2)*r1/49 + 9*r1/49, 4*r1*(11 - 6*sqrt(2))/49, 2*r1*(23/49 + 32*sqrt(2)/49))
   @printf("r1 = %.5f;  r2 = %.5f;  r3 = %.5f;  x = %.5f;  y = %.5f\n", r1, r2, r3, x, y)
   plot([x/2, x/2, -x/2, -x/2, x/2], [-y/2, y/2, y/2, -y/2, -y/2], color=:black, lw=0.5)
   circle(0, 0, r1)
   circle4(x/2 - r2, y/2 - r2, r2, :blue)
   circle4(x/2 - r3, r3, r3, :green)
   if more == true
       point(0, 0, " 大円", :red, :left, :bottom)
       point(x/2 - r2, y/2 - r2, "中円", :blue, :left)
       point(x/2 - r3, r3, "小円", :green, :left)
       point(x/2, y/2, "(x/2,y/2)", :black, :right, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その247）

算額（その247）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（229)
長野県諏訪市下諏訪 諏訪大社下社 慶応4年(1868)

菱形の中に日円，月円，星円がそれぞれ 1 個，3 個，2 個入っている。
日円，月円の径を知って，星円の径を求めよ。

菱形の中心を原点とする。x 軸上の頂点の x 座標を x，y 軸上の頂点の y 座標を y とする。
日円の半径，中心座標を r1, (0, y1) とする。
真ん中の月円の半径，中心座標を r2, (0, y2) とする。
右の月円の半径，中心座標を r3, (x3, y3) とする。r3 = r2 である。
星円の半径，中心座標を r4, (x4, r4) とする。

以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, y1::positive, r2::positive, y2::negative,
   r3::positive, x3::positive, y3::positive,
   r4::positive, x4::positive, y4::positive,
   x::positive, y::positive;

r3 = r2
eq1 = x4^2 + (y1 - y4)^2 - (r1 + r4)^2
eq2 = (x3 - x4)^2 + (y3 - y4)^2 - (r3 + r4)^2
eq3 = x3^2 + (y1 - y3)^2 - (r1 + r3)^2
eq4 = x3^2 + (y2 - y3)^2 - (r2 + r3)^2
eq5 = distance(0, y, x, 0, 0, y1) - r1^2
eq6 = distance(0, y, x, 0, x4, y4) - r4^2
eq7 = distance(0, -y, x, 0, 0, y2) - r2^2
eq8 = distance(0, -y, x, 0, x3, y3) - r3^2
eq9 = (sqrt(x^2 + y^2)/x + 1) * (r1 + r2) - 2y;

# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9])

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
println(eq5, ",  # eq5")
println(eq6, ",  # eq6")
println(eq7, ",  # eq7")
println(eq8, ",  # eq8")
println(eq9, ",  # eq9")

   x4^2 - (r1 + r4)^2 + (y1 - y4)^2,  # eq1
   -(r2 + r4)^2 + (x3 - x4)^2 + (y3 - y4)^2,  # eq2
   x3^2 - (r1 + r2)^2 + (y1 - y3)^2,  # eq3
   -4*r2^2 + x3^2 + (y2 - y3)^2,  # eq4
   -r1^2 + x^2*y^2*(y - y1)^2/(x^2 + y^2)^2 + (-y*(x^2 + y*y1)/(x^2 + y^2) + y1)^2,  # eq5
   -r4^2 + (-x*(x*x4 + y^2 - y*y4)/(x^2 + y^2) + x4)^2 + (-y*(x^2 - x*x4 + y*y4)/(x^2 + y^2) + y4)^2,  # eq6
   -r2^2 + x^2*y^2*(y + y2)^2/(x^2 + y^2)^2 + (-y*(-x^2 + y*y2)/(x^2 + y^2) + y2)^2,  # eq7
   -r2^2 + (-x*(x*x3 + y^2 + y*y3)/(x^2 + y^2) + x3)^2 + (-y*(-x^2 + x*x3 + y*y3)/(x^2 + y^2) + y3)^2,  # eq8
   -2*y + (1 + sqrt(x^2 + y^2)/x)*(r1 + r2),  # eq9

日円，月円の半径を 5, 4 とする。

using Printf
using Plots

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)
   (r1, r2) = (5, 4)
   (r4, y1, y2, x3, y3, x4, y4, y, x) = u
   return [
       x4^2 - (r1 + r4)^2 + (y1 - y4)^2,  # eq1
       -(r2 + r4)^2 + (x3 - x4)^2 + (y3 - y4)^2,  # eq2
       x3^2 - (r1 + r2)^2 + (y1 - y3)^2,  # eq3
       -4*r2^2 + x3^2 + (y2 - y3)^2,  # eq4
       -r1^2 + x^2*y^2*(y - y1)^2/(x^2 + y^2)^2 + (-y*(x^2 + y*y1)/(x^2 + y^2) + y1)^2,  # eq5
       -r4^2 + (-x*(x*x4 + y^2 - y*y4)/(x^2 + y^2) + x4)^2 + (-y*(x^2 - x*x4 + y*y4)/(x^2 + y^2) + y4)^2,  # eq6
       -r2^2 + x^2*y^2*(y + y2)^2/(x^2 + y^2)^2 + (-y*(-x^2 + y*y2)/(x^2 + y^2) + y2)^2,  # eq7
       -r2^2 + (-x*(x*x3 + y^2 + y*y3)/(x^2 + y^2) + x3)^2 + (-y*(-x^2 + x*x3 + y*y3)/(x^2 + y^2) + y3)^2,  # eq8
       -2*y + (1 + sqrt(x^2 + y^2)/x)*(r1 + r2),  # eq9
   ]
end;

iniv = [big"1.0", 1.9, -3.1, 3.1, -1.3, 4, 1.9, 5.4, 9.3]
res = nls(H, ini=iniv);
println(res[2] ? "収束した" : "収束していない")

   収束した

names = ["r4", "y1", "y2", "x3", "y3", "x4", "y4", "y ", "x "]
for i in 1:length(names)
   println("$(names[i]) = $(Float64(res[1][i]))")
end

   r4 = 1.7470608602667996
   y1 = 3.9418436943485617
   y2 = -5.058156305651438
   x3 = 7.166451331820933
   y3 = -1.5026007500958825
   x4 = 6.740960675000237
   y4 = 4.228687607033678
   y  = 9.523406750862943
   x  = 19.195039966835868

星円の半径は 1.7470608602667996。

using Plots
using Printf

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r4, y1, y2, x3, y3, x4, y4, y, x) = res[1]
   (r1, r2) = (5, 4)
   r3 = r2
   # y = 5/2 + 5*sqrt(3)/3
   # x = sqrt(3) * y
   @printf("r1 = %.5f;  r2 = %.5f;  \nr4 = %.5f;  y1 = %.5f;  y2 = %.5f;  \nx3 = %.5f;  y3 = %.5f;  x4 = %.5f;  \ny4 = %.5f;  y = %.5f;  x = %.5f\n", r1, r2, r4, y1, y2, x3, y3, x4, y4, y, x)
   plot([x, 0, -x, 0, x], [0, y, 0, -y, 0], color=:black, lw=0.5)
   circle(0, y1, r1)
   circle(0, y2, r2, :blue)
   circle(x3, y3, r3, :blue)
   circle(-x3, y3, r3, :blue)
   circle(x4, y4, r4, :green)
   circle(-x4, y4, r4, :green)
   if more == true
       point(0, y1, " y1\n 日円", :red, :left, :bottom)
       point(0, y2, " y2\n 月円", :blue, :left, :bottom)
       point(x3, y3, "(x3,y3)\n月円", :blue, :left, :center)
       point(x4, y4, "(x4,y4)\n星円", :green, :left, :center)
       point(0, y, " y", :black, :left, :bottom)
       point(x, 0, "x", :black)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その246）

算額（その246）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（229)
長野県諏訪市下諏訪 諏訪大社下社 慶応4年(1868)

部分円の中に，大円，小円が 2 個ずつ入っている。大円の径が 147 寸，元の長さが 294 寸のとき，小円の径を求めよ。

弦を x 軸上に取り，その中点を原点とする。弦の右端の x 座標を (x, 0) とする。x = 294/2 である。
部分円の半径を R 中心座標を (0, y) とする。
右の大円の半径，中心座標を r2, (r2, r2)　とする。r2 = 147/2 である。
右の小円の半径，中心座標を r1, (x1, r1) とする。

以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms x::positive, y::positive, r1::positive, x1::positive, r2::positive, R::positive;

x = 294 // 2
r2 = 147 // 2
eq1 = x^2 + y^2 - R^2
eq2 = r2^2 + (r2 - y)^2 - (R - r2)^2
eq3 = x1^2 + (y - r1)^2 - (R - r1)^2
eq4 = (x1 - r2)^2 + (r2 - r1)^2 - (r1 + r2)^2;

res = solve([eq1, eq2, eq3, eq4], (r1, x1, R, y))

   1-element Vector{NTuple{4, Sym}}:
    (27/2, 273/2, 1225/8, 343/8)

r1 = 13.50000;  x1 = 136.50000;  R = 153.12500;  y = 42.87500
小円の半径は 13.5。元の単位で直径は 27 寸である。

eq1 ~ eq4 を x, r2 も未知数（変数）として解くと，r1 は以下の式で求めることができる。

x = 294 // 2
r2 = 147 // 2
num = r2 * (x - r2)^2 * (x + r2)^2
den = (3r2^2 + x^2)^2
r1 = num / den
println(r1)

   27//2

using Plots
using Printf

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   # (r1, x1, R, y) = res[1]
   (r1, x1, R, y) = (27/2, 273/2, 1225/8, 343/8)
   x = 294 // 2
   r2 = 147 // 2
  @printf("r1 = %.5f;  x1 = %.5f;  R = %.5f;  y = %.5f\n", r1, x1, R, y)
   θ = atan(y/x)/pi*180
   plot()
   circle(0, y, R, :blue, beginangle=-θ, endangle=180+θ, n=200)
   circle(x1, r1, r1)
   circle(-x1, r1, r1)
   circle(r2, r2, r2, :green)
   circle(-r2, r2, r2, :green)
   segment(x, 0, -x, 0, :blue)
   if more == true
       point(0, y, " y", :blue)
       point(x1, r1, "小円      \n(x1,r1)   ", :red, :right, :bottom)
       point(r2, r2, " 大円\n (r2,r2)", :green, :left, :bottom)
       point(x, 0, "x", :blue)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5, ylims=(-12, 200))
   else
       plot!(showaxis=false)
   end
end;