算額（その726）

算額（その726）

群馬県桐生市梅田町 日枝神社 明治12年
https://gunmawasan.web.fc2.com/files/sangak-corner/hiejin.htm

外円内に弦を引き，菱形と天円と地円を入れる。地円，天円の直径がそれぞれ 7 寸，2 寸のとき，外円の直径はいかほどか。

菱形の対角線の長い方と短い方の長さをそれぞれ 2a, 2b とする。
外円の半径と中心座標を R, (0, 0)
天円の半径と中心座標を r1, (x1, y1)
地円の半径と中心座標を r2, (a, R - 2b + r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, a::positive, b::positive, d,
     r1::positive, x1::positive, y1::positive, r2::positive
eq1 = x1^2 + y1^2 - (R - r1)^2
eq2 = dist(0, R, a, R - b, a, R - 2b + r2) - r2^2
eq2 = div(numerator(apart(eq2, d)), b)  # 分子だけを評価する
eq3 = 4(R^2 - (R - 2r1)^2) - (a^2 + b^2)
eq4 = a^2 + (R - b)^2 - R^2
eq5 = b/a - x1/y1;
res = solve([eq1, eq2, eq3, eq4, eq5], (R, a, b, x1, y1));

12 組の解が得られるが，すべてのパラメータの符号を考慮すると 7 番目のものが適解である。

res[7]   

   (8*r1^2/(4*r1 - r2), r2*sqrt(4*r1 + r2)*sqrt(1/(4*r1 - r2)), 4*r1 + r2, (r1 + r2/4)*sqrt(16*r1^2 - r2^2)/(4*r1 - r2), r2*(4*r1 + r2)/(4*(4*r1 - r2)))

外円の半径は「天円の半径を二乗したものを，天円半径の4倍から地円の半径を引いたもので割り，8 倍する」ことで得られる。
天円，地円の半径がそれぞれ 1 寸, 3.5 寸 のとき，外円の半径は 16 寸である。
天円，地円の直径がそれぞれ 2 寸, 7 寸 のとき，外円の直径は 32 寸である。

res[7][1](r1 => 2//2, r2 => 7//2) |> println

   16

その他のパラメータは以下のとおりである。

R = 16;  a = 13.5554;  b = 7.5;  x1 = 7.26184;  y1 = 13.125

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (2, 7) .// 2
   t = 4*r1 - r2
   u = 4*r1 + r2
   (R, a, b, x1, y1) = (8r1^2/t, r2*sqrt(u/t), u, u*sqrt(t*u)/4t, r2*u/(4t))
   @printf("天円の直径 = %g，地円の直径 = %g のとき，外円の直径 = %g\n", r1, r2, 2R)
   @printf("R = %g;  a = %g;  b = %g;  x1 = %g;  y1 = %g\n", R, a, b, x1, y1)
   plot([a, 0, -a, 0, a], [R - b, R, R - b, R - 2b, R - b], color=:green, lw=0.25)
   circle(0, 0, R, :orange)
   circle(x1, y1, r1)
   circle(-x1, y1, r1)
   circle(a, R - 2b + r2, r2, :blue, beginangle=90, endangle=270)
   circle(-a, R - 2b + r2, r2, :blue, beginangle=270, endangle=450)
   y0 = R - 2b
   x0 = sqrt(R^2 - y0^2)
   segment(-x0, y0, x0, y0, :magenta)
   segment(a, y0, a, y0+2r2, :blue)
   segment(-a, y0, -a, y0+2r2, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
   end
end;

算額（その725）

算額（その725）

群馬県桐生市梅田町 日枝神社 明治12年
https://gunmawasan.web.fc2.com/files/sangak-corner/hiejin.htm

外円内に弦と斜線を 2 本引き，区画された領域に等円を 4 個入れる。
外円の直径が与えられたとき，等円の直径を求める方法を問う。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x, R - 3r), (0, r - R), (0, R - r)
斜線と円の交点座標を (x0, y0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, x::positive,
     x0::positive, y0::negative, d
eq0 = x0^2 + y0^2 - R^2
eq1 = x^2 + (R - 3r)^2 - (R - r)^2
eq2 = dist(0, R - 2r, x0, -sqrt(R^2 - x0^2), x, R - 3r) - r^2
eq2 = dist(0, R - 2r, x0, y0, x, R - 3r) - r^2
eq2 = numerator(apart(eq2, d)) |> simplify  # 分子 = 0 で十分
eq3 = dist(0, R - 2r, x0, -sqrt(R^2 - x0^2), 0, r - R) - r^2;
eq3 = dist(0, R - 2r, x0, y0, 0, r - R) - r^2;
eq3 = numerator(apart(eq3, d)) |> simplify  # 分子 = 0 で十分
res = solve([eq0, eq1, eq2, eq3], (r, x, x0, y0))

   1-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (-sqrt(17)*R/34 + R/2, sqrt(2)*R*sqrt(-53 + 13*sqrt(17))/4 + 5*sqrt(34)*R*sqrt(-53 + 13*sqrt(17))/68, R*sqrt(-424/17 + 104*sqrt(17)/17), R*(-4 + 13*sqrt(17)/17))

等円の直径は外円の直径の (17 - √17)/34 倍である。

「術」では，「五を開平して...」と言っているので，間違いであろう。

res[1][1] |> simplify |> println

   R*(17 - sqrt(17))/34

その他のパラメータは以下のとおりである。

   r = 46.5841;  x = 74.5571;  x0 = 65.3787;  y0 = -104.186

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 123
   (r, x, x0, y0) = (R*(17 - sqrt(17))/34, sqrt(2)*R*sqrt(-53 + 13*sqrt(17))/4 + 5*sqrt(34)*R*sqrt(-53 + 13*sqrt(17))/68, R*sqrt(-424/17 + 104*sqrt(17)/17), R*(-4 + 13*sqrt(17)/17))
   @printf("r = %g;  x = %g;  x0 = %g;  y0 = %g\n", r, x, x0, y0)
   @show(R*(17 - sqrt(17))/34)
   plot()
   circle(0, 0, R, :blue)
   circle(0, R - r, r)
   circle(x, R - 3r, r)
   circle(-x, R - 3r, r)
   circle(0, r - R, r)
   y1 = R - 2r
   x1 = sqrt(R^2 - y1^2)
   segment(-x1, y1, x1, y1, :green)
   segment(0, y1, x0, y0, :green)
   segment(0, y1, -x0, y0, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(x0, y0, "(x0,y0)", :green, :left, delta=-delta/2)
       point(x, R - 3r, "等円:r,(x,R-3r)", :red, :center, delta=-delta)
       point(0, R - r, " R-r", :red, :left, :vcenter)
       point(0, R - 2r, "R-2r", :red, :center, :bottom, delta=delta/2)
       point(0, r - R, " r-R", :red, :left, :vcenter)
   end
end;

算額（その724）

算額（その724）

群馬県桐生市梅田町 日枝神社 明治12年
https://gunmawasan.web.fc2.com/files/sangak-corner/hiejin.htm

長方形と輪（円）が交わってできる領域に 3 個の小円がある。
また，下部に 2 ほんの斜線を入れ，正方形を 3 個入れる。左右の正方形の 2 つの頂点は輪の円周上と斜線の上にある。中央の正方形の上辺は輪と接している。
小円の著系が 1 寸のとき，正方形の一辺の長さはいかほどか。

長方形の短辺と長辺を 2a, b とする。
輪の半径と中心座標を r1, (0, 2c + r1); a = r1
小円の半径と中心座標を r2, (a - r2, b - r2)
正方形の一辺の長さを 2c
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive, r1::positive, r2::positive, y0::positive
a = r1
r2 = 1//2
eq1 = 2c + 2r1 - 2r2 - b
eq2 = y0/2c - 2c/(a - c)
eq3 = (a - 2c)^2 + (y0 + 2c - 2c - r1)^2 - r1^2
eq4 = ( a - r2)^2 + (b - r2 - 2c - r1)^2 - (r1 + r2)^2;
res = solve([eq1, eq2, eq3, eq4], (b, c, r1, y0))

   4-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (2*r2/5, r2/5, r2, r2/5)
    (r2, r2/2, r2, 2*r2)
    (98*r2/5, 9*r2/5, 9*r2, 9*r2/5)
    (25*r2, 9*r2/2, 9*r2, 18*r2)

4 組の解が得られるが，3 番目のものが適解である。

r2 = 1/2 のとき，
(b, c, r1, y0) =(98*r2/5, 9*r2/5, 9*r2, 9*r2/5) = (9.8, 0.9, 4.5, 0.9) である。
すなわち，正方形の一辺の長さは小円の直径の 1.8 倍，すなわち 1 寸 8 分である。

r2 = 1/2
(98*r2/5, 9*r2/5, 9*r2, 9*r2/5)

   (9.8, 0.9, 4.5, 0.9)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1//2
   (b, c, r1, y0) = (98*r2/5, 9*r2/5, 9*r2, 9*r2/5)
   @printf("正方形の一辺の長さ = %g\n", 2c)
   @printf("b = %g;  c = %g;  r1 = %g;  y0 = %g\n", b, c, r1, y0)
   a = r1
   plot([a, a, -a, -a, a], [0, b, b, 0, 0], color=:black, lw=0.5)
   circle(0, 2c + r1, r1)
   circle(0, b + r2, r2, :blue)
   circle(a - r2, b - r2, r2, :blue)
   circle(r2 - a, b - r2, r2, :blue)
   rect(-c, 0, c, 2c, :green)
   rect(a - 2c, y0, a, y0 + 2c, :green) 
   rect(2c - a, y0, -a, y0 + 2c, :green)
   segment(a, 0, c, 2c)
   segment(-a, 0, -c, 2c)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(a - 2c, y0 + 2c, " (a-2c,y0+2c)", delta=-delta/2)
       point(a - 2c, y0, " (a-2c,y0)", :green, :left, :bottom, delta=delta/2)
       point(c, 0, " c", :green, :left, :bottom, delta=delta/2)
       point(0, 2c, " 2c", :green, :left, :bottom, delta=delta/2)
       point(0, 2c + r1, " 輪:r1,(0,2c+r1)", :red, :left, :vcenter)
       point(a - r2, b - r2, "小円:r2,(a-r2,b-r2)    ", :blue, :right, :vcenter)
       point(0, b + r2, " b+r2", :black, :center, :bottom, delta=delta/2)
       point(a, b, "(a,b)", :black, :center, :bottom, delta=delta/2)
   end
end;

算額（その723）

算額（その723）

群馬県桐生市天神町 桐生天満宮 明治11年
https://gunmawasan.web.fc2.com/files/sangak-corner/tenman-g1.htm

外円内に甲円，乙円，丙円，丁円，戊円が互いに接して入っている。乙円，丙円，戊円の直径がそれぞれ 15寸，6 寸，4 寸のとき，外円，甲円，丁円の直径はいかほどか。

算額（その127）を解き直した。
https://blog.goo.ne.jp/r-de-r/e/401c8ab1043169fd4660aec699096fc4

パラメータを減らすために，甲円の中心が y 軸上にあるように回転して考える。
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (x4, y4)
戊円の半径と中心座標を r5, (x5, y5)
とおき，以下の連立方程式を解く。
なお，丙円と丁円の外接条件として eq6 を用いようとしたが，感度不足で思うような収束解が得られなかった。
算法助術に活路を見出そうとして，公式77から eq12（または eq13)を得て，やっと解が得られた。

今回の算額もまた，図と実際の解は全く異なるものである。最も，問と答を見れば，図における円の比率がおかしいのはすぐに分かるのではあるが。

include("julia-source.txt");

using SymPy
@syms R, r1, r2, x2, y2, r3, x3, y3, r4, x4, y4, r5, x5, y5, t
eq1 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
eq2 = x3^2 + (R - r1 - y3)^2 - (r1 + r3)^2
eq3 = x5^2 + (R - r1 - y5)^2 - (r1 + r5)^2
eq4 = (x2 - x4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq5 = (x2 - x5)^2 + (y2 - y5)^2 - (r2 + r5)^2
# eq6 = (x3 - x4)^2 - (y3 - y4)^2 - (r3 + r4)^2  # 感度不足
eq7 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2
eq8 = (x4 - x5)^2 + (y4 - y5)^2 - (r4 + r5)^2
eq9 = x2^2 + y2^2 - (R - r2)^2
eq10 = x3^2 + y3^2 - (R - r3)^2
eq11 = x4^2 + y4^2 - (R - r4)^2;
# 算法助術 公式77 より　　eq12 か eq13 いずれか
eq12 = r4*r5*t^2 + r1*r5*t^2 + 2r1*r4*t^2 - 2r1*r4*r5*(r2 + r3) -4r2*r1*r3*r4;
# eq13 = -R*r4*t^2 - R*r1*t^2 + 2r1*r4*t^2 +2R*r1*r4*(r2 + r3) - 4r2*r1*r3*r4;

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

function H(u)
   (R, r1, x2, y2, x3, y3, r4, x4, y4, x5, y5) = u
   return [
       x2^2 - (r1 + r2)^2 + (R - r1 - y2)^2,  # eq1
       x3^2 - (r1 + r3)^2 + (R - r1 - y3)^2,  # eq2
       x5^2 - (r1 + r5)^2 + (R - r1 - y5)^2,  # eq3
       -(r2 + r4)^2 + (x2 - x4)^2 + (y2 - y4)^2,  # eq4
       -(r2 + r5)^2 + (x2 - x5)^2 + (y2 - y5)^2,  # eq5
       -(r3 + r5)^2 + (x3 - x5)^2 + (y3 - y5)^2,  # eq7
       -(r4 + r5)^2 + (x4 - x5)^2 + (y4 - y5)^2,  # eq8
       x2^2 + y2^2 - (R - r2)^2,  # eq9
       x3^2 + y3^2 - (R - r3)^2,  # eq10
       x4^2 + y4^2 - (R - r4)^2,  # eq11
       -4*r1*r2*r3*r4 - 2*r1*r4*r5*(r2 + r3) + 2*r1*r4*t^2 + r1*r5*t^2 + r4*r5*t^2,  # eq12
   ]
end;

(r2, r3, r5) = (15, 6, 4) .// 2
t = sqrt(Sym(2)r2*r3)
iniv = BigFloat[31, 15.5, 21, 7.4, 17, 21, 2.7, 21.5, 17.5, 17, 16]
res = nls(H, ini=iniv)

   ([30.0, 15.0, 21.213203435596427, 7.5, 16.97056274847714, 21.0, 2.5, 21.213203435596427, 17.5, 16.97056274847714, 16.0], true)

外円の直径 = 60;  甲円の直径 = 30;  丁円の直径 = 5 である。
なお，乙円はx 軸に接する。すなわち，y2 = r2 であった。

その他のパラメータは以下のとおりである。

r2 = 7.5;  r3 = 3;  r5 = 2
外円の直径 = 60;  甲円の直径 = 30;  丁円の直径 = 5
R = 30;  r1 = 15;  x2 = 21.2132;  y2 = 7.5;  x3 = 16.9706;  y3 = 21;  r4 = 2.5;  x4 = 21.2132;  y4 = 17.5;  x5 = 16.9706;  y5 = 16

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3, r5) = (15, 6, 4) .// 2
   (R, r1, x2, y2, x3, y3, r4, x4, y4, x5, y5) = res[1]
   @printf("r2 = %g;  r3 = %g;  r5 = %g\n", r2, r3, r5)
   @printf("外円の直径 = %g;  甲円の直径 = %g;  丁円の直径 = %g\n", 2R, 2r1, 2r4)
   @printf("R = %g;  r1 = %g;  x2 = %g;  y2 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  x4 = %g;  y4 = %g;  x5 = %g;  y5 = %g\n", R, r1, x2, y2, x3, y3, r4, x4, y4, x5, y5)
   plot()
   circle(0, 0, R, :orange)
   circle(0, R - r1, r1)
   circle(x2, y2, r2, :blue)
   circle(x3, y3, r3, :magenta)
   circle(x4, y4, r4, :green)
   circle(x5, y5, r5, :brown)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, R - r1, "甲円:r1,(0,R-r1)", :red, :center, delta=-delta/2)
       point(x2, y2, "乙円:r2,(x2,y2)", :blue, :center, delta=-delta/2)
       point(x3, y3, "丙円:r3,(x3,y3) ", :black, :right, :vcenter)
       point(x4, y4, "丁円:r4,(x4,y4) ", :black, :right, :bottom, delta=delta/2)
       point(x5, y5, "戊円:r5,(x5,y5) ", :black, :right, :vcenter)
   end
end;