算額（その600）

算額（その600）

長崎市 鎮西大社諏訪神社 明治20年(1887)
米光丁: 長崎県の和算の概説
http://hyonemitsu.web.fc2.com/Nagasakiwasan.pdf

問題 13. 外円の中に大円，中円，小円を入れる。大円の直径が 36 寸のとき，小円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, r1 - R)
中円の半径と中心座標を r2, (0, R - r2)
小円の半径と中心座標を r3, x3, y3)
斜線と外円の交点座標を (x, -sqrt(R^2 - x^2))
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms x::positive, R::positive, r1::positive, r2::positive, r3::positive, x3::positive, y3::positive
eq1 = r1 + r2 - R
eq2 =x3^2 + (y3 - R + r2)^2 - (r2 - r3)^2
eq3 = (y3 - R + r2)/x3*(-R -sqrt(R^2 - x^2))/x + 1
eq4 = (r2 - 2r3)/(R - r2) - x/sqrt(x^2 +(sqrt(R^2 - x^2) + R)^2)
eq5 = distance(0, R, x, -sqrt(R^2 - x^2), x3, y3) - r3^2
eq6 = distance(0, R, x, -sqrt(R^2 - x^2), 0, r1 - R) - r1^2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6], (x, R, r2, r3, x3, y3))

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 v, r.f_converged
end;

function H(u)
   (x, R, r2, r3, x3, y3) = u
   return [
       -R + r1 + r2,  # eq1
       x3^2 - (r2 - r3)^2 + (-R + r2 + y3)^2,  # eq2
       1 + (-R - sqrt(R^2 - x^2))*(-R + r2 + y3)/(x*x3),  # eq3
       -x/sqrt(x^2 + (R + sqrt(R^2 - x^2))^2) + (r2 - 2*r3)/(R - r2),  # eq4
       -r3^2 + (x3 - x*(R^2 - R*y3 + R*sqrt(R^2 - x^2) + x*x3 - y3*sqrt(R^2 - x^2))/(2*R*(R + sqrt(R^2 - x^2))))^2 + (y3 - (R^2 + R*y3 - R*sqrt(R^2 - x^2) - x*x3 + y3*sqrt(R^2 - x^2))/(2*R))^2,  # eq5
       -r1^2 + (-x + r1*x/(2*R))^2 + (-R + r1 - (R^2 - (R + sqrt(R^2 - x^2))*(2*R - r1)/2)/R)^2,  # eq6
   ]
end;

r1 = 36//2
iniv = BigFloat[22, 35, 17, 5, 10, 20]
res = nls(H, ini=iniv)

   (BigFloat[22.62741699796952078082701958735516925711475000603116912031953235176946800128799, 36.00000000000000000000000000000000000000000000000000023997217352842064306535952, 18.0000000000000000000000000000000000000000000000000002399721735284206430653598, 6.000000000000000000000000000000000000000000000000000298124281531485670596315888, 11.3137084989847603904135097936775846285573750030155858282654603937949442060392, 22.00000000000000000000000000000000000000000000000000175488066987747432400063655], true)

小円の直径は 12 寸である。

その他のパラメータは以下の通り。
x = 22.6274;  R = 36;  r2 = 18;  r3 = 6;  x3 = 11.3137;  y3 = 22

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 36//2
   (x, R, r2, r3, x3, y3) = res[1]
   @printf("小円の直径 = %g;  x = %g;  R = %g;  r2 = %g;  r3 = %g;  x3 = %g;  y3 = %g\n", 2r3, x, R, r2, r3, x3, y3)
   plot()
   circle(0, 0, R, :blue)
   circle(0, r1 - R, r1, :green)
   circle(0, R - r2, r2, :orange)
   circle(x3, y3, r3)
   segment(0, R, x, -sqrt(R^2 - x^2))
   segment(0, R, -x, -sqrt(R^2 - x^2))
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, R - r2, " R-r2")
       point(x3, y3, "(x3,y3)")
   end
end;

算額（その599）

算額（その599）

七六 北埼玉郡騎西町中種足 雷神社 明治5年(1872)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額　明治13年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg

東京都中央区 福徳神社 令和2年(2020)
http://www.wasan.jp/tokyo/hukutoku1.html

キーワード：円4個，正三角形，斜線2本
#Julia, #SymPy, #算額, #和算

正三角形を 2 本の線で 4 つの部分に分け，その中に大小の円をそれぞれ 2 個入れる。大円の直径が 1 のとき，小円の直径はいくらか。

正三角形の一辺の長さを 2a とする。
線分と斜辺の交点座標を (±x, (a - x)√3) とする。
大円の半径と中心座標を r1, (0, y1), (0, r1)
小円の半径と中心座標を r2, (x2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, x::positive,
     r1::positive, x1::positive, y1::positive,
     r2::positive, x2::positive, y2::positive

eq1 = distance(0, sqrt(Sym(3))a, a, 0, 0, y1) - r1^2
eq2 = distance(0, sqrt(Sym(3))a, a, 0, x2, y2) - r2^2
eq3 = distance(-x, (a-x)sqrt(Sym(3)), a, 0, 0, y1) - r1^2
eq4 = distance(-x, (a-x)sqrt(Sym(3)), a, 0, x2, y2) - r2^2
eq5 = distance(x, (a-x)sqrt(Sym(3)), -a, 0, 0, r1) - r1^2
eq6 = distance(x, (a-x)sqrt(Sym(3)), -a, 0, x2, y2) - r2^2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6])

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 v, r.f_converged
end;

function H(u)
   (a, x, y1, r2, x2, y2) = u
   return [
3*a^2/4 - sqrt(3)*a*y1/2 - r1^2 + y1^2/4,  # eq1
-r2^2 + (-3*a + 3*x2 + sqrt(3)*y2)^2/16 + (-sqrt(3)*a + sqrt(3)*x2 + y2)^2/16,  # eq2
(3*a^4/4 - 3*a^3*x/2 - sqrt(3)*a^3*y1/2 - a^2*r1^2 + 3*a^2*x^2/4 + a^2*y1^2/4 + a*r1^2*x + sqrt(3)*a*x^2*y1/2 + a*x*y1^2/2 - r1^2*x^2 + x^2*y1^2/4)/(a^2 - a*x + x^2),  # eq3
-r2^2 + (a + x)^2*(-sqrt(3)*a^2 + sqrt(3)*a*x + sqrt(3)*a*x2 + a*y2 - sqrt(3)*x*x2 + x*y2)^2/(16*(a^2 - a*x + x^2)^2) + (-3*a^3 + 6*a^2*x + 3*a^2*x2 + sqrt(3)*a^2*y2 - 3*a*x^2 - 6*a*x*x2 + 3*x^2*x2 - sqrt(3)*x^2*y2)^2/(16*(a^2 - a*x + x^2)^2),  # eq4
(3*a^4 - 2*sqrt(3)*a^3*r1 - 6*a^3*x - 3*a^2*r1^2 + 3*a^2*x^2 + 6*a*r1^2*x + 2*sqrt(3)*a*r1*x^2 - 3*r1^2*x^2)/(4*(a^2 - a*x + x^2)),  # eq5
-r2^2 + (-sqrt(3)*a^3 - sqrt(3)*a^2*x2 + a^2*y2 + sqrt(3)*a*x^2 + 2*a*x*y2 + sqrt(3)*x^2*x2 + x^2*y2)^2/(16*(a^2 - a*x + x^2)^2) + (3*a^3 - 6*a^2*x + 3*a^2*x2 - sqrt(3)*a^2*y2 + 3*a*x^2 - 6*a*x*x2 + 3*x^2*x2 + sqrt(3)*x^2*y2)^2/(16*(a^2 - a*x + x^2)^2),  # eq6
   ]
end;

r1 = 1/2
iniv = BigFloat[1.6, 0.66, 1.74, 0.322, 0.56, 1.1] .* (2r1)
res = nls(H, ini=iniv)

   (BigFloat[1.573132184970985912579820468331052169654464272601005512809483077205430728737341, 0.6610365985079724914150555424853615260166663839137217848002503187911426334707547, 1.724744871391588873493812097618873490781894353477567544930381227134074172251389, 0.3224744871391588398986467839454922576529915481447599591870862696864391756755701, 0.5585421958697395716072717462463142074739205600259983862524414059477805245308117, 1.112372435695794479712354981538013592707116186365025213663076513709077532563648], true)

大円の直径が 1 のとき，小円の直径は 0.644949 である。

その他のパラメータは以下の通り。
小円の直径 = 0.644949;  r1 = 0.5;  a = 1.57313;  x = 0.661037;  y1 = 1.72474;  r2 = 0.322474;  x2 = 0.558542;  y2 = 1.11237

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1/2
   (a, x, y1, r2, x2, y2) = res[1]
   @printf("小円の直径 = %g;  r1 = %g;  a = %g;  x = %g;  y1 = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", 2r2, r1, a, x, y1, r2, x2, y2)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:magenta, lw=0.5)
   circle(0, y1, r1)
   circle(0, r1, r1)
   circle(x2, y2, r2, :blue)
   circle(-x2, y2, r2, :blue)
   segment(-a, 0, x, (a - x)√3)
   segment(a, 0, -x, (a - x)√3)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, y1, " 大円:r1\n (0,y1)", :red, :left, :vcenter)
       point(0, r1, " 大円:r1\n (0,r1)", :red, :left, :vcenter)
       point(x2, y2, "小円:r2\n(x2,y2)", :blue, :center, delta=-delta)
       point(a, 0, "a", :magenta, :left, :bottom, delta=delta/2)
       point(0, √3a, " √3a", :magenta, :left, :vcenter)
       point(x, √3(a-x), " (x,√3(a-x))", :magenta, :left, :vcenter)
   end
end;

算額（その598）

算額（その598）

八五 北葛飾郡三郷町上彦川戸 香取神社 明治13年(1880)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

東京都中央区 福徳神社 令和2年(2020)
http://www.wasan.jp/tokyo/hukutoku1.html

高さ 30，底辺 42 の直角三角形を 2 本の線分で分け，面積を 3 等分する。このとき長い方の線分の長さはいくらか。

2 本の線分と斜辺の交点座標を (x1, y1), (x2, y2)，高さを「鈎」，底辺を「股」とおく。

弦を底辺として考えると，赤，白，緑の三角形の高さは同じなので，面積がおなじになるためには弦は三等分しなければならない。
そうすれば，
x1,y1 は 42*(2/3), 30*(1/3) で 28, 10
x2,y2 は 42*(1/3), 30*(2/3) で 14, 20
と暗算でわかる。
線分の長さは暗算ではわからないが。√(28^2 + 10^2) と √(14^2 + 20^2) の長いほうが答え 29.732137494637012 である。

√(28^2 + 10^2), √(14^2 + 20^2) 

   (29.732137494637012, 24.413111231467404)

連立方程式を解くなら，以下のように。

include("julia-source.txt");

using SymPy

@syms 鈎::positive, 股::positive, x1::positive, y1::positive,
     x2::positive, y2::positive

S = 股*鈎/2
eq1 = 股*y1/2 - S/3
eq2 = 鈎*x2/2 - S/3
eq3 = y1/(股 - x1) - 鈎/股
eq4 = (鈎 - y2)/x2 - 鈎/股
(x1, y1, x2, y2) = solve([eq1, eq2, eq3, eq4], (x1, y1, x2, y2))

   Dict{Any, Any} with 4 entries:
     x1 => 2*股/3
     y1 => 鈎/3
     x2 => 股/3
     y2 => 2*鈎/3

(鈎, 股) = (30, 42)

   (30, 42)

(x1, y1) = (2*股/3, 鈎/3)

   (28.0, 10.0)

(x2, y2) = (股/3, 2*鈎/3)

   (14.0, 20.0)

(緑の面積, 赤の面積, 白の面積) = (股*y1/2, 鈎*x2/2, 股*鈎/2 - 股*y1/2 - 鈎*x2/2)

   (210.0, 210.0, 210.0)

線分の長さ =  (sqrt(x1^2 + y1^2), sqrt(x2^2 + y2^2))

   (29.732137494637012, 24.413111231467404)

長い方の線分の長さ = maximum(線分の長さ)

   29.732137494637012

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   @printf("(x1, y1) = (%g, %g);  (x2, y2) = (%g, %g)\n", x1, y1, x2, y2)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   plot!([0, 股, x1, 0],  [0, 0, y1, 0], color=:green, lw=0.5, seriestype=:shape, fillcolor=:green)
   plot!([0, x2, 0, 0], [0, y2, 鈎, 0], color=:red, lw=0.5, seriestype=:shape, fillcolor=:red)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x1, y1, "　(x1, y1)", :green, :left, :vcenter)
       point(x2, y2, "　(x2, y2)", :red, :left, :vcenter)
   end
end;