算額（その1328）

算額（その1328）

七八 加須市大字外野 棘脱地蔵堂 明治9年(1876)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円7個，正方形
#Julia, #SymPy, #算額, #和算

正方形の中に大円，中円，小円を容れる。大円の直径が 10 寸，小円の直径が 1.2 寸のとき，中円の直径はいかほどか。

注：この「問」には難点が多い。(1) 図には小円が見当たらない（図に示したところにある）。(2) 小円の直径が 1.2 寸，「答」の中円の直径が 5.7 寸有奇というのも不適切な数値である。そこで，「大円の直径が 10 寸」のみを条件として正しい解を求める。

正方形の一辺の長さは，大円の直径と同じである。
正方形の一辺の長さを a
大円の半径と中心座標を r1, (r1, r1); a = 2r1
中円の半径と中心座標を r2, (a - r2, r2), (r2, a - r2)
小円の半径と中心座標を r3, (x3, r3), (r3, x3), (a - x3, a - r3), (a - r3, a - x3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive,
     r3::positive, x3::positive;
a = 2r1
eq1 = (r1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq2 = r3/x3 - r2/(a - r2)
eq3 = dist2(0, 0, a, a, a - r2, r2, r2)
res = solve([eq1, eq2, eq3], (r2, r3, x3))[2]

   (r1*(2 - sqrt(2)), r1*(-3*sqrt(2) - 2*sqrt(4 - 2*sqrt(2)) + 2*sqrt(2 - sqrt(2)) + 5), r1*(-2*sqrt(2 - sqrt(2)) - 1 + 2*sqrt(2)))

大円の直径（正方形の一辺の長さ）が 10 のとき，中円の直径は 5.85786437626905，小円の直径は 1.23308641756286 である。

res[1](r1 => 10/2).evalf() * 2 |> println
res[2](r1 => 10/2).evalf() * 2 |> println

   5.85786437626905
   1.23308641756286

function draw(r1, r3, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 2r1
   (r2, r3, x3) = (r1*(2 - sqrt(2)), r1*(-3*sqrt(2) - 2*sqrt(4 - 2*sqrt(2)) + 2*sqrt(2 - sqrt(2)) + 5), r1*(-2*sqrt(2 - sqrt(2)) - 1 + 2*sqrt(2)))
   @printf("大円の直径（正方形の一辺の長さ）が %g のとき，中円の直径は %g，小円の直径は %g である。\n", 2r1, 2r2, 2r3)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
   circle(r1, r1, r1, :blue)
   circle(a - r2, r2, r2)
   circle(r2, a - r2, r2)
   circle(x3, r3, r3, :magenta)
   circle(r3, x3, r3, :magenta)
   circle(a - x3, a - r3, r3, :magenta)
   circle(a - r3, a - x3, r3, :magenta)
   segment(0, 0, a, a, :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(r1, r1, "大円:r1,(r1,r1)", :blue, :left, delta=-delta/2)
       point(a - r2, r2, "中円:r2,(a-r2,r2)", :red, :center, delta=-delta/2)
       point(x3, r3, " 小円:r3,(x3,r3)", :magenta, :left, :vcenter)
       point(2r1, 0, " 2r1", :green, :left, :bottom, delta=delta/2)
       point(0, 2r1, " 2r1", :green, :left, :bottom, delta=delta/2)
   end  
end;

draw(10/2, 1.2/2, true)

算額（その1327）

算額（その1327）

六三 羽生市須影 八幡神社 慶応元年(1865)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円5個，外円，菱形4個
#Julia, #SymPy, #算額, #和算

外円の中に，合同な菱形 4 個と等円 4 個を容れる。等円の直径が 1.5 寸のとき，外円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x, y); y = R/2
菱形の頂点座標を (a, b), (a, 3b); a = sqrt(R^2 - (R - b)^2), b = R/4
とおき，以下の連立方程式を解く。

include("julia-source.txt");

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

   (14*r/3, 2*sqrt(2)*r)

外円の半径 R は，等円の半径 r の14/3 倍である。
等円の直径が 1.5 寸のとき，外円の直径は 1.5*14/3 = 7 寸である。

すべてのパラメータは以下のとおりである。

   r = 0.75;  R = 3.5;  x = 2.12132;  y = 1.75;  a = 2.31503;  b = 0.875

function draw(r, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, x) = (14*r/3, 2*sqrt(2)*r)
   b = R/4
   a = sqrt(R^2 - (R - b)^2)
   y = R/2
   @printf("等円の直径が %g のとき，外円の直径は %g である。\n", 2r, 2R)
   @printf("r = %g;  R = %g;  x = %g;  y = %g;  a = %g;  b = %g\n", r, R, x, y, a, b)
   plot([0, a, -a, a, -a, 0, a, -a, a, -a, 0],
       [R, 3b, b, -b, -3b, -R, -3b, -1b, b, 3b, R], color=:green, lw=0.5)
   circle(0, 0, R, :blue)
   circle4(x, y, r)
   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, b, " (a,b)", :green, :left, :vcenter)
       point(a, 3b, " (a,3b)", :green, :left, :vcenter)
       point(x, 2b, "等円:r,(x,2b)", :red, :center, delta=-delta/2)
       point(0, R, "R", :blue, :center, :bottom, delta=delta/2)
   end  
end;

draw(1.5/2, true)

 

算額（その1326）

算額（その1326）

六三 羽生市須影 八幡神社 慶応元年(1865)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：正方形，斜線2本
#Julia, #SymPy, #算額, #和算

正方形の相対する頂点から対辺へ 2 本の等斜を引く。斜線と正方形の辺でできる 2 個の三角形の面積（黒積）が最大になるのはどのようなときか。

正方形の一辺の長さを a
等斜の両端の座標を [(0, a), (a, b)], [(a, 0), (b, a)]
等斜の交点を (x, y) = (a^2/(2*a - b), a^2/(2*a - b))
面積を S = (a - x)*b
とおく。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, x::positive, y::positive;
(x, y) = intersection(a, 0, b, a, 0, a, a, b)

   (a^2/(2*a - b), a^2/(2*a - b))

面積 S は，b により変化する。0 ≦ b ≦ 10 で，b が 6 前後のとき，最大値 15 程度になる。

S = (a - x)*b
S |> println

   b*(-a^2/(2*a - b) + a)

pyplot(size=(250, 120), grid=false, aspectratio=:none, label="")plot(S(a => 10), xlims=(0,10), xlabel="b", ylabel="S")

S を b で偏微分し，接線の傾きを表す式を得る。

diff_S = diff(S, b)
diff_S |> println

   -a^2*b/(2*a - b)^2 - a^2/(2*a - b) + a

S が極大値を取るとき，接線が 0 になるので，方程式を解く。
b = a*(2 - √2) のときに接線の傾きが 0 になり， S に代入すると最大値が得られる。
b = 5.85786437626905 のとき，面積が最大値 17.1572875253810 になる。

max_at = solve(diff_S, b)[1]
max_at |> println

   a*(2 - sqrt(2))

b0 = max_at(a => 10).evalf()
b0 |> println

   5.85786437626905

S(a =>10, b => b0) |> println

   17.1572875253810

function draw(a, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = a*(2 - √2)
   x = a^2/(2*a - b)
   S = (a - x)*b
   @printf("正方形の一辺の長さが %g のとき，黒積は最大値 %g になる。\n", a, S)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
   plot!([0, a, a, b, 0], [a, b, 0, a, a], color=:gray70, seriestype=:shape, fillcolor=:gray70, lw=0.5)
   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, a, " a", :green, :left, :bottom, delta=delta/2)
       point(a, b, " (a,b)", :green, :left, :vcenter)
       point(b, a, " (b,a)", :green, :center, :bottom, delta=delta/2)
       point(a^2/(2*a - b), a^2/(2*a - b), "(a^2/(2*a-b),a^2/(2*a-b))  ", :green, :right, :vcenter)
   end  
end;

draw(10, true)

算額（その2125）

算額（その2125）

一八 大里郡岡部村岡 稲荷社 文化14年(1817)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円5個，外円，弦3本
#Julia, #SymPy, #算額, #和算

外円の中に 3 本の弦を引き，区画された領域に甲円 1 個，乙円 1 個，丙円 2 個を容れる。甲円の直径が 5 寸，丙円の直径は甲円の直径の 1/3 のとき，乙円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, r1 - R)
乙円の半径と中心座標を r2, (0, R - r2)
丙円の半径と中心座標を r3, (x3, R - 2r2 + r3)
弦と外円の交点座標を (x01, sqrt(R^2 - x01^2)), (x02, sqrt(R^2 - x02^2)), (x03, sqrt(R^2 - x03^2))
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive,
     r3::positive, x3::positive, x01::positive, x02::positive;
R = r1 + r2
r3 = r1/3
eq1 = dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), 0, R - r2) - r2^2
eq2 = dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), 0, r1 - R) - r1^2
eq3 = dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), x3, R - 2r2 + r3) - r3^2
eq4 = x3^2 + (R - 2r2 + r3)^2 - (R - r3)^2;
# res = solve([eq1, eq2, eq3, eq4], (r2, x3, x01, x02))

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)
   (r2, x3, x01, x02) = u
   return [
       -r2^2 + (-x01 - (-x01 + x02)*(-x01*(-x01 + x02) + (r1 - sqrt(-x01^2 + (r1 + r2)^2))*(-sqrt(-x01^2 + (r1 + r2)^2) - sqrt(-x02^2 + (r1 + r2)^2)))/((-x01 + x02)^2 + (-sqrt(-x01^2 + (r1 + r2)^2) - sqrt(-x02^2 + (r1 + r2)^2))^2))^2 + (r1 - sqrt(-x01^2 + (r1 + r2)^2) - (-x01*(-x01 + x02) + (r1 - sqrt(-x01^2 + (r1 + r2)^2))*(-sqrt(-x01^2 + (r1 + r2)^2) - sqrt(-x02^2 + (r1 + r2)^2)))*(-sqrt(-x01^2 + (r1 + r2)^2) - sqrt(-x02^2 + (r1 + r2)^2))/((-x01 + x02)^2 + (-sqrt(-x01^2 + (r1 + r2)^2) - sqrt(-x02^2 + (r1 + r2)^2))^2))^2,  # eq1
       -r1^2 + (-x01 - (-x01 + x02)*(-x01*(-x01 + x02) + (-r2 - sqrt(-x01^2 + (r1 + r2)^2))*(-sqrt(-x01^2 + (r1 + r2)^2) - sqrt(-x02^2 + (r1 + r2)^2)))/((-x01 + x02)^2 + (-sqrt(-x01^2 + (r1 + r2)^2) - sqrt(-x02^2 + (r1 + r2)^2))^2))^2 + (-r2 - sqrt(-x01^2 + (r1 + r2)^2) - (-x01*(-x01 + x02) + (-r2 - sqrt(-x01^2 + (r1 + r2)^2))*(-sqrt(-x01^2 + (r1 + r2)^2) - sqrt(-x02^2 + (r1 + r2)^2)))*(-sqrt(-x01^2 + (r1 + r2)^2) - sqrt(-x02^2 + (r1 + r2)^2))/((-x01 + x02)^2 + (-sqrt(-x01^2 + (r1 + r2)^2) - sqrt(-x02^2 + (r1 + r2)^2))^2))^2,  # eq2
       -r1^2/9 + (-x01 + x3 - (-x01 + x02)*((-x01 + x02)*(-x01 + x3) + (-sqrt(-x01^2 + (r1 + r2)^2) - sqrt(-x02^2 + (r1 + r2)^2))*(4*r1/3 - r2 - sqrt(-x01^2 + (r1 + r2)^2)))/((-x01 + x02)^2 + (-sqrt(-x01^2 + (r1 + r2)^2) - sqrt(-x02^2 + (r1 + r2)^2))^2))^2 + (4*r1/3 - r2 - sqrt(-x01^2 + (r1 + r2)^2) - ((-x01 + x02)*(-x01 + x3) + (-sqrt(-x01^2 + (r1 + r2)^2) - sqrt(-x02^2 + (r1 + r2)^2))*(4*r1/3 - r2 - sqrt(-x01^2 + (r1 + r2)^2)))*(-sqrt(-x01^2 + (r1 + r2)^2) - sqrt(-x02^2 + (r1 + r2)^2))/((-x01 + x02)^2 + (-sqrt(-x01^2 + (r1 + r2)^2) - sqrt(-x02^2 + (r1 + r2)^2))^2))^2,  # eq3
       x3^2 - (2*r1/3 + r2)^2 + (4*r1/3 - r2)^2,  # eq4
   ]
end;

r1 = 5/2
R = r1 + r2
r3 = r1/3
iniv = BigFloat[1.5, 2.89418, 1.2121, 2.90904]
res = nls(H, ini=iniv)

   ([1.5, 2.581988897471611, 1.2103072956898178, 2.9047375096555625], true)

甲円の直径が 5 寸のとき，乙円の直径は 3 寸である。

すべてのパラメータは以下のとおりである。

   r1 = 2.5;  r2 = 1.5;  r3 = 0.833333;  x3 = 2.58199;  x01 = 1.21031;  x02 = 2.90474;  R = 4

function draw(r1, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = r1/3
   (r2, x3, x01, x02) = res[1]
   R = r1 + r2
   @printf("甲円の直径が %g のとき，乙円の直径は %g である。\n", 2r1, 2r2)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  x01 = %g;  x02 = %g;  R = %g\n", r1, r2, r3, x3, x01, x02, R)
   plot()
   circle(0, 0, R)
   circle(0, R - r2, r2, :blue)
   circle(0, r1 - R, r1, :magenta)
   circle2(x3, R - 2r2 + r3, r3, :orange)
   y01 = sqrt(R^2 - x01^2)
   y02 = -sqrt(R^2 - x02^2)
   segment(x01, y01, x02, y02, :green)
   segment(-x01, y01, -x02, y02, :green)
   y03 = R - 2r2
   x03 = sqrt(R^2 - y03^2)
   segment(-x03, y03, x03, y03, :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(0, R, "R", :red, :center, :bottom, delta=delta/2)
       point(x01, y01, "(x01,y01)", :green, :left, :bottom, delta=delta/2)
       point(x02, y02, "(x02,y02)", :green, :left, delta=-delta/2)
       point(x03, y03, "(x03,y03) ", :green, :right, delta=-delta/2)
       point(0, r1 - R, "甲円:r1,(0,r1-R)", :magenta, :center, delta=-delta/2)
       point(0, R - r2, "乙円:r2,(0,R-r2)", :blue, :center, delta=-delta/2)
       point(x3, R - 2r2 + r3, "丙円:r3\n(x3,R-2r2+r3)", :orange, :center, :bottom, delta=delta/2)
   end  
end;

draw(5/2, true)

 

算額（その1324）

算額（その1324）

一八 大里郡岡部村岡 稲荷社 文化14年(1817)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円7個，斜線2本，直線
#Julia, #SymPy, #算額, #和算

水平線と 2 本の斜線で区切られた領域に，等円 6 個と容円 1 個を容れる。等円の直径が 1.6 寸のとき，容円の直径はいかほどか。

等円の半径と中心座標を r1, (-2r1, -r1), (0, -r1), (2r1, -r1), (x11, r1), (x12, y12), (x13, r1)
容円の半径と中心座標を r2, (x2, r2)
3 直線の交点座標を ((a, 0), (b, 0), (c, d)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, x11::negative, x12::positive,
     y12::positive, x13::positive, r2::positive, x2::positive,
     a::positive, b::negative, c::negative, d::positive;;
eq1 = dist2(a, 0, c, d, x11, r1, r1)
eq2 = dist2(a, 0, c, d, x12, y12, r1)
# eq3 = dist2(a, 0, c, d, x13, r1, r1)
eq3 = (x11 - x12)^2 + (y12 - r1)^2 - 4r1^2;
eq4 = dist2(a, 0, c, d, 2r1, -r1, r1)
eq5 = dist2(a, 0, c, d, x2, r2, r2)
eq6 = dist2(b, 0, c, d, -2r1, -r1, r1)
eq7 = dist2(b, 0, c, d, x13, r1, r1)
eq8 = dist2(b, 0, c, d, x2, r2, r2)
eq9 = dist2(b, 0, c, d, x12, y12, r1)
eq10 = r2*(a - x13) - r1*(a - x2);
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10], (x11, x12, y12, x13, r2, x2, a, b, c, d))

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)
   (x11, x12, y12, x13, r2, x2, a, b, c, d) = u
   return [
       d*(a^2*d - 2*a^2*r1 + 2*a*c*r1 - 2*a*d*x11 + 2*a*r1*x11 - 2*c*r1*x11 - d*r1^2 + d*x11^2),  # eq1
       a^2*d^2 - 2*a^2*d*y12 - a^2*r1^2 + a^2*y12^2 + 2*a*c*d*y12 + 2*a*c*r1^2 - 2*a*c*y12^2 - 2*a*d^2*x12 + 2*a*d*x12*y12 - c^2*r1^2 + c^2*y12^2 - 2*c*d*x12*y12 - d^2*r1^2 + d^2*x12^2,  # eq2
       -4*r1^2 + (-r1 + y12)^2 + (x11 - x12)^2,  # eq3
       d*(a^2*d + 2*a^2*r1 - 2*a*c*r1 - 4*a*d*r1 - 4*a*r1^2 + 4*c*r1^2 + 3*d*r1^2),  # eq4
       d*(a^2*d - 2*a^2*r2 + 2*a*c*r2 - 2*a*d*x2 + 2*a*r2*x2 - 2*c*r2*x2 - d*r2^2 + d*x2^2),  # eq5
       d*(b^2*d + 2*b^2*r1 - 2*b*c*r1 + 4*b*d*r1 + 4*b*r1^2 - 4*c*r1^2 + 3*d*r1^2),  # eq6
       d*(b^2*d - 2*b^2*r1 + 2*b*c*r1 - 2*b*d*x13 + 2*b*r1*x13 - 2*c*r1*x13 - d*r1^2 + d*x13^2),  # eq7
       d*(b^2*d - 2*b^2*r2 + 2*b*c*r2 - 2*b*d*x2 + 2*b*r2*x2 - 2*c*r2*x2 - d*r2^2 + d*x2^2),  # eq8
       b^2*d^2 - 2*b^2*d*y12 - b^2*r1^2 + b^2*y12^2 + 2*b*c*d*y12 + 2*b*c*r1^2 - 2*b*c*y12^2 - 2*b*d^2*x12 + 2*b*d*x12*y12 - c^2*r1^2 + c^2*y12^2 - 2*c*d*x12*y12 - d^2*r1^2 + d^2*x12^2,  # eq9
       -r1*(a - x2) + r2*(a - x13),  # eq10
       ]
end;

r1 = 1.6/2
iniv = BigFloat[1.91073, 0.42645, 1.38001, -2.32044, 0.5, -0.81688, 1.73992, -1.95295, -0.951, 1.09764]
res = nls(H, ini=iniv)

   ([1.915367495237846, 0.43503904608801747, 1.4071471671989124, -2.3010731608400214, 0.504131448185704, -0.8, 1.757683747618923, -1.9505365804200108, -0.933017057376002, 1.1035735835994562], true)

等円の直径が 1.6 寸のとき，容円の直径は 1.00826 寸である。

function draw(r1, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x11, x12, y12, x13, r2, x2, a, b, c, d) = res[1]
   @printf("等円の直径が %g のとき，容円の直径は %g である。\n", 2r1, 2r2)
   plot()
   circle(-2r1, -r1, r1)
   circle(0, -r1, r1)
   circle(2r1, -r1, r1)
   circle(x11, r1, r1)
   circle(x12, y12, r1)
   circle(x13, r1, r1)
   circle(x2, r2, r2, :blue)
   abline(b, 0, d/(c - b), 1.4b, 0.8r1, :green)
   abline(a, 0, d/(c - a), 1.5a, 1.5b, :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)
       segment(x13 - r1, 0, x11 + r1, 0, :green)
       point(a, 0, " a", :green, :left, :bottom, delta=0.1delta)
       point(b, 0, "  b", :green, :left, :bottom, delta=0.1delta)
       point(c, d, "(c,d)", :green, :left, :vcenter, delta=delta, deltax=2delta)
       point(-2r1, -r1, "等円:r1\n(-2r1,-r1)", :red, :center, delta=-delta)
       point(0, -r1, "等円:r1\n(0,-r1)", :red, :center, delta=-delta)
       point(2r1, -r1, "等円:r1\n(2r1,-r1)", :red, :center, delta=-delta)
       point(x11, r1, "等円:r1\n(x11,r1)", :red, :center, delta=-delta)
       point(x12, y12, "等円:r1\n(x12,y12)", :red, :center, delta=-delta)
       point(x13, r1, "等円:r1\n(x13,r1)", :red, :center, delta=-delta)
       point(x2, r2, "容円:r2\n(x2,r2)", :blue, :center, delta=-delta)
   end  
end;

draw(1.6/2, true)

 

算額（その1323）

算額（その1323）

一八 大里郡岡部村岡 稲荷社 文化14年(1817)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円3個，正方形，等脚台形
#Julia, #SymPy, #算額, #和算

正方形の中に等脚台形と甲円，乙円，丙円を容れる。甲円の直径が 3 寸，丙円の直径が 1 寸のとき，乙円の直径はいかほどか。

正方形の一辺の長さを a
台形の頂点座標を (0, 0), (0, b), (c, a), (a, d)
甲円の半径と中心座標を r1, (a - r1, a - r1)
乙円の半径と中心座標を r2, (a - r2, r2)
丙円の半径と中心座標を r3, (r3, a - r3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive,
     d::positive, r1::positive, r2::positive, r3::positive;
eq1 = (a - c)^2 + (a - d)^2 - b^2
eq2 = ((a - b) + c - 2r3)^2 - ((a - b)^2 + c^2)
eq3 = ((a - d) + (a - c) - 2r1)^2 - ((a - d)^2 + (a - c)^2)
eq4 = (a + d - 2r2)^2 - (a^2 + d^2)
eq5 = (a - b)/c - d/a;

まず，eq1, eq2, eq5 から a, b, d を求める。

res = solve([eq1, eq2, eq5], (a, b, d))[2]  # 2 of 2

   (c*(c^2 - 2*r3^2)/(c^2 - 4*c*r3 + 2*r3^2), (c^2 - 2*c*r3 + 2*r3^2)^2/((c - 2*r3)*(c^2 - 4*c*r3 + 2*r3^2)), 2*r3*(c^3 - c^2*r3 - 2*c*r3^2 + 2*r3^3)/(c^3 - 6*c^2*r3 + 10*c*r3^2 - 4*r3^3))

eq3 の a, d に代入して c を求める。

eq13 = eq3(a => res[1], d => res[3]) |> simplify |> numerator |> factor
ans_c = solve(eq13, c)[1]
ans_c |> println
ans_c(r1 => 3/2, r3 => 1/2) |> println

   2*r3*(r1 - r3)/(r1 - 2*r3)
   2.00000000000000

eq4 に，a, b, d, c の順に代入し，r2 を求める。

eq14 = eq4(a => res[1], b => res[2], d => res[3])(c => ans_c) |> simplify |> numerator |> factor
ans_r2 = solve(eq14, r2)[2]  # 2 of 2
ans_r2 |> println
ans_r2(r1 => 3/2, r3 => 1/2) |> println

   r3*(-r1^2 + 2*r3^2)/(r1^2 - 4*r1*r3 + 2*r3^2)
   3.50000000000000

r2 は r1, r3 の関数である。甲円の直径が 3 寸，丙円の直径が 1 寸のとき，乙円の直径は 7 寸である。

function draw(r1, r3, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   c = 2*r3*(r1 - r3)/(r1 - 2*r3)
   r2 = r3*(-r1^2 + 2*r3^2)/(r1^2 - 4*r1*r3 + 2*r3^2)
   (a, b, d) = (c*(c^2 - 2*r3^2)/(c^2 - 4*c*r3 + 2*r3^2), (c^2 - 2*c*r3 + 2*r3^2)^2/((c - 2*r3)*(c^2 - 4*c*r3 + 2*r3^2)), 2*r3*(c^3 - c^2*r3 - 2*c*r3^2 + 2*r3^3)/(c^3 - 6*c^2*r3 + 10*c*r3^2 - 4*r3^3))
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
   plot!([0, a, c, 0, 0], [0, d, a, b, 0], color=:orange, lw=0.5)
   circle(a - r1, a - r1, r1)
   circle(a - r2, r2, r2, :blue)
   circle(r3, a - r3, r3, :magenta)
   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, a, "(a,a) ", :green, :right, :bottom, delta=delta/2)
       point(a, d, "(a,d)  ", :green, :right, :vcenter)
       point(c, a, "(c,a)", :green, :center, :bottom, delta=delta/2)
       point(0, b, "b ", :green, :right, :vcenter)
       point(a - r1, a - r1, "甲円:r1\n(a-r1,a-r1)", :red, :center, :bottom, delta=delta/2)
       point(a - r2, r2, "乙円:r2,(a-r2,r2)", :blue, :center, :bottom, delta=delta/2)
       point(r3, a - r3, " 丙円:r3,(r3,a-r3)", :magenta, :left, delta=-3delta)
       xlims!(-5delta, a + 3delta)
   end  
end;

draw(3/2, 1/2, true)

 

算額（その1322）

算額（その1322）

一七 大里郡岡部村岡 稲荷社 文化13年(1816)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円5個，斜線2本
#Julia, #SymPy, #算額, #和算

直線上に，2 本の接線を共有する大円 3 個，中円 1 個，小円 1 個を描く。大円と小円の直径の差が 2 寸のとき，中円の直径はいかほどか。

上の図は，大円と小円の直径の差が 5 寸のときのものである。問題の通り大円と小円の直径の差が 2 寸のときのものは下の方に示すが，似ても似つかない図になる。

大円と小円の直径の差を K
大円の半径と中心座標を r1, (x1, r1), (0, y1)
中円の半径と中心座標を r2, (0, r2)
小円の半径と中心座標を r3, (0, y1 - r1 - r3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive, y1::positive,
     r2::positive, r3::positive, a::positive,
     K::positive;
eq1 = dist2(-a, 0, 0, r1, x1, r1, r1)
eq2 = dist2(-a, 0, 0, r1, 0, y1, r1)
eq3 = dist2(-a, 0, 0, r1, 0, r2, r2)
eq4 = dist2(-a, 0, 0, r1, 0, y1 - r1 - r3, r3)
eq5 = r2*(a + x1) - r1*a  # r2/a - r1/(a + x1)
eq6 = (2r1 - 2r3) - K;

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)
   (r1, x1, y1, r2, r3, a) = u
   return [
       r1^2*(-a^2 - r1^2 + x1^2),  # eq1
       -2*a^2*r1*y1 + a^2*y1^2 - r1^4,  # eq2
       r1*(a^2*r1 - 2*a^2*r2 - r1*r2^2),  # eq3
       4*a^2*r1^2 + 4*a^2*r1*r3 - 4*a^2*r1*y1 - 2*a^2*r3*y1 + a^2*y1^2 - r1^2*r3^2,  # eq4
       -a*r1 + r2*(a + x1),  # eq5
       -K + 2*r1 - 2*r3,  # eq6
   ]
end;
K = 2.0
iniv = BigFloat[1.5, 6.3, 3.1, 0.75, 0.023, 6.1]
#iniv = BigFloat[3, 4.3, 7.2, 1.25, 0.5, 3.1]  # K = 5 のときの初期値
res = nls(H, ini=iniv)

   ([1.0034409020547763, 8.597223980953418, 2.0137872878330065, 0.5, 0.003440902054776372, 8.538463944689585], true)

大円と小円円の直径の差が 2 寸のとき，中円の直径は 1 寸である。

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

   K = 2;  r1 = 1.00296;  x1 = 9.25092;  y1 = 2.01188;  r2 = 0.5;  r3 = 0.00296477;  a = 9.19639

function draw(K, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, x1, y1, r2, r3, a) = res[1]
   @printf("大円と小円円の直径の差が %g のとき，中円の直径は %g である。\n", K, 2r2)
   @printf("K = %g;  r1 = %g;  x1 = %g;  y1 = %g;  r2 = %g;  r3 = %g;  a = %g\n", K, r1, x1, y1, r2, r3, a) 
   plot()
   circle2(x1, r1, r1)
   circle(0, y1, r1)
   circle(0, r2, r2, :blue)
   circle(0, y1 - r1 - r3, r3, :green)
   abline(-a, 0, r1/a, -a, 1.2x1, :magenta)
   abline(a, 0, -r1/a, a, -1.2x1, :magenta)
   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, y1, "大円:r1,(0,y1)", :red, :center, delta=-delta/2)
       point(x1, r1, "大円:r1,(x1,r1)", :red, :center, delta=-delta/2)
       point(0, r2, " 中円:r2\n(0,r2)", :blue, :center, delta=-delta/2)
       point(0, y1 - r1 - r3, " 小円:r3,(0,y1-r1-r3)", :black, :left, :vcenter)
   end  
end;

draw(2, false)

算額（その1321）

算額（その1321）

一七 大里郡岡部村岡 稲荷社 文化13年(1816)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円4個，直角三角形，円の個数
#Julia, #SymPy, #算額, #和算

直角三角形の中に数個の等円を容れる。等円の直径と赤積が与えられたとき，等円の個数を求めよ。

注：図がカラーではないし，赤積が何なのかも明示されていないので，「赤積　= 直角三角形の面積から等円の面積（n 個分）を引いたもの」とする。しかし，そのようにしても，与える赤積の精度により解が求まらないこともあるかもしれない。そこで，まず，等円の個数が n 個のときに，鈎，股を求める一般解を求める方法から始める。

直角三角形の直角を挟む二辺の短い方を「鈎」，長い方を「股」
等円の個数を n
等円の半径と中心座標を r, (r, r), (r, 3r), ..., (x, r)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms n, r::positive,  x::positive,  鈎::positive, 股::positive;
eq1 = 4r^2 + (x - r)^2 - ((n - 2)*2r)^2
eq2 = 2r/(x - r) - 鈎/股
eq3 = dist2(0, 鈎, 股, 0, x, r, r)
res = solve([eq1, eq2, eq3], (x, 鈎, 股))[3]  # 3 of 4

   (r*(6*n^2 + 10*n*sqrt(n^2 - 4*n + 3) - 19*n - 25*sqrt(n^2 - 4*n + 3) + 4)/(5*n + 3*sqrt(n^2 - 4*n + 3) - 14), 2*r*(3*n^2 + 5*n*sqrt(n^2 - 4*n + 3) - 12*n - 14*sqrt(n^2 - 4*n + 3) + 9)/(n^2 + 3*n*sqrt(n^2 - 4*n + 3) - 4*n - 9*sqrt(n^2 - 4*n + 3) + 3), r*(n + 3*sqrt(n^2 - 4*n + 3) - 1))

鈎，股は以下の通りである。

res[2] |> println
res[3] |> println

   2*r*(3*n^2 + 5*n*sqrt(n^2 - 4*n + 3) - 12*n - 14*sqrt(n^2 - 4*n + 3) + 9)/(n^2 + 3*n*sqrt(n^2 - 4*n + 3) - 4*n - 9*sqrt(n^2 - 4*n + 3) + 3)
   r*(n + 3*sqrt(n^2 - 4*n + 3) - 1)

n = 4, 5, ... と変化させて，鈎，股，赤積（と後に導く n の推定値）を書き出す。

r = 1/2
for n = 4:20 # n = 4:449 にすれば全ての結果が表示される
   鈎 = 2*r*(3*n^2 + 5*n*sqrt(n^2 - 4*n + 3) - 12*n - 14*sqrt(n^2 - 4*n + 3) + 9)/(n^2 + 3*n*sqrt(n^2 - 4*n + 3) - 4*n - 9*sqrt(n^2 - 4*n + 3) + 3)
   股 = r*(n + 3*sqrt(n^2 - 4*n + 3) - 1)
   赤積 = 鈎*股/2 - n*pi*r^2
   推定されたn = round(赤積*0.8223 + 2.522, digits=0)
   @printf("n = %2d;  鈎 = %g;  股 = %g;  赤積 = %g;  推定値 n = %d\n", n, 鈎, 股, 赤積, 推定されたn)
   round(赤積*0.8223 + 2.522, digits=0) != n &&        println(n)
end

   n =  4;  鈎 = 2.36603;  股 = 4.09808;  赤積 = 1.70648;  推定値 n = 4
   n =  5;  鈎 = 2.20711;  股 = 6.24264;  赤積 = 2.9621;  推定値 n = 5
   n =  6;  鈎 = 2.1455;  股 = 8.30948;  赤積 = 4.20159;  推定値 n = 6
   n =  7;  鈎 = 2.11237;  股 = 10.3485;  赤積 = 5.43212;  推定値 n = 7
   n =  8;  鈎 = 2.09161;  股 = 12.3741;  赤積 = 6.65772;  推定値 n = 8
   n =  9;  鈎 = 2.07735;  股 = 14.3923;  赤積 = 7.88035;  推定値 n = 9
   n = 10;  鈎 = 2.06695;  股 = 16.4059;  赤積 = 9.10106;  推定値 n = 10
   n = 11;  鈎 = 2.05902;  股 = 18.4164;  赤積 = 10.3205;  推定値 n = 11
   n = 12;  鈎 = 2.05277;  股 = 20.4248;  赤積 = 11.539;  推定値 n = 12
   n = 13;  鈎 = 2.04772;  股 = 22.4317;  赤積 = 12.7567;  推定値 n = 13
   n = 14;  鈎 = 2.04356;  股 = 24.4374;  赤積 = 13.974;  推定値 n = 14
   n = 15;  鈎 = 2.04006;  股 = 26.4422;  赤積 = 15.1909;  推定値 n = 15
   n = 16;  鈎 = 2.03709;  股 = 28.4464;  赤積 = 16.4075;  推定値 n = 16
   n = 17;  鈎 = 2.03452;  股 = 30.4499;  赤積 = 17.6238;  推定値 n = 17
   n = 18;  鈎 = 2.03229;  股 = 32.4531;  赤積 = 18.8399;  推定値 n = 18
   n = 19;  鈎 = 2.03033;  股 = 34.4558;  赤積 = 20.0558;  推定値 n = 19
   n = 20;  鈎 = 2.02859;  股 = 36.4583;  赤積 = 21.2716;  推定値 n = 20

以上より，「等円の直径が 1 で，赤積 < 542.3553863336211 であれば，赤積を 0.8223 倍して 2.522 を加え，小数点以下を四捨五入すると等円の個数が得られる」
なお，赤積の二乗項まで使えば，予測可能領域は広がる。

n = 449
鈎 = 2*r*(3*n^2 + 5*n*sqrt(n^2 - 4*n + 3) - 12*n - 14*sqrt(n^2 - 4*n + 3) + 9)/(n^2 + 3*n*sqrt(n^2 - 4*n + 3) - 4*n - 9*sqrt(n^2 - 4*n + 3) + 3)
股 = r*(n + 3*sqrt(n^2 - 4*n + 3) - 1)
赤積 = 鈎*股/2 - n*pi*r^2
#@printf("n = %2d;  鈎 = %g;  股 = %g;  赤積 = %g\n", n, 鈎, 股, 赤積)
println(赤積, ";  ", round(赤積*0.8223 + 2.522, digits=0))

   542.3553863336211;  449.0

以下の描画プログラムは，精度の良くない赤積の値を与えると，それに最も近い赤積を与える n を求める。

r = 1/2, 赤積 = 9 を与えると「与えられた赤積 = 9;  等円の直径 = 2;  正確な赤積 = 6.82593;  n = 4」を表示する。

もし，等円の半径が 1/2 以外のときは，面積は半径の二乗に比例するので，換算して計算する。
r = 4, n = 10 のときは赤積は r = 1/2 のときの 8^2 倍なので，8^2 * 9.10106 = 582 を与えると，「与えられた赤積 = 582;  等円の直径 = 8;  正確な赤積 = 582.468;  n = 10」を表示する。

function draw(r0, 赤積0, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   f = 2r0
   赤積 = 赤積0 / f^2
   n = round(赤積*0.8223 + 2.522, digits=0)
   (x, 鈎, 股) = (r*(6*n^2 + 10*n*sqrt(n^2 - 4*n + 3) - 19*n - 25*sqrt(n^2 - 4*n + 3) + 4)/(5*n + 3*sqrt(n^2 - 4*n + 3) - 14), 2*r*(3*n^2 + 5*n*sqrt(n^2 - 4*n + 3) - 12*n - 14*sqrt(n^2 - 4*n + 3) + 9)/(n^2 + 3*n*sqrt(n^2 - 4*n + 3) - 4*n - 9*sqrt(n^2 - 4*n + 3) + 3), r*(n + 3*sqrt(n^2 - 4*n + 3) - 1))
   赤積 = 鈎*股/2 - n*pi*r^2
   string = @sprintf("与えられた赤積 = %g; 等円の直径 = %g\n正確な赤積 = %g;  n = %d", 赤積0, 2r0, 赤積*f^2, n)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:green, lw=0.5)
   circle(r, r, r)
   for i = 0:n - 2
       circle(r + (x - r)/(n - 2)*i, 3r - 2r/(n - 2)*i, r)
   end
   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, 鈎, string, :black, :left, :bottom, delta=delta, mark=false)
   end
end;

draw(4, 582, true)

 

算額（その1320）

算額（その1320）

九 武州崎玉郡騎西町 久伊豆神社 文化4年(1807)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円8個，接線4本，直線上
#Julia, #SymPy, #算額, #和算

直線上に接線 4 本で仕切られている区画に 8 個の円を容れる。乙円と丙円の直径がそれぞれ 44 寸，20 寸のとき，丁円の直径はいかほどか。

甲円の半径と中心座標を r1,(0, y1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (0, r3)
丁円の半径と中心座標を r4, (x4, r4), (x42, r2)
接線が通る 2 点の座標対を [(-a, 0), (x01, y01)], [(b, 0), (x02, y02)]
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms r1::positive, y1::positive, r2::positive, x2::positive,
     r3::positive, r4::positive, x4::positive, x42::positive,
     a::positive, b::positive, x01::positive, y01::positive,
     x02::positive, y02::positive;
eq1 = dist2(-a, 0, x01, y01, 0, y1, r1)
eq2 = dist2(-a, 0, x01, y01, x2, r2, r2)
# eq3 = dist2(-a, 0, x01, y01, 0, r3, r3)
# eq4 = dist2(-a, 0, x01, y01, x42, r2, r4)
eq3 = r3/a - r2/(a + x2)
eq4 = r4/(a - x4) - r3/a
eq5 = dist2(b, 0, x02, y02, 0, y1, r1)
eq6 = dist2(b, 0, x02, y02, x2, r2, r2)
eq7 = dist2(b, 0, x02, y02, 0, r3, r3)
eq8 = dist2(b, 0, x02, y02, x4, r4, r4)
eq9 = dist2(b, 0, x02, y02, x42, r2, r4)
eq10 = dist2(a, 0, -x01, y01, x2, r2, r2)
eq11 = dist2(a, 0, -x01, y01, x4, r4, r4)
eq12 = dist2(a, 0, -x01, y01, x42, r2, 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)
   (r1, y1, x2, r4, x4, x42, a, b, x01, y01, x02, y02) = u
   return [
-a^2*r1^2 + a^2*y01^2 - 2*a^2*y01*y1 + a^2*y1^2 - 2*a*r1^2*x01 - 2*a*x01*y01*y1 + 2*a*x01*y1^2 - r1^2*x01^2 - r1^2*y01^2 + x01^2*y1^2,  # eq1
y01*(-2*a^2*r2 + a^2*y01 - 2*a*r2*x01 - 2*a*r2*x2 + 2*a*x2*y01 - r2^2*y01 - 2*r2*x01*x2 + x2^2*y01),  # eq2
-r2/(a + x2) + r3/a,  # eq3
r4/(a - x4) - r3/a,  # eq4
-b^2*r1^2 + b^2*y02^2 - 2*b^2*y02*y1 + b^2*y1^2 + 2*b*r1^2*x02 + 2*b*x02*y02*y1 - 2*b*x02*y1^2 - r1^2*x02^2 - r1^2*y02^2 + x02^2*y1^2,  # eq5
y02*(-2*b^2*r2 + b^2*y02 + 2*b*r2*x02 + 2*b*r2*x2 - 2*b*x2*y02 - r2^2*y02 - 2*r2*x02*x2 + x2^2*y02),  # eq6
y02*(-2*b^2*r3 + b^2*y02 + 2*b*r3*x02 - r3^2*y02),  # eq7
y02*(-2*b^2*r4 + b^2*y02 + 2*b*r4*x02 + 2*b*r4*x4 - 2*b*x4*y02 - r4^2*y02 - 2*r4*x02*x4 + x4^2*y02),  # eq8
b^2*r2^2 - 2*b^2*r2*y02 - b^2*r4^2 + b^2*y02^2 - 2*b*r2^2*x02 + 2*b*r2*x02*y02 + 2*b*r2*x42*y02 + 2*b*r4^2*x02 - 2*b*x42*y02^2 + r2^2*x02^2 - 2*r2*x02*x42*y02 - r4^2*x02^2 - r4^2*y02^2 + x42^2*y02^2,  # eq9
y01*(-2*a^2*r2 + a^2*y01 - 2*a*r2*x01 + 2*a*r2*x2 - 2*a*x2*y01 - r2^2*y01 + 2*r2*x01*x2 + x2^2*y01),  # eq10
y01*(-2*a^2*r4 + a^2*y01 - 2*a*r4*x01 + 2*a*r4*x4 - 2*a*x4*y01 - r4^2*y01 + 2*r4*x01*x4 + x4^2*y01),  # eq11
a^2*r2^2 - 2*a^2*r2*y01 - a^2*r4^2 + a^2*y01^2 + 2*a*r2^2*x01 - 2*a*r2*x01*y01 + 2*a*r2*x42*y01 - 2*a*r4^2*x01 - 2*a*x42*y01^2 + r2^2*x01^2 + 2*r2*x01*x42*y01 - r4^2*x01^2 - r4^2*y01^2 + x42^2*y01^2,  # eq12
   ]
end;

(r2, r3) = (44/2, 20/2)
iniv = BigFloat[28, 55, 40, 5.5, 15, 10, 33, 6.6, 56, 59, 39, 77]
res = nls(H, ini=iniv)

   ([27.5, 55.0, 39.7994974842648, 5.5, 14.9248115565993, 9.9498743710662, 33.166247903554, 6.6332495807108, 54.652326389566944, 58.25225211084647, 39.43354494756409, 77.70448053191785], true)

乙円と丙円の直径がそれぞれ 44 寸，20 寸のとき，丁円の直径は 2*5.5 = 11 寸である。

すべてのパラメータは以下のとおりである。

    r2 = 22;  r3 = 10;  r1 = 27.5;  y1 = 55;  x2 = 39.7995;  r4 = 5.5;  x4 = 14.9248;  x42 = 9.94987;  a = 33.1662;  b = 6.63325;  x01 = 54.6523;  y01 = 58.2523;  x02 = 39.4335;  y02 = 77.7045

function draw(r2, r3, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, y1, x2, r4, x4, x42, a, b, x01, y01, x02, y02) = [27.5, 55.0, 39.7994974842648, 5.5, 14.9248115565993, 9.9498743710662, 33.166247903554, 6.6332495807108, 54.652326389566944, 58.25225211084647, 39.43354494756409, 77.70448053191785]
   plot()
   circle(0, y1, r1)
   circle2(x2, r2, r2, :blue)
   circle(0, r3, r3, :green)
   circle2(x4, r4, r4, :magenta)
   circle2(x42, r2, r4, :magenta)
   segment(-a, 0, x01, y01)
   segment(b, 0, x02, y02)
   segment(a, 0, -x01, y01)
   segment(-b, 0, -x02, y02)
   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(x01, y01, " (x01,y01)", :black, :left, :vcenter)
       point(x02, y02, " (x02,y02)", :black, :left, :vcenter)
       point(-a, 0, "-a", :black, :center, delta=-1.5delta)
       point(b, 0, "b", :black, :center, delta=-1.5delta)
       point(0, y1, "甲円:r1,(0,y1)", :red, :center, delta=-delta/2)
       point(x2, r2, "乙円:r2,(x2,r2)", :blue, :center, delta=-delta/2)
       point(0, r3, "丙円:r3\n(0,r3)", :green, :center, delta=-delta/2)
       point(x4, r4, "丁円:r4,(x4,r4)", :black, :left, delta=-delta/2)
       point(x42, r2, "丁円:r4,(x42,r2)", :black, :left, :bottom, delta=delta)
       plot!(xlims=(-(x2 + r2 + 3delta), x2 + r2 + 13delta), ylims=(-5delta, y1 + r1 + 2delta))
   end
end;

draw(44/2, 20/2, true)

 

算額（その1319）

算額（その1319）

八 武州足立郡桶川宿 不動堂 文化3年(1806)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円8個
#Julia, #SymPy, #算額, #和算

大円 4 個が交わっている隙間に天円 2 個，等円 2 個を容れる。天円の直径が 15 寸のとき，等円の直径はいかほどか。

大円の半径と中心座標を r1, (0, r1/2), (√3r1/2, 0)
天円の半径と中心座標を r2, (0, 3r1/2 - r2)
等円の半径と中心座標を r3, (√3r1/2 + r1 - r3, 0)
とおき，以下の方程式を解く。
まず eq1 を解いて r1 を求め，解を eq2 に代入して解いて r3 を求める。
（連立して解くと SymPy では解を求めることができない）

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive;
eq1 = (√Sym(3)*r1/2)^2 + (3r1/2 - r2)^2 - (r1 + r2)^2
eq2 = (√Sym(3)*r1/2 + r1 - r3)^2 + (r1/2)^2 - (r1 + r3)^2;

大円の半径 r1 は，天円の半径 r2 の 5/2 倍である。
天円の直径が 15 寸のとき，大円の直径は 15*5/2 = 57.5 である。

ans_r1 = solve(eq1, r1)[1]
ans_r1 |> println
ans_r1(r2 => 15/2) * 2 |> println

   5*r2/2
   37.5000000000000

等円の半径 r3 は，天円の半径 r2 の 5(1 + 3√3)/26 倍である。
天円の直径が 15 寸のとき，等円の直径は 15*5(1 + 3√3)/26 = 17.873516603961438 である。

eq12 = eq2(r1 => ans_r1)
ans_r3 = solve(eq12, r3)[1] |> factor
ans_r3 |> println
ans_r3(r2 => 15/2).evalf() * 2 |> println

   5*r2*(1 + 3*sqrt(3))/26
   17.8735166039614

「術」は以下のようになっており，簡約化すると同じ関係式になる。

@syms 天径, d
子 = √Sym(3) + 1
等円径 = 子/2(子 + 3)*5天径
等円径 |> println
等円径(天径 => 15).evalf() |> println

   5*天径*(1 + sqrt(3))/(2*sqrt(3) + 8)
   17.8735166039614

apart(等円径, d) |> simplify |> println

   5*天径/26 + 15*sqrt(3)*天径/26

function draw(r2, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 5r2/2
   r3 = 5r2*(1 + 3√3)/26
   @printf("天円の直径が %g のとき，等円の直径は %g，大円の直径は %g である。\n", 2r2, 2r3, 2r1)
   plot()
   circle2(√3r1/2, 0, r1)
   circle22(0, r1/2, r1)
   circle22(0, 3r1/2 - r2, r2, :green)
   circle2(√3r1/2 + r1 - r3, 0, r3, :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)
       point(0, r1/2, " r1/2", :red, :left, :bottom, delta=delta/4)
       point(√3r1/2, 0, "√3r1/2 ", :red, :right, delta=-delta/2)
       point(0, 3r1/2 - r2, "天円:r2\n(0,3r1/2-r2)", :green, :center, delta=delta/2)
       point(√3r1/2 + r1 - r3, 0, "等円:r3\n(√3r1/2+r1-r3, 0)", :blue, :center, :bottom, delta=delta/2)
   end
end;

draw(15/2, true)

算額（その1318）

算額（その1318）

八六 加須市多聞寺 愛宕神社 明治13年(1880)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円7個，外円，円弧，最大値
#Julia, #SymPy, #算額, #和算

外円の中に，外円と同じ半径を持つ円弧 1 個と円弧に接する弦と，甲円，乙円，丙円を 2 個ずつ容れる。乙円の直径が 1 寸のとき，丙円の直径が取りうる最大値はいかほどか。
注：弦の位置（弦と y 軸の交点の y 座標）により丙円の大きさが変化する。

弦と y 軸の交点の y 座標を y
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (r1, 0); r1 = R/2
乙円の半径と中心座標を r2, (0, R - r2)
丙円の半径と中心座標を r3, (x3, y - r3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms y, R::positive, r1::positive,
     r2::positive, r3::positive, x3::positive;
R = 2r1
eq1 = r1^2 + (R - r2)^2 - (r1 + r2)^2
eq2 = x3^2 + (R - r3)^2 - (r3 + R)^2
eq3 = x3^2 + (y - r3)^2 - (R - r3)^2 |> expand
res = solve([eq1, eq2, eq3], (r3, x3, r1))[1]

   ((-3*r2 + y)*(3*r2 + y)/(2*(-9*r2 + y)), -sqrt(6)*sqrt(r2)*sqrt((-3*r2 + y)*(3*r2 + y)/(-9*r2 + y)), 3*r2/2)

丙円の半径は r3 = (-3*r2 + y)*(3*r2 + y)/(2*(-9*r2 + y)) となり，y と r2 の関数である。
r2 は事前に与えられるので，y の取る値で r3 が変化する。

たとえば，乙円の直径が 1 寸のとき(注2)，-1.5 ≦ y ≦ 1.5 で，y が 0.0 〜 0.5 の区間でr3 が最大値を取る。

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="")
plot(res[1](r2 => 1/2), xlims=(-1.5,1.5), xlabel="y", ylabel="r3")

r3 が最大値を取るときの y，および，そのときの最大値を求める。

r3_dash = diff(res[1], y) |> factor;  # 導関数
r3_dash |> println

   (9*r2^2 - 18*r2*y + y^2)/(2*(-9*r2 + y)^2)

r3_dash = 0 を解いて y を求める。

ans_y = solve(r3_dash, y)[1]
ans_y |> println

   3*r2*(3 - 2*sqrt(2))

甲円の半径 r2 が 1/2 のときの丙円の半径 r3 を求める。

ans_r3 = res[1](y => ans_y) |> simplify
ans_r3 |> println
ans_r3(r2 => 1/2).evalf() |> println

   3*r2*(3 - 2*sqrt(2))
   0.257359312880715

丙円の中心の x 座標を求める。
res[2] は負の値を取るが，左右いずれの甲円でも問題を解く上では無関係である。

ans_x3 = res[2](y => ans_y) |> simplify
ans_x3 |> println
ans_x3(r2 => 1/2).evalf() |> println

   -3*2^(3/4)*r2*sqrt(-4 + 3*sqrt(2))
   -1.24264068711929

最後に甲円の半径を求める。

ans_r1 = res[3]
ans_r1 |> println
ans_r1(r2 => 1/2) |>  println

   3*r2/2
   0.750000000000000

図を描くために，甲円と円弧の交点座標を求める。

@syms, x0, y0
eq4 = x0^2 +  y0^2 - R^2
eq5 = x0^2 + (y0 - y + R)^2 - R^2
res2 = solve([eq4, eq5], (x0, y0))[1]

   (-sqrt(-(-6*r1 + y)*(2*r1 + y))/2, (-2*r1 + y)/2)

まとめ

乙円の半径 r2 が与えられれば，乙円の半径が最大になるときの弦の位置 y が決まる。y = 3r2*(3 - 2√2)
丙円の半径 r3 は r3 = 3r2*(3 - 2√2) である。丙円の半径だけを求めるのならば，y は求める必要はない。

乙円の直径が 1 寸のとき，丙円の直径は 3(3 - 2√2) = 0.5147186257614291 寸である。

なお，下に示す算額の図は，丙円が最大になったときの図ではない。問題を解こうとする人を欺くためにわざとこのような図を描いたのではないかと邪推する。

function draw(r2, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   y = 3*r2*(3 - 2*sqrt(2))
   r3 = 3*r2*(3 - 2*sqrt(2))
   x3 = 3*2^(3/4)*r2*sqrt(-4 + 3*sqrt(2))
   r1 = 3*r2/2
   R = 2r1
   x = sqrt(R^2 - y^2)
   (x0, y0) = (-sqrt(-(-6*r1 + y)*(2*r1 + y))/2, (-2*r1 + y)/2)
   @printf("乙円の直径が %g，y が %g のとき，丙円の直径は最大値 %g になる。\n", 2r2, y, 2r3)
   plot()
   circle(0, 0, R, :green)
   circle2(r1, 0, r1)
   circle22(0, R - r2, r2, :blue)
   circle2(x3, y - r3, r3, :magenta)
   segment(-x, y, x, y)
   θ = atand(y0, x0)
   circle(0, y - R, R, :green, beginangle=-θ, endangle=180+θ)
   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, "R", :green, :center, :bottom, delta=delta)
       point(0, y, "y", :black, :center, :bottom, delta=delta)
       point(r1, 0, "甲円:r1,(r1,0)", :red, :right, delta=-delta, deltax=3delta)
       point(0, R - r2, "乙円:r2\n(0,R-r2)", :blue, :center, delta=-delta)
       point(0, r2 - R, "乙円:r2\n(0,r2-R)", :blue, :center, :bottom, delta=delta)
       point(x3, y - r3, "丙円:r3,(y-r3,0)", :magenta, :right, :bottom, delta=delta, deltax=3delta)
       point(0, y - R, "円弧:R,(0,y-R)", :green, :center, delta=-delta)
   end
end;

draw(1/2, true)

算額（その1317）

算額（その1317）

百四十七 群馬県甘楽郡妙義町下高田 高太神社 大正12年(1923)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，等脚台形
#Julia, #SymPy, #算額, #和算

等脚台形の中に大円 2 個，小円 2 個を容れる。台形の高さが 20 寸，大円の直径が 14.4 寸のとき，小円の直径はいかほどか。

等脚台形の上底，下底の長さをそれぞれ 2b, 2a 
等脚台形の斜辺を延長してできる二等辺三角形の高さを h0
大円の半径と中心座標を r1, (r1, r1)
小円の半径と中心座標を r2, (r2, h - r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, h::positive, h0::positive,
     r1::positive, r2::positive, y2::positive;
eq1 = r1/(h0 - r1) - r2/(h0 - (h - r2))
eq2 = b/(h0 - h) - a/h0
eq3 = dist2(0, h0, a, 0, r1, r1, r1)
eq4 = (r1 - r2)^2 + (h - r2 - r1)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3, eq4], (r2, a, b, h0))[2]

   ((-2*h^(5/2)*sqrt(r1) + 12*h^(3/2)*r1^(3/2) - 16*sqrt(h)*r1^(5/2) + h^3 - 5*h^2*r1 + 2*h*r1^2 + 8*r1^3)/(h^2 - 6*h*r1 + 8*r1^2), 4*(h^(17/2)*r1^(3/2) - 26*h^(15/2)*r1^(5/2) + 249*h^(13/2)*r1^(7/2) - 1192*h^(11/2)*r1^(9/2) + 3160*h^(9/2)*r1^(11/2) - 4768*h^(7/2)*r1^(13/2) + 3984*h^(5/2)*r1^(15/2) - 1664*h^(3/2)*r1^(17/2) + 256*sqrt(h)*r1^(19/2) - h^8*r1^2 + 26*h^7*r1^3 - 249*h^6*r1^4 + 1192*h^5*r1^5 - 3160*h^4*r1^6 + 4768*h^3*r1^7 - 3984*h^2*r1^8 + 1664*h*r1^9 - 256*r1^10)/(h^9 - 28*h^8*r1 + 301*h^7*r1^2 - 1690*h^6*r1^3 + 5544*h^5*r1^4 - 11088*h^4*r1^5 + 13520*h^3*r1^6 - 9632*h^2*r1^7 + 3584*h*r1^8 - 512*r1^9), 4*(-h^(9/2)*sqrt(r1) + 13*h^(7/2)*r1^(3/2) - 52*h^(5/2)*r1^(5/2) + 68*h^(3/2)*r1^(7/2) - 16*sqrt(h)*r1^(9/2) - h^4*r1 + h^3*r1^2 + 24*h^2*r1^3 - 52*h*r1^4 + 16*r1^5)/(h^4 - 18*h^3*r1 + 84*h^2*r1^2 - 120*h*r1^3 + 32*r1^4), (sqrt(h)*r1*(h - 6*r1) + 2*r1^(3/2)*(-h + 2*r1))/(sqrt(h)*(h - 4*r1)))

res[1] |> println

   (-2*h^(5/2)*sqrt(r1) + 12*h^(3/2)*r1^(3/2) - 16*sqrt(h)*r1^(5/2) + h^3 - 5*h^2*r1 + 2*h*r1^2 + 8*r1^3)/(h^2 - 6*h*r1 + 8*r1^2)

小円の半径 r2 は，h, r1 を用いてい以下のように計算できる。

f(h, r1) = (-2*h^(5/2)*sqrt(r1) + 12*h^(3/2)*r1^(3/2) - 16*sqrt(h)*r1^(5/2) + h^3 - 5*h^2*r1 + 2*h*r1^2 + 8*r1^3)/(h^2 - 6*h*r1 + 8*r1^2)
f(20, 14.4/2) * 2

   6.399999999999852

台形の高さが 20 寸，大円の直径が 14.4 寸のとき，小円の直径は 6.4 寸である。

すべてのパラメータは以下のとおりである。

   h = 20;  r1 = 7.2;  r2 = 3.2;  a = 24.6857;  b = 4.51765;  h0 = 24.48

function draw(h, r1, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, a, b, h0) = ((-2*h^(5/2)*sqrt(r1) + 12*h^(3/2)*r1^(3/2) - 16*sqrt(h)*r1^(5/2) + h^3 - 5*h^2*r1 + 2*h*r1^2 + 8*r1^3)/(h^2 - 6*h*r1 + 8*r1^2), 4*(h^(17/2)*r1^(3/2) - 26*h^(15/2)*r1^(5/2) + 249*h^(13/2)*r1^(7/2) - 1192*h^(11/2)*r1^(9/2) + 3160*h^(9/2)*r1^(11/2) - 4768*h^(7/2)*r1^(13/2) + 3984*h^(5/2)*r1^(15/2) - 1664*h^(3/2)*r1^(17/2) + 256*sqrt(h)*r1^(19/2) - h^8*r1^2 + 26*h^7*r1^3 - 249*h^6*r1^4 + 1192*h^5*r1^5 - 3160*h^4*r1^6 + 4768*h^3*r1^7 - 3984*h^2*r1^8 + 1664*h*r1^9 - 256*r1^10)/(h^9 - 28*h^8*r1 + 301*h^7*r1^2 - 1690*h^6*r1^3 + 5544*h^5*r1^4 - 11088*h^4*r1^5 + 13520*h^3*r1^6 - 9632*h^2*r1^7 + 3584*h*r1^8 - 512*r1^9), 4*(-h^(9/2)*sqrt(r1) + 13*h^(7/2)*r1^(3/2) - 52*h^(5/2)*r1^(5/2) + 68*h^(3/2)*r1^(7/2) - 16*sqrt(h)*r1^(9/2) - h^4*r1 + h^3*r1^2 + 24*h^2*r1^3 - 52*h*r1^4 + 16*r1^5)/(h^4 - 18*h^3*r1 + 84*h^2*r1^2 - 120*h*r1^3 + 32*r1^4), (sqrt(h)*r1*(h - 6*r1) + 2*r1^(3/2)*(-h + 2*r1))/(sqrt(h)*(h - 4*r1)))
   @printf("台形の高さが %g，大円の直径が %g のとき，小円の直径は %g である。\n", h, 2r1, 2r2)
   @printf("h = %g;  r1 = %g;  r2 = %g;  a = %g;  b = %g;  h0 = %g\n", h, r1, r2, a, b, h0)
   plot([a, 0, -a, a], [0, h0, 0, 0], color=:green, lw=0.5)
   segment(-b, h, b, h, :green)
   circle2(r1, r1, r1)
   circle2(r2, h - r2, r2, :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)
       point(r1, r1, "大円:r1,(r1,r1)", :red, :center, delta=-delta/2)
       point(r2, h - r2, " 小円:r2,(r2,h-r2)", :blue, :left, :vcenter)
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
       point(0, h, " h", :green, :left, :bottom, delta=delta/2)
       point(0, h0, " h0", :green, :left, :bottom, delta=delta/2)
       point(b, h, " (b,h)", :green, :left, :vcenter)
   end
end;

draw(20, 14.4/2, true)

 

算額（その1316）

算額（その1316）

百四十七 群馬県甘楽郡妙義町下高田 高太神社 大正12年(1923)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，外円，正三角形
#Julia, #SymPy, #算額, #和算

外円の中に大円 2 個，小円 1 個，直角三角形 1 個を容れる。正三角形の一辺の長さが 7.5 寸，小円の直径が 4.755 寸のとき，大円の直径はいかほどか。

正三角形の一辺の長さを 2a，頂点と円の接点座標を (a, x0, y0);x0 = a, y0 = -sqrt(R^2 - a^2)
外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (x1, y1); x1 = r1, y1 = √Sym(3)a + y0
小円の半径と中心座標を r2, (0, R - r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, R::positive, r1::positive,
     x1::positive, y1::positive, r2::positive;
x0 = a
y0 = -sqrt(R^2 - a^2)
x1 = r1
y1 = √Sym(3)a + y0
eq1 = x1^2 + y1^2 - (R - r1)^2
eq2 = R - 2r2 - √Sym(3)a - y0
res = solve([eq1, eq2], (r1, R))[1]

   (a*(-3*a^3 - 4*sqrt(3)*a^2*r2 - 4*a*r2^2 + 3*a*sqrt(a^4 + 4*sqrt(3)*a^3*r2 + 16*a^2*r2^2 + 8*sqrt(3)*a*r2^3 + 4*r2^4) + 2*r2*sqrt(3*a^4 + 12*sqrt(3)*a^3*r2 + 48*a^2*r2^2 + 24*sqrt(3)*a*r2^3 + 12*r2^4))/(2*(sqrt(3)*a^3 + 5*a^2*r2 + 3*sqrt(3)*a*r2^2 + 2*r2^3)), 2*(a^2 + sqrt(3)*a*r2 + r2^2)/(sqrt(3)*a + 2*r2))

正三角形の一辺の長さが 7.5 寸，小円の直径が 4.755 寸のとき，大円の直径は 5.89244501321687 寸である。

「答」には「大円径六寸二分五厘」と書いているが，それは多分，外円の半径である。

res[1](a => 7.5/2, r2 => 4.755/2).evalf() * 2 |> println

   5.89244501321687

f(a, r2) = a*(-3*a^3 - 4*sqrt(3)*a^2*r2 - 4*a*r2^2 + 3*a*sqrt(a^4 + 4*sqrt(3)*a^3*r2 + 16*a^2*r2^2 + 8*sqrt(3)*a*r2^3 + 4*r2^4) + 2*r2*sqrt(3*a^4 + 12*sqrt(3)*a^3*r2 + 48*a^2*r2^2 + 24*sqrt(3)*a*r2^3 + 12*r2^4))/(2*(sqrt(3)*a^3 + 5*a^2*r2 + 3*sqrt(3)*a*r2^2 + 2*r2^3))
f(7.5/2, 4.755/2) * 2

   5.892445013216866

function draw(a, r2, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, R) = (a*(-3*a^3 - 4*sqrt(3)*a^2*r2 - 4*a*r2^2 + 3*a*sqrt(a^4 + 4*sqrt(3)*a^3*r2 + 16*a^2*r2^2 + 8*sqrt(3)*a*r2^3 + 4*r2^4) + 2*r2*sqrt(3*a^4 + 12*sqrt(3)*a^3*r2 + 48*a^2*r2^2 + 24*sqrt(3)*a*r2^3 + 12*r2^4))/(2*(sqrt(3)*a^3 + 5*a^2*r2 + 3*sqrt(3)*a*r2^2 + 2*r2^3)), 2*(a^2 + sqrt(3)*a*r2 + r2^2)/(sqrt(3)*a + 2*r2))
   x1 = r1
   y1 = √3a - sqrt(R^2 - a^2)
   x0 = a
   y0 = -sqrt(R^2 - x0^2)
   @printf("正三角形の一辺の長さが %g，小円の直径が %g のとき，大円の直径は %g である。\n", 2a, 2r2, 2r1)
   @printf("a = %g;  r2 = %g;  R = %g;  r1 = %g;  x1 = %g;  y1 = %g;  x0 = %g;  y0 = %g\n", a, r2, R, r1, x1, y1, x0, y0)
   plot([a, 0, -a, a], [y0, √3a + y0, y0, y0], color=:green, lw=0.5)
   circle(0, 0, R, :magenta)
   circle2(x1, y1, r1)
   circle(0, R - r2, r2, :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)
       point(0, R, "R", :red, :center, :bottom, delta=delta/2)
       point(0, √3a + y0, "√3a+y0  ", :green, :right, delta=-delta/2)
       point(x1, y1, "大円:r1,(x1,y1)", :red, :center, delta=-delta)
       point(0, R - r2, "小円:r2,(0,R-r2)", :blue, :center, :bottom, delta=2delta)
       point(x0, y0, " (x0,y0)", :green, :left, :vcenter)
   end
end;

draw(7.5/2, 4.755/2, true)

算額（その1315）

算額（その1315）

百四十三 群馬県榛名町榛名山 榛名神社 明治33年(1900)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円15個，外円，菱形3個
#Julia, #SymPy, #算額, #和算

外円の中に菱形 3 個，等円 4 個を容れる。等円の直径が 1 寸のとき，外円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x, y)
菱形の頂点と円の接点座標を (x0, y0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, x::positive,
     y::positive, x0::positive, y0::positive;
y0 = sqrt(R^2 - x0^2)
eq1 = x^2 + y^2 - (R - r)^2
eq2 = dist2(-R, 0, x0, y0, x, y, r)
eq3 = y /R - (R - y0)/y0
eq4 = 2y0 - y - R;

eq1, eq3, eq4 を連立させて解き，x, y, x0 を求める。

res = solve([eq1, eq3, eq4], (x, y, x0))[1]

   (sqrt(-(sqrt(2)*R - r)*(-2*R + sqrt(2)*R + r)), R*(-1 + sqrt(2)), sqrt(2)*R/2)

eq2 の x, y, x0 に代入し，simplify する。

eq12 = eq2(x => res[1], y => res[2], x0 => res[3]) |> simplify
eq12 |> println

   -2*R^2 + 3*sqrt(2)*R^2/2 - R*r + sqrt(2)*R*r/2 - r^2/2 - sqrt(2)*r^2/4

R について解く。

ans_R = solve(eq12, R)[1] |> factor |> println

   r*(4 + 3*sqrt(2))/2

外円の半径 R は，等円の半径 r の (4 + 3√2)/2 倍である。
等円の直径が 1 寸のとき，外円の直径は (4 + 3√2)/2 = 4.121320343559643 寸である。

function draw(r, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = r*(4 + 3*sqrt(2))/2
   (x, y, x0) = (sqrt(-(sqrt(2)*R - r)*(-2*R + sqrt(2)*R + r)), R*(-1 + sqrt(2)), sqrt(2)*R/2)
   y0 = sqrt(R^2 - x0^2)
   plot([0, x0, -R, x0, 0, -x0, R, -x0, 0],
        [R, y0, 0, -y0, -R, -y0, 0, y0, R], color=:green, lw=0.5)
   circle(0, 0, R)
   circle4(x, y, r, :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)
       point(0, R, "R", :red, :center, :bottom, delta=delta/2)
       point(x0, y0, " (x0,y0)", :green, :left, :vcenter)
       point(x, y, "等円:r,(x,y)", :blue, :center, delta=-delta/2)
   end
end;

draw(1/2, true)

算額（その1314）

算額（その1314）

百四十 群馬県甘楽郡妙義町下高田 高太神社 明治22年(1888)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円16個，菱形
#Julia, #SymPy, #算額, #和算

菱形の中に等円を 16 個容れる。菱形の対角線の長いほうが 80 寸，菱形の一辺の長さが 50 寸のとき，等円の直径はいかほどか。

本問は，算額（その928）の類題である。

菱形の対角線を 2a, 2b;  a > b
等円の半径と中心座標を r, (x, r), (r, y), ((2r + x)/3, (r + 2y)/3), ((r + 2x)/3, (2r + y)/3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive,
     r::positive, x::positive, y::positive
b = sqrt(c^2 - a^2) 
eq1 = dist2(a, 0, 0, b, r, y, r)
eq2 = dist2(a, 0, 0, b, x, r, r)
eq3 = (x - r)^2 + (y - r)^2 - 36r^2;

eq1, eq2 から y, x を求める。

ans_y = solve(eq1, y)[1]
ans_y |> println
ans_x = solve(eq2, x)[1]
ans_x |> println

   (-c*r + (a - r)*sqrt(-a^2 + c^2))/a
   (a^2 - a*c + r*sqrt(-a^2 + c^2))/(a - c)

y, x を eq3 に代入する。

eq13 = eq3(x => ans_x, y => ans_y) |> expand |> simplify |> numerator
eq13 |> println

   a^3*c^2 - 36*a^3*r^2 - a^2*c^3 - 2*a^2*c^2*r + 36*a^2*c*r^2 + 2*a*c^3*r + 2*a*c^2*r*sqrt(-a^2 + c^2) - 2*c^3*r^2 - 2*c^2*r^2*sqrt(-a^2 + c^2)

r について解く。

ans_r = solve(eq13, r)[2]  # 2 of 2
ans_r |> println

   a*c*(6*a*(a - c) + c*(-a + c + sqrt(-a^2 + c^2)))/(2*(18*a^3 - 18*a^2*c + c^3 + c^2*sqrt(-a^2 + c^2)))

等円の直径は 2*4.54545454545454 = 9.09090909090908 = 100/11 寸である。

ans_r(a => 80/2, c=> 50, b => 30) |> println

   4.54545454545454

ans_x, ans_y は a, c の他に r を含む式であるが，上で得られた ans_r を代入すればよい。

ans_x(a => 80/2, c=> 50, r => 50/11) |> println  # 290/11
ans_y(a => 80/2, c=> 50, r => 50/11) |> println  # 230/11

   26.3636363636364
   20.9090909090909

function draw(a, c, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = sqrt(c^2 - a^2)
   r = a*c*(6*a*(a - c) + c*(-a + c + sqrt(-a^2 + c^2)))/(2*(18*a^3 - 18*a^2*c + c^3 + c^2*sqrt(-a^2 + c^2)))
   y = (-c*r + (a - r)*sqrt(-a^2 + c^2))/a
   x = (a^2 - a*c + r*sqrt(-a^2 + c^2))/(a - c)  
   @printf("a = %g;  b = %g;  c = %g\n", a, b, c)
   @printf("r = %g;  x = %g;  y = %g\n", r, x, y)
   @printf("等円の直径は %g 寸である。\n", 2r)
   plot([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:green, lw=0.5)
   circle4(r, y, r)
   circle4((2r + x)/3, (r + 2y)/3, r)
   circle4((r + 2x)/3, (2r + y)/3, r)
   circle4(x, r, r)
   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(r, y, "(r,y)", :red, :right, delta=-delta, deltax=4delta)
       point(x, r, "(x,r)", :red, :right, delta=-delta, deltax=4delta)
       point(0, b, " b", :green, :left, :bottom, delta=delta/2)
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
   end
end;

draw(80/2, 50, true)