算額（その299）

算額（その299）

愛媛県 伊佐爾波神社 明治6年
高阪金次郎峻則
http://www.wasan.jp/ehime/isaniwa22.html
https://www.city.matsuyama.ehime.jp/kanko/kankoguide/rekishibunka/bunkazai/ken/isaniwa_sangaku.html
https://isaniwa.official.jp/2012/12/28/%E7%AE%97%E9%A1%8D%E3%81%AB%E6%8C%91%E6%88%A6/

扇面（中心角 120 度）に東円 1 個，西円，南円，北円が 2 個ずつ入っている。南円の直径を知って北円の直径を求めよ。

以下の方程式を解けばよいが，一挙に解くのが難しい。

方程式は [eq4, eq7] と残りの方程式は(r0 を除いて）互いに独立である。

まず [eq4, eq7] を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive, r3::positive, r4::positive,
     x2::positive, x3::positive, x4::positive, y4::positive;

r1 = r0/4
eq1 = x2^2 + (r0 - r1 - r0/2 - r2)^2 - (r1 + r2)^2
eq2 = x4^2 + (y4 - r0 + r1)^2 - (r1 + r4)^2
eq3 = (x2 - x4)^2 + (y4 - r0/2 - r2)^2 - (r2 + r4)^2
eq4 = x3^2 + (r0/2 - r3)^2 - (r3 + r0//2)^2
eq5 = x2^2 + (r0/2 + r2)^2 - (r0 - r2)^2
eq6 = x4^2 + y4^2 - (r0 - r4)^2
eq7 = r3/(sqrt(Sym(3))*r0/2 - x3) - tan(PI/12);

まず，eq4, eq7 は 未知数 r3, x3 を r0 で表して解くことができる。

res = solve([eq4, eq7], (r3, x3))

   1-element Vector{Tuple{Sym, Sym}}:
    (3*r0*(7 - 4*sqrt(3))/2, r0*(-3 + 2*sqrt(3)))

次に，残りの方程式のセット [eq1, eq2, eq3, eq5, eq6] から (r2, x2, r4, x4, y4) を求める（r0 を記号として含む）。

東円の半径は扇の半径の 1/4 である。

@syms r0::positive, r1::positive, r2::positive, r3::positive, r4::positive,
     x2::positive, x3::positive, x4::positive, y4::positive;

r1 = r0//4
eq1 = x2^2 + (r0 - r1 - r0//2 - r2)^2 - (r1 + r2)^2
eq2 = x4^2 + (y4 - r0 + r1)^2 - (r1 + r4)^2
eq3 = (x2 - x4)^2 + (y4 - r0//2 - r2)^2 - (r2 + r4)^2
eq5 = x2^2 + (r0//2 + r2)^2 - (r0 - r2)^2
eq6 = x4^2 + y4^2 - (r0 - r4)^2;

res = solve([eq1, eq2, eq3, eq5, eq6], (r2, x2, r4, x4, y4))

   1-element Vector{NTuple{5, Sym}}:
    (3*r0/16, r0*(-3 + 14*sqrt(3))*sqrt(28*sqrt(3) + 199)/772, 3*r0*(25 - 12*sqrt(3))/193, r0*sqrt(336*sqrt(3)/37249 + 2388/37249), r0*(68/193 + 60*sqrt(3)/193))

res[1][1] |> simplify |> println  # r2 西円の半径

   3*r0/16

res[1][3] |> simplify |> println  # r4 北円の半径

   3*r0*(25 - 12*sqrt(3))/193

以上の結果から，問題への答えは，「北円の半径は南円の半径の (32√3 + 62)/193 ≒ 0.6084229318 倍」である。

(3*r0*(25 - 12*sqrt(Sym(3)))/193) / (3*r0*(7 - 4*sqrt(Sym(3)))/2) |> simplify |> println

   32*sqrt(3)/193 + 62/193

扇形の半径が 20 のとき，各円の直径を求める。

   南円の直径 = 4.307806
   北円の直径 = 2.620968
   東円の直径 = 10.000000
   西円の直径 = 7.500000
   扇形の半径 = 20.000000
   北円の直径は南円の直径の 0.6084229318 倍

一挙に数値解を求めることもできる。

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)
   (r0, r2, x2, x3, r4, x4, y4) = u
   return [
x2^2 + (r0/4 - r2)^2 - (r0/4 + r2)^2,  # eq1
x4^2 + (-3*r0/4 + y4)^2 - (r0/4 + r4)^2,  # eq2
-(r2 + r4)^2 + (x2 - x4)^2 + (-r0/2 - r2 + y4)^2,  # eq3
x3^2 + (r0/2 - r3)^2 - (r0/2 + r3)^2,  # eq4
x2^2 + (r0/2 + r2)^2 - (r0 - r2)^2,  # eq5
x4^2 + y4^2 - (r0 - r4)^2,  # eq6
r3/(sqrt(3)*r0/2 - x3) - 2 + sqrt(3),  # eq7
   ]
end;

r3 = 1
iniv = [big"157.0", 29.1, 68.2, 73.2, 10.3, 45.2, 140]
iniv = [big"9.0", 2, 4, 4, 0.5, 3, 8]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[9.285468820183671312972329900930808394772674258716870318199501414526902889954248, 1.74102540378443837117861943355628490397136248947194920156030216985869867358414, 4.020725942163689539353471755593597407534879935141037673297496272009488982864332, 4.309401076758502716985195569314300770820850779232595561630515863568051039690857, 0.6084229318248914707848103579780154432784296686583510598363298988001496181698688, 2.621938437530752623591216631288396179019582647806301923444247679455871388075398, 8.271430600475518861666281929830511546667832242439599342793407445063645710020953], true)

@printf("r3 = %.6f\n", r3)
names = ("r0", "r2", "x2", "x3", "r4", "x4", "y4")
for (i, name) in enumerate(names)
   @printf("%2s = %.6f\n", name, res[1][i])
end

   r3 = 1.000000
   r0 = 9.285469
   r2 = 1.741025
   x2 = 4.020726
   x3 = 4.309401
   r4 = 0.608423
   x4 = 2.621938
   y4 = 8.271431

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 20
   r1 = r0/4
   (r3, x3) = (3*r0*(7 - 4*sqrt(3))/2, r0*(-3 + 2*sqrt(3)))
   (r2, x2, r4, x4, y4) = (3*r0/16, r0*(-3 + 14*sqrt(3))*sqrt(28*sqrt(3) + 199)/772, 3*r0*(25 - 12*sqrt(3))/193, r0*sqrt(336*sqrt(3)/37249 + 2388/37249), r0*(68/193 + 60*sqrt(3)/193))
   @printf("南円の直径 = %.6f\n", 2r3)
   @printf("北円の直径 = %.6f\n", 2r4)
   @printf("東円の直径 = %.6f\n", 2r1)
   @printf("西円の直径 = %.6f\n", 2r2)
   @printf("扇形の半径 = %.6f\n", r0)
   @printf("北円の直径は南円の直径の %.10f 倍\n", (32√3 + 62)/193)
   plot()
   circle(0, 0, r0, :black, beginangle=30, endangle=150)
   circle(0, 0, r0/2, :black, beginangle=30, endangle=150)
   circle(x3, r0/2 - r3, r3, :blue)
   circle(-x3, r0/2 - r3, r3, :blue)
   circle(0, r0 - r1, r1)
   circle(x2, r0/2 + r2, r2, :green)
   circle(-x2, r0/2 + r2, r2, :green)
   circle(x4, y4, r4, :orange)
   circle(-x4, y4, r4, :orange)
   segment(-√3r0/2, r0/2, √3r0/2, r0/2)
   segment(-√3r0/2, r0/2, 0, 0)
   segment(√3r0/2, r0/2, 0, 0)
   if more
       point(0, r0 - r1, " 東円:r1\n r0-r1", :red)
       point(x2, r0/2 + r2, " 西円:r2\n(x2,r0/2+r2)", :green, :center)
       point(x3, r0/2 - r3, "   南円:r3\n(x3,r0/2-r3)", :blue, :left, :bottom)
       point(x4, y4, "  北円:r4,(x4,y4)\n", :orange, :left, :bottom)
       point(0, r0/2, " r0/2", :black)
       point(0, r0, " r0", :black, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;