算額（その585）

算額（その585）

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

問題 19. 外円の中に長方形，その中に 4 個の円を入れる。甲円，乙円の直径が 15 寸，11 寸のとき，外円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
長方形の第1象限にある頂点の座標を (x, r1); x = sqrt(R^2 - r1^2)
甲円の半径と中心座標を r1, (x - r1, r1)
乙円の半径と中心座標を r2, (0, r1 - r2)
とおき，以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, r2::positive, x::positive

eq = (sqrt(R^2 -r1^2) - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
res = solve(eq, R)
res |> println

   Sym[sqrt(-4*r1^(3/2)*sqrt(r2) + 2*r1^2 + 4*r1*r2), sqrt(4*r1^(3/2)*sqrt(r2) + 2*r1^2 + 4*r1*r2)]

解は 2 個得られるが，2 番目のものが適解である。

甲円，乙円の直径が 15 寸，11 寸のとき，外円の半径は 21.6835995323649，直径は 43.3671990647299 である。

res[2](r1 => 15/2, r2 => 11/2) |> println
2res[2](r1 => 15/2, r2 => 11/2) |> println

   21.6835995323649
   43.3671990647299

算額の答えは「外円の直径は四十三寸五八」となっているが何らかの誤りであろう。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (15, 11) .// 2
   R = sqrt(4*r1^(3/2)*sqrt(r2) + 2*r1^2 + 4*r1*r2)
   x = sqrt(R^2 - r1^2)
   @printf("外円の直径 = %g;  R = %g; r1 = %g;  r2 = %g;  x = %g\n", 2R, R, r1, r2, x)
   plot([x, x, -x, -x, x], [-r1, r1, r1, -r1, -r1], color=:blue, lw=0.5)
   circle(0, 0, R)
   circle(x - r1, 0, r1, :magenta)
   circle(r1 - x, 0, r1, :magenta)
   circle(0, r1 - r2, r2, :green)
   circle(0, r2 - r1, r2, :green)
   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(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(x, r1, "(x,r1)  ", :red, :right, :bottom, delta=delta)
       point(x - r1, 0, "甲円:r1,(x-r1,0)", :magenta, :center, :bottom, delta=delta)
       point(0, r1 - r2, "乙円:r2,(0,r1-r2) ", :black, :right, :vcenter)
   end
end;

算額（その584）

算額（その584）

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

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

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, R - r1)
中円の半径と中心座標を r2, (0, r2 + r3)
小円の半径と中心座標を r3, (0, 0), (0, 2r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, r2::positive, r3::positive
eq1 = R - 2r1 +r3
eq2 = (R - r1)^2 + (2r2)^2 - (r1 + r3)^2
res = solve([eq1, eq2], (R, r1))

   1-element Vector{Tuple{Sym, Sym}}:
    ((2*r2^2 - r3^2)/r3, r2^2/r3)

外円の半径は (2r2^2 - r3^2)/r3
直径は中円，小円の直径を R2, R3 とすると 2R2^2/R3 - R3

中円，小円の直径が 10 寸，5 寸のとき，外円の直径は 35 寸。

ちなみに，大円の直径は 20 寸である。

@syms R2, R3
r2 = R2/2
r3 = R3/2
2(2r2^2 - r3^2)/r3 |> simplify |> println

   2*R2^2/R3 - R3

R2 = 10; R3 = 5; 2R2^2/R3 - R3

   35.0

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3) = (10, 5) .// 2
   (R, r1) = (35//2,10)
   @printf("外円の直径 = %g;  R = %g; r1 = %g;  r2 = %g;  r3 = %g\n", 2R, R, r1, r2, r3)
   plot()
   circle(0, 0, R)
   circle(R - r1, 0, r1, :magenta)
   circle(r1 - R, 0, r1, :magenta)
   circle(0, r3 + r2, r2, :green)
   circle(0, -r3 - r2, r2, :green)
   circle(0, 0, r3, :blue)
   circle(0, 2r2, r3, :blue)
   circle(0, -2r2, r3, :blue)
   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(R - r1, 0, "大円:r1,(R-r1,0)", :magenta, :center, delta=-delta)
       point(0, r3 + r2, "  中円:r2,(0,r3+r2)", :green, :left, :vcenter)
       point(0, 2r2, "  小円:r3,(0,2r2)", :blue, :left, :vcenter)
       point(R, 0, "R ", :red, :right, :bottom, delta=delta/2)
       point(r3, 0, " r3", :blue, :left, :bottom, delta=delta/2)
   end
end;

算額（その583）

算額（その583）

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

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

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, r1 - R)
小円の半径と中心座標を r2, y + r2)
ただし，y は弦と y 軸の交点の y 座標値
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, r2::positive, y::positive
eq1 = 2r1 - y - R
eq2 = r2^2 + (y + r2)^2 - (R - r2)^2
res = solve([eq1, eq2], (R, y))

   1-element Vector{Tuple{Sym, Sym}}:
    ((2*r1 + r2)^2/(4*r1), (4*r1^2 - 4*r1*r2 - r2^2)/(4*r1))

外円の直径は (2r1 + r2)^2/(2r1) = (大円の直径 + 小円の半径)^2/大円の直径 である。
大円，小円の直径が 36 寸，24 寸のとき，64 寸である。

(36 + 12)^2/36

   64.0

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (36, 24) .// 2
   (R, y) = ((2r1 + r2)^2/4r1, r1 - r2 - r2^2/4r1)
   @printf("外円の直径 = %g;  R = %g;  y = %g;  r1 = %g;  r2 = %g\n", 2R, R, y, r1, r2)
   plot()
   circle(0, 0, R)
   circle(0, r1 - R, r1, :magenta)
   circle(r2, y + r2, r2, :blue)
   circle(-r2, y + r2, r2, :blue)
   x = sqrt(R^2 - y^2)
   segment(-x, y, x, y, :green)
   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, y, " y", :green, :left, :bottom, delta=delta)
       point(r2, y + r2, "小円:r2,(r2,y+r2)", :blue, :center, :top, delta=-delta)
       point(0, r1 - R, " 大円:r1,(0,r1-R)", :magenta, :left, :vcenter)
       point(R, 0, "R ", :red, :right, :bottom, delta=delta)
   end
end;

算額（その582）

算額（その582）

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

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

算額（その334）で小生が最初に考えたものと図形的には同じである。
これも算額にはよくあることであるが，算額の図と実際に得られた解に基づいて描いた図の乖離に驚く。

外円の半径と中心座標を r, (0, 0)
小円の半径と中心座標を r1, (r1, y1)
中円の半径と中心座標を r2, (r2, y2)
とおき，以下の連立方程式を解く。

最初，y2 を negative としていたので，解が得られない。実に，y2 は 0 である。
図を描いてみれば，問題文を見ただけで解けたのにと悔しくなる。

include("julia-source.txt");

using SymPy

@syms r::positive, r1::positive, y1::positive, r2::positive, y2;

(r1, r2) = (20, 40) .// 2
eq1 = (r1 - r2)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq2 = r1^2 + y1^2 - (r - r1)^2
eq3 = r2^2 + y2^2 - (r - r2)^2

res1 = solve([eq1, eq2, eq3], (r, y1, y2))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (40, 20*sqrt(2), 0)

外円の直径は中円の直径の2倍，80 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (20//2, 40//2)
   (r, y1, y2) = (40, 20*sqrt(2), 0)
   @printf("外円の直径 = %g\n", 2r)
   plot()
   circle(0, 0, r)
   circle(r1, y1, r1, :blue)
   circle(-r1, y1, r1, :blue)
   circle(r2, y2, r2, :orange)
   circle(-r2, y2, r2, :orange)
   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(r1, y1, "小円:r1\n(r1,y1)", :blue, :center, :top, delta=-delta/2)
       point(r2, y2, "中円:r2,(r2,y2)", :black, :center, :top, delta=-delta/2)
       point(r, 0, "r ", :red, :right)
   end
end;

算額（その581）

算額（その581）

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

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

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (x1, y - r1)
小円の半径と中心座標を r2,（x2, y - r2）, (0, y + r2)
弦と y 軸の交点座標を (0, y)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, y::positive, r1::positive, x1::negative,
     r2::positive, x2::positive

eq1 = (x2 - x1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x1^2 + (y - r1)^2 - (R - r1)^2
eq3 = x2^2 + (y - r2)^2 - (R - r2)^2
eq4 = y + 2r2 - R
res = solve([eq1, eq2, eq3, eq4], (y, x1, r2, x2))

   1-element Vector{NTuple{4, Sym}}:
    ((-4*R^(11/2)*sqrt(r1) + 32*R^(9/2)*r1^(3/2) + 128*R^(7/2)*r1^(5/2) - 2048*R^(5/2)*r1^(7/2) + 7168*R^(3/2)*r1^(9/2) - 8192*sqrt(R)*r1^(11/2) + R^6 - 4*R^5*r1 - 64*R^4*r1^2 + 384*R^3*r1^3 + 256*R^2*r1^4 - 5120*R*r1^5 + 8192*r1^6)/(R^5 - 8*R^4*r1 - 32*R^3*r1^2 + 512*R^2*r1^3 - 1792*R*r1^4 + 2048*r1^5), 2*sqrt(2)*(-R^(11/2)*sqrt(r1) - R^(9/2)*r1^(3/2) + 68*R^(7/2)*r1^(5/2) - 80*R^(5/2)*r1^(7/2) - 1088*R^(3/2)*r1^(9/2) + 2560*sqrt(R)*r1^(11/2) + 4*R^5*r1 - 14*R^4*r1^2 - 200*R^3*r1^3 + 1184*R^2*r1^4 - 1408*R*r1^5 - 1024*r1^6)/(sqrt(R^(3/2)*sqrt(r1) + 8*sqrt(R)*r1^(3/2) - 5*R*r1 - 4*r1^2)*(-R^4 + 4*R^3*r1 + 48*R^2*r1^2 - 320*R*r1^3 + 512*r1^4)), 2*(R^(11/2)*sqrt(r1) - 8*R^(9/2)*r1^(3/2) - 32*R^(7/2)*r1^(5/2) + 512*R^(5/2)*r1^(7/2) - 1792*R^(3/2)*r1^(9/2) + 2048*sqrt(R)*r1^(11/2) - R^5*r1 + 8*R^4*r1^2 + 32*R^3*r1^3 - 512*R^2*r1^4 + 1792*R*r1^5 - 2048*r1^6)/(R^5 - 8*R^4*r1 - 32*R^3*r1^2 + 512*R^2*r1^3 - 1792*R*r1^4 + 2048*r1^5), sqrt(8*R^(3/2)*sqrt(r1) + 64*sqrt(R)*r1^(3/2) - 40*R*r1 - 32*r1^2))

小円の半径を表す式は長いが，問で与えられた数数値を代入すると，小円の半径は6，直径は 12 寸である。

res[1][3](R => 32//2, r1 => 18//2) |> println

   6

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

R = 16;  r1 = 9;  x1 = -4.89898;  r2 = 6;  x2 = 9.79796;  y = 4

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1) = (32, 18) .// 2
   (y, x1, r2, x2) = ((-4*R^(11/2)*sqrt(r1) + 32*R^(9/2)*r1^(3/2) + 128*R^(7/2)*r1^(5/2) - 2048*R^(5/2)*r1^(7/2) + 7168*R^(3/2)*r1^(9/2) - 8192*sqrt(R)*r1^(11/2) + R^6 - 4*R^5*r1 - 64*R^4*r1^2 + 384*R^3*r1^3 + 256*R^2*r1^4 - 5120*R*r1^5 + 8192*r1^6)/(R^5 - 8*R^4*r1 - 32*R^3*r1^2 + 512*R^2*r1^3 - 1792*R*r1^4 + 2048*r1^5), 2*sqrt(2)*(-R^(11/2)*sqrt(r1) - R^(9/2)*r1^(3/2) + 68*R^(7/2)*r1^(5/2) - 80*R^(5/2)*r1^(7/2) - 1088*R^(3/2)*r1^(9/2) + 2560*sqrt(R)*r1^(11/2) + 4*R^5*r1 - 14*R^4*r1^2 - 200*R^3*r1^3 + 1184*R^2*r1^4 - 1408*R*r1^5 - 1024*r1^6)/(sqrt(R^(3/2)*sqrt(r1) + 8*sqrt(R)*r1^(3/2) - 5*R*r1 - 4*r1^2)*(-R^4 + 4*R^3*r1 + 48*R^2*r1^2 - 320*R*r1^3 + 512*r1^4)), 2*(R^(11/2)*sqrt(r1) - 8*R^(9/2)*r1^(3/2) - 32*R^(7/2)*r1^(5/2) + 512*R^(5/2)*r1^(7/2) - 1792*R^(3/2)*r1^(9/2) + 2048*sqrt(R)*r1^(11/2) - R^5*r1 + 8*R^4*r1^2 + 32*R^3*r1^3 - 512*R^2*r1^4 + 1792*R*r1^5 - 2048*r1^6)/(R^5 - 8*R^4*r1 - 32*R^3*r1^2 + 512*R^2*r1^3 - 1792*R*r1^4 + 2048*r1^5), sqrt(8*R^(3/2)*sqrt(r1) + 64*sqrt(R)*r1^(3/2) - 40*R*r1 - 32*r1^2))
   @printf("小円の直径 = %g;  R = %g;  r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g;  y = %g\n",
      2r2, R, r1, x1, r2, x2, y)
   plot()
   circle(0, 0, R, :magenta)
   circle(x1, y - r1, r1)
   circle(x2, y - r2, r2, :blue)
   circle(0, y + r2, r2, :blue)
   x = sqrt(R^2 - y^2)
   segment(-x, y, x, y)
   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, y - r1, "大円:r1,(x1,y-r1)", :red, :center, :top, delta=-delta/2)
       point(x2, y - r2, "小円:r2,(x2,y-r2)", :blue, :center, :top, delta=-delta/2)
       point(0, y + r2, " 小円:r2\n (x2,y+r2)", :blue, :left, :vcenter)
       point(0, y, " y", :black, :left, :top, delta=-delta/2)
       point(0, R, " R", :magenta, :left, :bottom, delta=delta/2)
   end
end;

算額（その580）

算額（その580）

群馬の算額 19-4 榛名神社 文化8年
http://takasakiwasan.web.fc2.com/gunnsann/g019-4.html

鈎股弦（直角三角形）の中に正方形，大円，小円を入れる。中鈎とは直角の頂点から弦へおろした垂線である。大円，小円の直径が 12 寸，5 寸のとき，正方形の一辺の長さを求めよ。

「鈎」，「股」の長さを，鈎，股とする。
中鈎と弦の交点座標を (x, y) とする。
正方形の一辺の長さを a とする。
大円の半径と中心座標を r1, (股 - a + r1, r1)
小円の半径と中心座標を r2, (股 - r2, a - r2)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, a::positive, r1::positive, x1::positive,
     r2::positive, x::positive, y::positive
(r1, r2) = (12, 5) .// 2
eq1 = a/(股 - a) - 鈎/股
eq2 = distance(股, 0, x, y, 股 - a + r1, r1) - r1^2
eq3 = distance(股, 0, x, y, 股 - r2, a - r2) - r2^2
eq4 = 股*y + 鈎*(股 - x) - 股*鈎
eq5 = a*(鈎 - a) + a*(股 - a) + 2*a^2 - 股*鈎;
# res = solve([eq1, eq2, eq3, eq4, eq5], (鈎, 股, a, x, y))

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, y) = u
   return [
       a/(-a + 股) - 鈎/股,  # eq1
       (-y*(-a*x + a*股 + 6*x + 6*y - 6*股)/(x^2 - 2*x*股 + y^2 + 股^2) + 6)^2 + (-a + 股 - (股*(x^2 - 2*x*股 + y^2 + 股^2) + (x - 股)*(-a*x + a*股 + 6*x + 6*y - 6*股))/(x^2 - 2*x*股 + y^2 + 股^2) + 6)^2 - 36,  # eq2
       (a - y*(2*a*y - 5*x - 5*y + 5*股)/(2*(x^2 - 2*x*股 + y^2 + 股^2)) - 5/2)^2 + (股 - (股*(x^2 - 2*x*股 + y^2 + 股^2) + (x - 股)*(2*a*y - 5*x - 5*y + 5*股)/2)/(x^2 - 2*x*股 + y^2 + 股^2) - 5/2)^2 - 25/4,  # eq3
       y*股 - 股*鈎 + 鈎*(-x + 股),  # eq4
       2*a^2 + a*(-a + 股) + a*(-a + 鈎) - 股*鈎,  # eq5
   ]
end;

(r1, r2) = (12, 5) .// 2
iniv = BigFloat[20.6, 50.1, 14.5, 42.1, 17.9]
res = nls(H, ini=iniv)

   (BigFloat[21.35733319823884343090635521344017153294802876982109796830915235365518705818788, 50.39219873866784319300311036791102016301992183128791550477294030211998940687325, 15.00000000000000000000000000000000000000000000000000000000000000000000000000028, 42.82892754606824024613300750612840676283894513302086625119773230436335822945181, 18.15185086223904707248824686827827216043434407584091820858049919461591482581305], true)

正方形の一辺の長さは 15 寸である。

他のパラメータは以下の通り。
鈎 = 21.3573;  股 = 50.3922;  a = 15;  x1 = 41.3922;  x = 42.8289;  y = 18.1519

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (12, 5) .// 2
   (鈎, 股, a, x1, x, y) = [98, 238, 69, 195, 200, 85] .* (12/57)
   (鈎, 股, a, x, y) = res[1]
   @printf("鈎 = %g;  股 = %g;  a = %g;  x1 = %g;  x = %g;  y = %g\n",
       鈎, 股, a, x1, x, y)
   plot([0, 股, 股, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(股 - a + r1, r1, r1)
   circle(股 - r2, a - r2, r2, :blue)
   segment(股, 0, x, y)
   rect(股 - a, 0, 股, a, :green)
   plot!(xlims=(-1, 57), ylims=(-2.5, 25))
   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(股 - a + r1, r1, "大円:r1\n(股-a+r1,r1)", :red, :center, :bottom, delta=delta)
       point(股 - r2, a - r2, "小円:r2,(股-r2,a-r2) ", :blue, :right, :vcenter)
       point(股, a, " (股,a)", :green, :left, :vcenter)
       point(x, y, "(x,y)", :black, :right, :bottom, delta=delta)
       point(股 - a, 0, "股-a", :green, :center, :top, delta=-delta/2)
       point(股, 0, "股", :green, :center, :top, delta=-delta/2)
       point(股, 鈎, "(股,鈎) ", :black, :right, :bottom, delta=delta/2)
   end
end;

算額（その579）

算額（その579）

群馬の算額 19-6 榛名神社 文化8年
http://takasakiwasan.web.fc2.com/gunnsann/g019-6.html

3 本の線分で 3 個の円を挟んでいる。大円，小円の直径が 4 寸，1 寸 4 分 4 厘 のとき，界斜（図の記号では点 (x,0) と (x0, y0) を結ぶ線分）の長さを求めよ。

界斜の両端の座標を (x,0) と (x0, y0) とする。
大円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r1, (x11, r1), (x12, r1)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, x11::positive, x12::positive, r2::positive, x2::positive,
     x::positive, x0::positive, y0::positive
(r1, r2) = (144//200, 4//2)
eq1 = (x2 - x12)^2 + (r2 - r1)^2 - (r1 + r2)^2
eq2 = distance(0, 0, x0, y0, x11, r1) - r1^2
eq3 = distance(0, 0, x0, y0, x2, r2) - r2^2    
eq4 = distance(x, 0, x0, y0, x11, r1) - r1^2
eq5 = distance(x, 0, x0, y0, x12, r1) - r1^2
eq6 = distance(x, 0, x0, y0, x2, r2) - r2^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (x11, x12, x2, x, x0, y0))

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)
   (x11, x12, x2, x, x0, y0) = u
   return [
       (-x12 + x2)^2 - 144/25,  # eq1
       (-x0*(25*x0*x11 + 18*y0)/(25*(x0^2 + y0^2)) + x11)^2 + (-y0*(25*x0*x11 + 18*y0)/(25*(x0^2 + y0^2)) + 18/25)^2 - 324/625,  # eq2
       (-x0*(x0*x2 + 2*y0)/(x0^2 + y0^2) + x2)^2 + (-y0*(x0*x2 + 2*y0)/(x0^2 + y0^2) + 2)^2 - 4,  # eq3
       (x11 - (x11*(x^2 - 2*x*x0 + x0^2 + y0^2) + y0*(25*x*y0 - 18*x + 18*x0 - 25*x11*y0)/25)/(x^2 - 2*x*x0 + x0^2 + y0^2))^2 + (-y0*(25*x^2 - 25*x*x0 - 25*x*x11 + 25*x0*x11 + 18*y0)/(25*(x^2 - 2*x*x0 + x0^2 + y0^2)) + 18/25)^2 - 324/625,  # eq4
       (x12 - (x^2*x12 - 2*x*x0*x12 + x*y0^2 - 18*x*y0/25 + x0^2*x12 + 18*x0*y0/25)/(x^2 - 2*x*x0 + x0^2 + y0^2))^2 + (-y0*(25*x^2 - 25*x*x0 - 25*x*x12 + 25*x0*x12 + 18*y0)/(25*(x^2 - 2*x*x0 + x0^2 + y0^2)) + 18/25)^2 - 324/625,  # eq5
       (x2 - (x^2*x2 - 2*x*x0*x2 + x*y0^2 - 2*x*y0 + x0^2*x2 + 2*x0*y0)/(x^2 - 2*x*x0 + x0^2 + y0^2))^2 + (-y0*(x^2 - x*x0 - x*x2 + x0*x2 + 2*y0)/(x^2 - 2*x*x0 + x0^2 + y0^2) + 2)^2 - 4,  # eq6
   ]
end;

(r1, r2) = (144//200, 4//2)
iniv = BigFloat[2.2, 4, 6.4, 2.6, 5.412, 3.8]
res = nls(H, ini=iniv)

   (BigFloat[2.32537499999999981452205370806705695806346869674534523102666679819954313421592, 4.05937499999999973306648353848901437850698923197181920887229922564784525145764, 6.459374999999999688657562553482752350696667705728683787732940545284447986400679, 2.709374999999999822647282592843187726735573892402228819331695263097710712289562, 5.012399221453286844773215617274192837885357156437953123416192597373948544218003, 3.433079584775086464352099968927128171012730897483676505879678311456937136579126], true)

界斜の長さは sqrt((x0 - x)^2 + y0^2) ≒ 4.134である。

x  = 2.7093749999999998
x0 = 5.0123992214532868
y0 = 3.4330795847750864
sqrt((x0 - x)^2 + y0^2)

   4.134

他のパラメータは以下の通り。
r1 = 0.72;  r2 = 2;  x11 = 2.32537;  x12 = 4.05937;  x2 = 6.45937

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (144//200, 4//2)
   (x11, x12, x2, x, x0, y0) = res[1]
   l = sqrt((x0 - x)^2 + y0^2)
   @printf("界斜の長さ = %g;  r1 = %g;  r2 = %g;  x11 = %g;  x12 = %g;  x2 = %g;  x = %g;  x0 = %g;  y0 = %g\n",
       l, r1, r2, x11, x12, x2, x, x0, y0)
   plot()
   circle(x11, r1, r1)
   circle(x12, r1, r1)
   circle(x2, r2, r2, :blue)
   abline(0, 0, y0/x0, 0, 6.5)
   segment(x, 0, x0, y0)
   plot!(ylims=(-0.3, 4.5))
   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(x11, r1, "小円:r1\n(x11,r2)", :red, :center, :bottom, delta=delta)
       point(x12, r1, "小円:r1\n(x12,r2)", :red, :center, :bottom, delta=delta)
       point(x2, r2, "大円:r2,(x2,r2)", :magenta, :center, :top, delta=-delta)
       point(x, 0, "x", :black, :center, :top, delta=-delta/2)
       point(x0, y0, "(x0,y0)", :black, :right, :bottom, delta=delta/2)
   end
end;

算額（その578）

算額（その578）

群馬の算額 17-5 八幡宮
http://takasakiwasan.web.fc2.com/gunnsann/g017-5.html

直線，垂線と斜線に大円，中円，小円が接している。大円，中円の直径が 5 寸，3 寸のとき，小円の直径を求めよ。

斜線の x 切片，y 切片を x, y とする。
大円の半径と中心座標を r1, (r1, r1)
中円の半径と中心座標を r2, (-r2, r2)
小円の半径と中心座標を r3, (r3, y3)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive, x::positive, y::positive

eq1 = distance(x, 0, 0, y, r1, r1) - r1^2
eq2 = distance(x, 0, 0, y, -r2, r2) - r2^2
eq3 = distance(x, 0, 0, y, r3, r3) - r3^2
res = solve([eq1, eq2, eq3], (r3, x, y))

   4-element Vector{Tuple{Sym, Sym, Sym}}:
    (r1, -2*r1*r2/(r1 - r2), r1 + r2)
    (r1, r1 - r2, 2*r1*r2/(r1 + r2))
    (-r2*(r1 + r2)/(r1 - r2), -2*r1*r2/(r1 - r2), r1 + r2)
    (r2*(r1 - r2)/(r1 + r2), r1 - r2, 2*r1*r2/(r1 + r2))

4 番目の解が適解である。

小円の直径は r2*(r1 - r2)/(r1 + r2) で，r1 = 5//2, r2 = 3//2 のときは，3/4 = 0.75 = 7 分 5 厘である。

その他のパラメータは以下の通り。
小円の直径 = 0.75;  r1 = 2.5;  r2 = 1.5;  r3 = 0.375;  x = 1;  y = 1.875

(r1, r2) = (5, 3) .// 2
r2*(r1 - r2)/(r1 + r2)

   3//8

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (5, 3) .// 2
   (r3, x, y) = (r2*(r1 - r2)/(r1 + r2), r1 - r2, 2*r1*r2/(r1 + r2))
   @printf("小円の直径 = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x = %g;  y = %g\n", 2r3, r1, r2, r3, x, y)
   plot()
   circle(r1, r1, r1)
   circle(-r2, r2, r2, :blue)
   circle(r3, r3, r3, :magenta)
   abline(0, y, -y/x, -5, 5)
   plot!(xlims=(-3.5, 5.2), ylims=(-0.5, 5))
   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(r1, r1, "大円:r1,(r1,r1)", :red, :center, :top, delta=-delta)
       point(-r2, r2, "中円:r2,(-r2,r2)", :blue, :center, :top, delta=-delta)
       point(r3, r3, "  小円:r3,(r3,r3)", :magenta, :left, :bottom, delta=5delta)
       point(x, 0, " x", :black, :left, :bottom, delta=delta/2)
       point(0, y, "  y", :black, :left, :vcenter)
   end
end;

算額（その577）

算額（その577）

群馬の算額 18-5 八幡宮
http://takasakiwasan.web.fc2.com/gunnsann/g018-5.html

2 本の直線と 1 本の垂線に接している甲円，乙円，丙円がある。乙円，丙円の直径が 12 寸，3 寸のとき，甲円の直径を求めよ。

甲円の半径と中心座標を r1, (-r1, r1)
乙円の半径と中心座標を r2, (r2, r2)
丙円の半径と中心座標を r3, (r3, y3)
とする。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive, y3::positive,
     a::negative, b::positive, x::positive, y::positive

(r2, r3) = (12, 3) .// 2

   (6//1, 3//2)

乙円と丙円が外接することから，丙円の中心座標 y3 が求まる。

eq1 = (r2 - r3)^2 + (y3 - r2)^2 - (r2 + r3)^2
y3 = solve(eq1)[1]
y3 |> println

   12

乙円と丙円の中心を通る直線の傾き a と y 切片 b を求める。

eq2 = r3*a + b - y3
eq3 = r2*a + b - r2
res = solve([eq2, eq3], (a, b))

   Dict{Any, Any} with 2 entries:
     b => 14
     a => -4/3

甲円と丙円の共通接線と乙円と丙円の中心を通る直線の y 切片は同じ 14 である。
甲円と乙円の中心座標と共通接線の距離についての式を立てる。共通接線は (0, 14) と (x, y) を通る。x は任意で(
例えば 10 とか）， y は x に応じて決まる。

eq5 = distance(0, 14, 10, y, -r1, r1) - r1^2
eq6 = distance(0, 14, 10, y, r2, r2) - r2^2
res = solve([eq5, eq6])

   1-element Vector{Dict{Any, Any}}:
    Dict(r1 => 8, y => 133/12)

甲円の半径は 8，直径は 16 寸である。

共通接線は (0, 14) と (10, 133/12) を通る。傾きは -7//24

(133//12 - 14)//10 

   -7//24

他のパラメータは以下の通り。

   甲円の直径 = 16;  r1 = 8;  r2 = 6;  r3 = 1.5;  y3 = 12;  y 切片 = 14;  傾き = -0.291667

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3) = (12, 3) .// 2
   y3 = 12
   r1 = 8
   intercept = 14
   slope = -7//24
   @printf("甲円の直径 = %g;  r1 = %g;  r2 = %g;  r3 = %g;  y3 = %g;  y 切片 = %g;  傾き = %g\n", 2r1, r1, r2, r3, y3, intercept, slope)
   plot()
   circle(r2, r2, r2)
   circle(r3, y3, r3, :blue)
   circle(-r1, r1, r1, :magenta)
   abline(0, intercept, slope, -18, 15)
   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(-r1, r1, "甲円:r1,(-r1,r1)", :magenta, :center, :top, delta=-delta)
       point(r2, r2, "乙円:r2,(r2,r2)", :red, :center, :top, delta=-delta)
       point(r3, y3, "丙円:r3,(r3,y3)", :blue, :left, :bottom, delta=7delta)
       point(0, 14, "14 ", :black, :right, :top)
       point(10, 133/12, "(10,133/12)", :black, :left, :bottom, delta=2delta)
   end
end;

 

算額（その576）

算額（その576）

群馬の算額 25-1 木部村鎮守社
http://takasakiwasan.web.fc2.com/gunnsann/g025-1.html

大円と小円が交わっている。大円と小円の直径の和が 21 寸 6 分，大小の矢の和が 3 寸 6 分，共通弦の長さが 8 寸のとき，大円の直径はいかほどか。

大円の半径と中心座標を r1, (0, 0)
小円の半径と中心座標を r2, (0, y2)
大円と小円の交点座標を (x, y)
とおき，以下の連立方程式を解く。
大矢と小矢は (r1 - y), (y - y2 + r2)，共通弦の長さは 2x である。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, y2::positive, x::positive, y::positive

x = 8//2
eq1 = r1 + r2 - 216/20
eq2 = (r1 - y) + (y - (y2 - r2)) - 36/10
eq3 = (r1^2 - y^2) - x^2
eq4 = (r2^2 - (y2 - y)^2) - x^2
res = solve([eq1, eq2, eq3, eq4], (r1, r2, y2, y))

   2-element Vector{NTuple{4, Sym}}:
    (5.00000000000000, 5.80000000000000, 7.20000000000000, 3.00000000000000)
    (5.80000000000000, 5.00000000000000, 7.20000000000000, 4.20000000000000)

r1 > r2 なので，2 番目の解が適解である。

   r1 = 5.8;  r2 = 5;  y2 = 7.2;  x = 4;  y = 4.2

大円の直径は 11 寸 6 分である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   x = 8//2
   (r1, r2, y2, y) = res[2]
   @printf("大円の直径 = %g;  r1 = %g;  r2 = %g;  y2 = %g;  x = %g;  y = %g\n", 2r1, r1, r2, y2, x, y)
   plot()
   circle(0, 0, r1)
   circle(0, y2, r2, :blue)
   segment(-x, y, x, y, :magenta)
   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(r1, 0, "大円:r1,(r1,0)  r1 ", :red, :right, :bottom, delta=delta/2)
       point(0, y2, " y2\n 小円:r2,(0,y2)", :blue, :left, :vcenter)
       point(x, y, " (x,y)", :black, :left, :vcenter)
       point(0, y2 - r2, " y2-r2", :blue, :left, :top, delta=-delta/2)
       point(0, r1, " r1", :red, :left, :top, delta=-delta/2)
   end
end;

算額（その575）

算額（その575）

福岡県久留米高良山者
藤田貞資(1789): 神壁算法
http://hyonemitsu.web.fc2.com/Korataisha.pdf

外円内に斜線を隔てて甲円 1 個，乙円 2 個，丙円 1 個がある。外円，甲円，丙円の直径が 60 寸，20 寸，28 寸のとき，乙円の直径はいかほどか。

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

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive, x::positive, y::positive,
     r3::positive, x1::positive, x2::positive

(R, r1, r3) = (60, 20, 28) .// 2
eq1 = distance(-x1, sqrt(R^2 - x1^2), x2, -sqrt(R^2 - x2^2), 0, R - r1)  - r1^2
eq2 = distance(-x1, sqrt(R^2 - x1^2), x2, -sqrt(R^2 - x2^2), x, y) - r2^2
eq3 = distance(-x1, sqrt(R^2 - x1^2), x2, -sqrt(R^2 - x2^2), 0, r3 - R) - r3^2
eq4 = distance(-x2, -sqrt(R^2 - x2^2), x1, sqrt(R^2 - x1^2), x, y) - r2^2
eq5 = x^2 + y^2 - (R - r2)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r2, x2, y2, x1, x2))

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)
   (r2, x, y, x1, x2) = u
   return [
       (20 - (-x1^2*sqrt(900 - x2^2) - 20*x1^2 + x1*x2*sqrt(900 - x1^2) - x1*x2*sqrt(900 - x2^2) + x2^2*sqrt(900 - x1^2) - 20*x2^2 + 40*sqrt(900 - x1^2)*sqrt(900 - x2^2) + 36000)/(2*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)))^2 + (-x1*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) + (x1 + x2)*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) - 20*sqrt(900 - x1^2) - 20*sqrt(900 - x2^2) + 900)/2)^2/(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)^2 - 100,  # eq1
       -r2^2 + (x - (-x1*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) + (x1 + x2)*(x*x1 + x*x2 + x1*x2 - y*sqrt(900 - x1^2) - y*sqrt(900 - x2^2) + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)/2)/(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900))^2 + (y - (2*sqrt(900 - x1^2)*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) - (sqrt(900 - x1^2) + sqrt(900 - x2^2))*(x*x1 + x*x2 + x1*x2 - y*sqrt(900 - x1^2) - y*sqrt(900 - x2^2) + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900))/(2*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)))^2,  # eq2
       (-16 - (-x1^2*sqrt(900 - x2^2) + 16*x1^2 + x1*x2*sqrt(900 - x1^2) - x1*x2*sqrt(900 - x2^2) + x2^2*sqrt(900 - x1^2) + 16*x2^2 - 32*sqrt(900 - x1^2)*sqrt(900 - x2^2) - 28800)/(2*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)))^2 + (-x1*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) + (x1 + x2)*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 16*sqrt(900 - x1^2) + 16*sqrt(900 - x2^2) + 900)/2)^2/(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)^2 - 196,  # eq3
       -r2^2 + (x - (-x2*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) + (x1 + x2)*(x*x1 + x*x2 + x1*x2 + y*sqrt(900 - x1^2) + y*sqrt(900 - x2^2) + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)/2)/(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900))^2 + (y - (-2*sqrt(900 - x2^2)*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900) + (sqrt(900 - x1^2) + sqrt(900 - x2^2))*(x*x1 + x*x2 + x1*x2 + y*sqrt(900 - x1^2) + y*sqrt(900 - x2^2) + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900))/(2*(x1*x2 + sqrt(900 - x1^2)*sqrt(900 - x2^2) + 900)))^2,  # eq4
       x^2 + y^2 - (30 - r2)^2,  # eq5
   ]
end;

iniv = BigFloat[29, 40, 12, 21, 28] .* (30/70)
res = nls(H, ini=iniv)

   (BigFloat[12.50000000000000000000000000000000000000000000000000000001913542531009364124636, 16.77050983124842272306880251548457176580463769708644293204781524730502299339952, 4.999999999999999999999999999999999999999999999999999999963738320809945557470043, 17.39163982499836430540468409013214849787147613031186674435362589343567963493901, 22.3606797749978969640917366873127623544061835961152572427143001079596068757832], true)

乙円の直径は 25 寸である。

その他のパラメータは以下の通り。
乙円の直径 = 25;  r2 = 12.5;  x = 16.7705;  y = 5;  x1 = 17.3916;  x2 = 22.3607

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1, r3) = (60, 20, 28) .// 2
   (r2, x, y, x1, x2) = res[1]
   @printf("乙円の直径 = %g;  r2 = %g;  x = %g;  y = %g;  x1 = %g;  x2 = %g\n", 2r2, r2, x, y, x1, x2)
   plot()
   circle(0, 0, R)
   circle(0, R - r1, r1, :blue)
   circle(x, y, r2, :magenta)
   circle(-x, y, r2, :magenta)
   circle(0, r3 - R, r3, :green)
   segment(-x1, sqrt(R^2 - x1^2), x2, -sqrt(R^2 - x2^2), :gray)
   segment(x1, sqrt(R^2 - x1^2), -x2, -sqrt(R^2 - x2^2), :gray)
   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 - r1, " 甲円:r1\n (0,R-r1)", :blue, :left, :vcenter)
       point(0, r3 - R, " 丙円:r1\n (0,r3-R)", :green, :left, :vcenter)
       point(x, y, "乙円:r2,(x,y)", :magenta, :center, :top, delta=-delta/2)
       point(x1, √(R^2 - x1^2), "(x1,√(R^2 - x1^2))", :black, :center, :bottom, delta=delta/2)
       point(x2, -√(R^2 - x2^2), "(x2,-√(R^2 - x2^2))", :black, :center, :bottom, delta=delta/2)
   end
end;

算額（その574）

算額（その574）

佐賀県神埼市 冠者神社（米光丁奉納） 平成21年
http://www.wasan.jp/saga/kanjya.html

大小の直角三角形の中に黃円，赤円が入っている。黃円の直径が 2 寸，甲斜が 4 寸，丁斜が 12 寸のとき，赤円の直径を求めよ。

乙斜，丙斜，戊斜の長さ を 「乙斜」，「丙斜」，「戊斜」
赤円の半径を r2
として，以下の eq1〜eq4 の連立方程式を解けば「乙斜，丙斜，戊斜, 赤円の半径」が求まる。この場合には座標軸の設定は不要である。
なお，図を描くために以下のように座標軸を設定し，追加パラメータとして，小さい方の直角三角形の頂点を (x, y)，黃円の中心座標を （x1，y1); y1 = 13 として eq1〜eq7 をまとめて解く。

算額の常套手段では一般解の式を求めて最後に具体的な数値を代入して解を求めるのであるが，3 変数を式に含めたままの一般式を求めるのはあまり現実的でないので，最初から具体的に数値を使用して連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms 甲斜::positive, 丁斜::positive, 乙斜::positive, 丙斜::positive, 戊斜::positive, r2::positive,
     x::positive, y::positive, r1::positive, x1::positive

(r1, 甲斜, 丁斜) = (1, 4, 12)
eq1 = 乙斜 + 甲斜 - 丙斜 - 2r1
eq2 = 丙斜 + 丁斜 - 戊斜 - 2r2
eq3 = 乙斜^2 + 甲斜^2 - 丙斜^2
eq4 = 丙斜^2 + 丁斜^2 - 戊斜^2
eq5 = x^2 + (y - 丁斜)^2 - 乙斜^2
eq6 = (丙斜 - x)^2 + (y - 丁斜)^2 - 甲斜^2
eq7 = distance(0, 丁斜, x, y, x1, 丁斜 + r1) - r1^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (乙斜, 丙斜, 戊斜, r2, x, y, x1))

   2-element Vector{NTuple{7, Sym}}:
    (3, 5, 13, 2, 9/5, 48/5, 1/2)
    (3, 5, 13, 2, 9/5, 72/5, 2)

2 通りの解が得られるが，2 番目のものが適解である。
黃円の直径が 2 寸，甲斜が 4 寸，丁斜が 12 寸のとき，赤円の直径は 4 寸である。

その他のパラメータは以下の通り。
乙斜 = 3;  丙斜 = 5;  戊斜 = 13;  r2 = 2;  x = 1.8;  y = 14.4;  x1 = 2

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, 甲斜, 丁斜) = (1, 4, 12)
   (a, b, c, r2, x, y, x1) = (3, 5, 13, 2, 9/5, 72/5, 2)
   @printf("赤円の直径 = %g;  乙斜 = %g;  丙斜 = %g;  戊斜 = %g;  r2 = %g;  x = %g;  y = %g;  x1 = %g\n", 2r2, a, b, c, r2, x, y, x1)
   plot([b, x, 0, b, 0, 0], [12, y, 12, 12, 0, 12], color=:blue, lw=0.5)
   circle(r2, 12- r2, r2)
   circle(x1, 13, 1, :green)
   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(x, y, "(x,y)", :blue, :center, :bottom, delta=delta)
       point(b, 12, " (b,12)", :blue, :left, :vcenter)
       point(x1, 12 + r1, "黃円:r1,(x1,12+r1)", :green, :center, delta=-delta/2)
       point(r2, 12 - r2, "赤円:r2\n(r2,12-r2)", :red, :center, delta=-delta/2)
       plot!(xlims=(-0.5, 5.5))
       point(3, 14, "甲斜", :black, mark=false)
       point(0.1, 14, "乙斜", :black, mark=false)
       point(3.5, 11.8, "丙斜", :black, mark=false)
       point(0.2, 7, "丁斜", :black, mark=false)
       point(3, 7, "戊斜", :black, mark=false)
   end
end;

算額（その573）

算額（その573）

七十二 群馬県富岡市一ノ宮 貫前神社 嘉永2年(1849)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

群馬の算額 72-6 貫前神社 嘉永2年
http://takasakiwasan.web.fc2.com/gunnsann/g072-6.html

楕円の中に 2 個の合同な菱形と，2 個の合同な円を入れる。円は楕円の長径の端で楕円に接する最大のものとする。楕円の短径が 1 寸のとき，菱形の一辺の長さを求めよ。

楕円の長径と短径を a, b とする。
円は曲率円なので，その半径は b^2/a である。
円の半径と中心座標を r, (0, a - r); r = b^2/a
第一象限にある菱形の頂点座標を (x, y); y = b/2
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, x::positive, y::positive,
     r::positive, c::positive
r = b^2/a
y = b/2
eq1 = x^2/a^2 + y^2/b^2 - 1
eq2 = y/sqrt(x^2 + y^2) - r/(a - r);
res = solve([eq1, eq2], (a, x))

   1-element Vector{Tuple{Sym, Sym}}:
    (sqrt(5)*b, sqrt(15)*b/2)

長径は短径の √5 倍である。
短径が 1/2 寸のとき，長径は a = (1/2)√5 ≒ 1.11803，r = (1/2)^2 / ((1/2)√5) ≒ 0.223607，x = (1/2)*√15/2 ≒ 0.968246，y = 1/4 である。

a = (1/2)^2/√5
r = (1/2)^2 / ((1/2)√5)
x = (1/2)*√15/2
y = 1/4
(a, r, x, y)

   (0.11180339887498948, 0.22360679774997896, 0.9682458365518543, 0.25)

菱形の一辺の長さは sqrt(x^2 + y^2) = 1 である。

sqrt(x^2 + y^2)

   1.0

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = 1//2
   (a, x) = (sqrt(5)*b, sqrt(15)*b/2)
   y = b/2
   r = b^2/a
   c = sqrt(x^2 + y^2)
   @printf("a = %g;  r = %g;  x = %g;  y = %g;  c = %g\n", a, r, x, y, c)
   plot()
   ellipse(0, 0, a, b, color=:red)
   circle(a - r, 0, r, :green)
   circle(r - a, 0, r, :green)
   plot!([0, x, 0, -x, x, 0, -x, 0], [0, y, b, y, -y, -b, -y, 0], color=:blue, lw=0.5)
   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(x, y, " (x,y)", :green, :left, :vcenter)
       point(a - r, 0, "a-r", :green, :center, :bottom, delta=delta/2)
       point(a, 0, "a ", :green, :right, :bottom, delta=delta/2)
       point(0, b, " b", :green, :left, :bottom, delta=delta/2)
   end
end;

 

算額（その572）

算額（その572）

群馬の算額 19-2 榛名神社 文化8年
http://takasakiwasan.web.fc2.com/gunnsann/g019-2.html

正方形の中に 3 本の斜線を引き，区画された領域に大円，中円，小円を入れる。正方形の一辺の長さが 3 寸のとき，大円の直径を求めよ。

正方形の一辺の長さを a とする。
大円の半径と中心座標を r1, (r1, a - r1)
中円の半径と中心座標を r2, (a - r2, r2)
小円の半径と中心座標を r3, (r3, r3), (a - r3, a - r3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive, r3::positive,
     b::positive, c::positive, d::positive, e::positive
eq1 = distance(d, 0, 0, e, r1, a - r1) - r1^2 |> expand |> simplify |> numerator
eq2 = distance(d, 0, 0, e, a - r2, r2) - r2^2 |> expand |> simplify |> numerator
eq3 = distance(d, 0, 0, e, r3, r3) - r3^2 |> expand |> simplify |> numerator
eq4 = distance(b, 0, a, c, r1, a - r1) -r1^2 |> expand |> simplify |> numerator
eq5 = distance(b, 0, a, c, a - r2, r2) -r2^2 |> expand |> simplify |> numerator
eq6 = distance(b, 0, a, c, r3, r3) -r3^2 |> expand |> simplify |> numerator
eq7 = distance(b, 0, a, c, a - r3, a - r3) -r3^2 |> expand |> simplify |> numerator;

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)
   (r1, r2, r3, b, c, d, e) = u
   return [
       a^2*d - 2*a*d*e - 2*a*d*r1 + 2*a*e*r1 + d*e^2 + 2*d*e*r1 - 2*e^2*r1 - 2*e*r1^2,  # eq1
       a^2*e - 2*a*d*e + 2*a*d*r2 - 2*a*e*r2 + d^2*e - 2*d^2*r2 + 2*d*e*r2 - 2*d*r2^2,  # eq2
       d*e - 2*d*r3 - 2*e*r3 + 2*r3^2,  # eq3
       a^4 - 2*a^3*b - 2*a^3*r1 + a^2*b^2 + 2*a^2*b*c + 4*a^2*b*r1 - 2*a^2*c*r1 - 2*a*b^2*c - 2*a*b^2*r1 + 2*a*c*r1^2 + b^2*c^2 + 2*b^2*c*r1 - 2*b*c^2*r1 - 2*b*c*r1^2,  # eq4
       a^2*c - 2*a^2*r2 - 2*a*b*c + 4*a*b*r2 - 2*a*c*r2 + 2*a*r2^2 + b^2*c - 2*b^2*r2 + 2*b*c*r2 - 2*b*r2^2,  # eq5
       2*a*b*r3 - 2*a*r3^2 + b^2*c - 2*b^2*r3 - 2*b*c*r3 + 2*b*r3^2,  # eq6
       a^4 - 2*a^3*b - 2*a^3*c - 2*a^3*r3 + a^2*b^2 + 4*a^2*b*c + 4*a^2*b*r3 + a^2*c^2 + 4*a^2*c*r3 - 2*a*b^2*c - 2*a*b^2*r3 - 2*a*b*c^2 - 6*a*b*c*r3 - 2*a*c^2*r3 - 2*a*c*r3^2 + b^2*c^2 + 2*b^2*c*r3 + 2*b*c^2*r3 + 2*b*c*r3^2,  # eq7
   ]
end;
a = 3
iniv = BigFloat[1.05, 0.7, 0.48, 0.7, 2.5, 2.0, 1.4]
res = nls(H, ini=iniv)
   (BigFloat[1.074953912118040074148850906669271748631882063090927359480694112255763934484609, 0.6824054007626747794460829207016340156591023107782284209892666737724326939646588, 0.4738479700017459009683547528030295288720148655840843288085599668979468150150135, 0.670122225679428548319738745786848395358404741740433798054273686137994090593031, 2.329877774320571451680261254213151604641595258259566201945726313849903192359496, 2.121320343559642573202533086314547117854507813065422109765019606987497365607011, 1.330325847824819753801482812391145847210905426272977123802622116991951746211069], true)

大円の直径は 2.14991 寸である。

その他のパラメータは以下の通り。

r1 = 1.07495;  r2 = 0.682405;  r3 = 0.473848;  b = 0.670122;  c = 2.32988;  d = 2.12132;  e = 1.33033

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, b, c, d, e) = [1.05, 0.7, 0.48, 0.7, 2.5, 2.0, 1.4]
   (r1, r2, r3, b, c, d, e) = res[1]
   @printf("大円の直径 = %g;  r1 = %g;  r2 = %g;  r3 = %g;  b = %g;  c = %g;  d = %g;  e = %g\n",
       2r1, r1, r2, r3, b, c, d, e)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:blue, lw=0.5)
   circle(r1, a - r1, r1)
   circle(a - r2, r2, r2, :green)
   circle(r3, r3, r3, :magenta)
   circle(a - r3, a - r3, r3, :magenta)
   segment(b, 0, a, c)
   segment(d, 0, 0, e)
   segment(a, a - d, a - e, a)
   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(b, 0, "b ", :black, :right, delta=3delta)
       point(d, 0, "d", :black, :left, delta=3delta)
       point(a, 0, "a ", :black, :right, delta=3delta)
       point(0, e, " e", :black, :left, delta=2delta)
       point(a, c, "(a,c) ", :black, :right, :vcenter)
       point(a, a - d, "(a,a-d) ", :black, :right, :vcenter)
       point(a - e, a, "(a-e,a)", :black, :center, :bottom, delta=delta/2)
       point(0, a, " a", :black, :left, :bottom, delta=delta/2)
       point(r1, a - r1, "大円:r1,(r1,a-r1)", :red, :center, delta=-delta/2)
       point(a - r2, r2, "中円:r2,(a-r2,r2)", :green, :center, delta=-delta/2)
       point(r3, r3, "小円:r3,(r3,r3)", :magenta, :center, delta=-delta/2)
       point(a - r3, a - r3, "小円:r3,(a-r3,a-r3)", :magenta, :center, delta=-delta/2)
   end
end;

算額（その571）

算額（その571）

群馬の算額 19-7 榛名神社 文化8年
http://takasakiwasan.web.fc2.com/gunnsann/g019-7.html

正三角形の中に円弧を描き，大円 1 個と黒円 2 個を入れる。正三角形の長さが 10 寸のとき，最小になる黒円の直径を求めよ。

黒円の直径が最小になるのは，弧が正三角形の頂点で接するときである。
正三角形の一辺の長さを 2a とする。
弧の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, R + r1)
黒円の半径と中心座標を r2, (x2, y2)
とおき，以下の連立方程式を解く。

まず R を求める。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, x2::positive, y2::positive,
     R::positive, y::positive, a::positive
eq0 = R^2 + 4a^2 - (sqrt(R^2 - a^2) + sqrt(Sym(3))a)^2
R = solve(eq0, R)  # 10/sqrt(Sym(3))
R |> println

   Sym[2*sqrt(3)*a/3]

次いで，R を既知として残りのパラメータを求める。

@syms r1::positive, r2::positive, x2::positive, y2::positive,
     R::positive, y::positive, a::positive
R = 2a/sqrt(Sym(3))
eq1 = x2^2 + (R + r1 - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + y2^2 - (R + r2)^2
eq3 = r1/(y + sqrt(Sym(3))a - (R + r1)) - 1//2
eq4 = a^2 + y^2 - R^2
eq5 = 2a*r1 + (R + r1 - y)a - sqrt(Sym(3))a^2
eq6 = (sqrt(Sym(3))a*x2 + a*(y2 - y) + 2a*r2) - sqrt(Sym(3))a^2;
eq7 = distance(0, y + sqrt(Sym(3))a, a, y, x2, y2) - r2^2
eq8 = distance(0, y + sqrt(Sym(3))a, a, y, 0, R + r1) - r1^2

res = solve([eq1, eq2, eq3, eq4, eq6], (r1, r2, x2, y2, y))

   1-element Vector{NTuple{5, Sym}}:
    (2*sqrt(3)*a/9, a*(-3 + 2*sqrt(3))/3, a*(-6 + 5*sqrt(3))/6, a*(-1/2 + sqrt(3)), sqrt(3)*a/3)

黒円の半径は a * (2√3 - 3)/3 である。黒円の直径は正三角形の一辺の長さの (2√3 - 3)/3 倍である。
正三角形の一辺の長さが 10 寸のとき，1.5470053837925146 寸である。

10 * (2√3 - 3)/3

   1.5470053837925146

その他のパラメータは以下の通り。

   黒円の直径 = 1.54701;  R = 5.7735;  r1 = 1.9245;  r2 = 0.773503;  x2 = 2.21688;  y2 = 6.16025;  R = 5.7735;  y = 2.88675

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 10//2
   R = 2a/sqrt(Sym(3))
   (r1, r2, x2, y2, y) = a .* (2√3/9, (2√3 - 3)/3, (5√3 - 6)/6, √3 - 1/2, √3/3)
   @printf("黒円の直径 = %g;  R = %g;  r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g;  R = %g;  y = %g\n", 2r2, R, r1, r2, x2, y2, R, y)
   plot([a, 0, -a, a], y .+ [0, √3a, 0, 0], color=:blue, lw=0.5)
   circle(0, 0, R, beginangle=0, endangle=180)
   circle(0, R + r1, r1, :magenta)
   circle(x2, y2, r2, :green)
   circle(-x2, y2, r2, :green)
   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, " R", :red, :left, :top, delta=-delta)
       point(0, y, " y", :blue, :left, :top, delta=-delta)
       point(0, R + r1, " 大円:r1\n (0,R+r1)")
       point(0, y + √3a, " y+√3a", :blue, :left, :vcenter)
       point(x2, y2, " 黒円:r2,(x2,y2)", :black, :left, :vcenter)
       point(a, y, "(a,y)", :blue, :right, :top, delta=-delta/2)
       point(0, 0, " 弧:R,(0,0)", :red, :left, :bottom, delta=delta/2)
   end
end;