算額（その351）

算額（その351）

岐阜県大垣市赤坂町 金生山明星輪寺 元治2年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

正三角形と外円との間に，緑円，赤円，白円を計 16 個をいれる。正三角形の内接円の半径が与えられたとき白円の直径を求めよ。

記述を簡潔にできるように元の図を上下反転させて考える。
緑円の半径と中心座標 r1, (x, y)
赤円の半径と中心座標 r2, (0, 0), (0, 2r2 + 4r3), (0, 4r2 + 4r3)
白円の半径と中心座標 r3, (0, r2 + r3), (0, r2 + 3r3)
外円の半径と中心座標 r0, (0, 0);  r0 = 5r2 + 4r3
正三角形に内接する円の半径は 3r2 + 4r3

外円の半径を 1 として他の変数の値を求め，r2, r3 を求め，r2 と 3r2 + 4r3 の比を求める。

include("julia-source.txt");

using SymPy

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

r0 = 1
y = r2 + 2r3
r1 = 2r2 + 2r3
eq1 = x^2 + y^2 - (r1 + r2)^2
eq2 = x^2 + (y - r2 - r3)^2 - (r1 + r3)^2
eq3 = 5r2 + 4r3 - r0

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

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

   r0 = 1; r1 = 0.428571; r2 = 0.142857; r3 = 0.0714286; x = 0.494872, y = 0.285714
   正三角形に内接する円の半径 = 0.714286
   白円の半径 = 0.0714286

赤円，白円の半径はそれぞれ 1/7, 1/14。
正三角形の内接円の半径は 3r2 + 4r3 = 3/7 + 4/14 = 10/14。
よって，白円の半径は正三角形の内接円の半径の 1/10 である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3, x) = res[1]
   r0 = 1
   y = r2 + 2r3
   r1 = 2r2 + 2r3
   @printf("r0 = %g; r1 = %g; r2 = %g; r3 = %g; x = %g, y = %g\n", r0, r1, r2, r3, x, y)
   @printf("正三角形に内接する円の半径 = %g\n", 3r2 + 4r3)
   @printf("白円の半径 = %g\n", r3)
   plot()
   circle(0, 0, r0, :black)
   circle(0, 0, r2)
   rotate(x, y, r1, :green)
   rotate(0, r2 + r3, r3, :blue)
   rotate(0, r2 + 3r3, r3, :blue)
   rotate(0, 2r2 + 4r3, r2)
   rotate(0, 4r2 + 4r3, r2)
   ic = 3r2 + 4r3
   plot!([√3ic, -√3ic, 0, √3ic], [ic, ic, -2ic, ic], color=:black, lw=0.5)
   if more
       point(x, y, "(x,y)")
       point(0, r2 + r3, "  r2+r3", :blue, :left, :vcenter)
       point(0, r2 + 3r3, "  r2+3r3", :blue, :left, :vcenter)
       point(0, 2r2 + 4r3, " 2r2+4r3", :red, :left, :vcenter)
       point(0, 3r2 + 4r3, "   3r2+4r3", :red, :left, :bottom)
       point(0, 4r2 + 4r3, " 4r2+4r3", :red, :left, :vcenter)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その350）

算額（その350）

岐阜県大垣市赤坂町 金生山明星輪寺 元治2年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

外円内に青，緑，赤，白の 12 個の円をそれぞれ接するように入れる。外円の直径が与えられたときに白円の直径を求めよ。

2 つの青円の中心を y 軸，2 つの緑円の中心を x 軸に置く。
x 軸，y 軸で点対称なので，第1象限の図形を対象とする。

外円の半径と中心座標 r0, (0, 0)
青円の半径と中心座標 r1, (0, r1)
緑円の半径と中心座標 r2, (x2, 0)
赤円の半径と中心座標 r3, (x3, y3)
白円の半径と中心座標 r4, (x4, y4)

以下の連立方程式を nlsolve() で解く。

include("julia-source.txt");

using SymPy

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

r0 = 1
r1 = r0//2
eq1 = x2^2 + (r0 - r1)^2 - (r1 + r2)^2 |> expand
eq2 = x3^2 + (r0 - r1 - y3)^2 - (r1 + r3)^2 |> expand
eq3 = (x3 - x2)^2 + y3^2 - (r2 + r3)^2 |> expand
eq4 = (x4 - x2)^2 + r4^2 - (r2 + r4)^2 |> expand
eq5 = (x4 - x3)^2 + (y3 - r4)^2 - (r3 + r4)^2 |> expand
eq6 = x3^2 + y3^2 - (r0 - r3)^2 |> expand
eq7 = x4^2 + r4^2 - (r0 - r4)^2 |> expand;

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7]) #, (r2, x2, r3, x3, y3, r4, x4))

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

iniv = [big"0.5", 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]
res = nls(H, ini=iniv);
println(res);

 (BigFloat[0.2307692307692307945960496972420115643146885404357311695540829279753122088871492, 0.5329387100211930608365475760066592887331159643976213583701227822622771799660959, 0.1999999999999999665985768137693666773172486596065619280957172091107786570697997, 0.6928203230275509687226737223425773030846798934032123402136212243156691574043034, 0.400000000000000100204269558691899968048254021180314215712850083601413507819824, 0.1249999999999999591818515058195745225114059751553107390846466073349432317897547, 0.8660254037844387219680669190319846379474904384664942387376085192902656644536619], true)

外円の半径が 1 のとき，緑円の半径 = 0.230769;  赤円の半径 = 0.2;  白円の半径 = 0.125

0.230769 = 3/13

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 1
   r1 = 1//2
   (r2, x2, r3, x3, y3, r4, x4) = res[1]
   @printf("外円の半径が %g のとき，緑円の半径 = %g;  赤円の半径 = %g;  白円の半径 = %g\n", r0, r2, r3, r4)
   plot()
   circle(0, 0, r0, :black)
   circle(0, r0 - r1, r1, :blue)
   circle(0, r1 - r0, r1, :blue)
   circle(x2, 0, r2, :green)
   circle(-x2, 0, r2, :green)
   circle4(x3, y3, r3)
   circle4(x4, r4, r4, :gray)
   if more
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;