算額（その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)