算額（その1173）

算額（その1173）

『算法用学精 巻之二』 文久2年(1862)
http://www.wasan.jp/terasima/terasima3.pdf
キーワード：円4個，外円，弦

外円の中に，水平な弦，大円 1 個，小円 2 個を容れる。外円の直径が 1 寸，大円の直径が 0.81 寸のとき，小円の直径を求めよ。

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

図形的には，算額（その750）と同等である。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, x1::positive,
     r2::positive, x2::negative;

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

SymPy では十分な簡約化ができないので，手動で簡約化する。

@syms tR, tr1
res[1](√R => tR, √r1 => tr1) |> simplify |> (x -> x(tR => √R, tr1 => √r1)) |> println

   2*sqrt(r1)*(sqrt(R) - sqrt(r1))

res[2](√R => tR, √r1 => tr1) |> simplify |> (x -> x(tR => √R, tr1 => √r1)) |> println

   2*sqrt(2)*r1^(1/4)*(R^(3/2) + 5*sqrt(R)*r1 - 4*R*sqrt(r1) - 2*r1^(3/2))/sqrt(R^(3/2) + 3*sqrt(R)*r1 - 3*R*sqrt(r1) - r1^(3/2))

res[3](√R => tR, √r1 => tr1) |> simplify |> (x -> x(tR => √R, tr1 => √r1)) |> println

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

r2 = 2*sqrt(r1)*(sqrt(R) - sqrt(r1))
t = sqrt(R^(3/2) + 3*sqrt(R)*r1 - 3*R*sqrt(r1) - r1^(3/2))
x2 = 2*sqrt(2)*r1^(1/4)*(R^(3/2) + 5*sqrt(R)*r1 - 4*R*sqrt(r1) - 2*r1^(3/2))/t
x1 = 2*sqrt(2)*r1^(1/4)*t

小円の半径 r2 は 2√r1*(√R - √r1) で求めることができる。
外円，大円の直径をそれぞれ 1 寸，0.81 寸のとき，小円の直径は 0.18 寸である。

R = 1/2
r1 = 0.81/2
r2 = 2√r1*(√R - √r1)

   0.0900000000000001

function draw(R, r1, more, plot_flag)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 2*sqrt(r1)*(sqrt(R) - sqrt(r1))
   t = sqrt(R^(3/2) + 3*sqrt(R)*r1 - 3*R*sqrt(r1) - r1^(3/2))
   u = 2*sqrt(2)*r1^(1/4)
   x1 = u*t
   x2 = u*(R^(3/2) + 5*sqrt(R)*r1 - 4*R*sqrt(r1) - 2*r1^(3/2))/t
   @printf("外円の直径が %g，大円の直径が %g のとき，小円の直径は %g である。\n", 2R, 2r1, 2r2)
   y = 2r2 - R
   x = sqrt(R^2 - y^2)
   plot()
   circle(0, 0, R, :magenta)
   circle(x1, y + r1, r1)
   circle(x2, y + r2, r2, :blue)
   circle(0, r2 - R, r2, :blue)
   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=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(R, 0, " R", :magenta, :left, :bottom, delta=delta/2)
       point(0, R, " R", :magenta, :left, :bottom, delta=delta/2)
       point(x1, y + r1, "大円:r1,(x1,2r2-R+r1)", :red, :center, delta=-delta/2)
       point(x2, y + r2, " 小円:r2,(x2,3r2-R)", :blue, :left, :vcenter)
       point(0, r2 - R, " 小円:0,(0,r2-R)", :blue, :left, :vcenter)
   end
end;

算額（その1172）

算額（その1172）

九八 鴻巣市三ツ木山王 三木神社 明治28年(1895)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円8個，外円

外円の中に，甲円 1 個，乙円 2 個，丙円 2 個，丁円 2 個を容れる。丁円の直径が 1 寸のとき，甲円の直径はいかほどか。

注：「問」には明示的に書かれておらず，『埼玉の算額』の図では下にある乙円の中心が甲円の円周上にないように見えるが，実際は「円周上にある」としないと（条件不足で）連立方程式が解けない。この制約がないと甲円の大きさに制限がつかない。甲円の大きさに制限をつけるのはおかしい気もする。甲円の大きさに制約をつけたとき，実は丙円と丁円が同じ大きさになる。「問」でも『埼玉の算額』の図でも，丙円と乙円の大きさは違うと示唆している。しかし，どの程度違うかを明示しないと甲円の大きさは決まらない。

以下では，条件不足を解消するために，丙円の大きさも指定する。

外円の半径と中心座標を R, (0, 0); R = 2r2
甲円の半径と中心座標を r1, (0, R - r1); 2r1 = 3r2
乙円の半径と中心座標を r2, (0, r2), (0, -r2)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (x4, y4); y4 = 0
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, r2::positive,
     r3::positive, x3::positive, y3::negative,
     r4::positive, x4::positive;

eq1 = x3^2 + y3^2 - (R - r3)^2
eq2 = x3^2 + (R - r1 - y3)^2 - (r1 + r3)^2
eq3 = x4^2 + (R - r1)^2 - (r1 - r4)^2
eq4 = x4^2 + r2^2 - (r2 + r4)^2
eq5 = x3^2 + (y3 + r2)^2 - (r2 + r3)^2
eq6 = R - 2r2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (R, r1, r2, x3, y3, x4))[2]

   (3*r3/2 + 3*r4/2 + sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2)/2, (3*r3 + 3*r4 + sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2))*(3*r3 + 5*r4 + sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2))/(4*(3*r3 + r4 + sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2))), 3*r3/4 + 3*r4/4 + sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2)/4, sqrt(2)*sqrt(r3)*sqrt(-r3 + 3*r4 + sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2)), 3*r3/2 - 3*r4/2 - sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2)/2, sqrt(2)*sqrt(r4)*sqrt(3*r3 + 5*r4 + sqrt(9*r3^2 - 2*r3*r4 + 9*r4^2))/2)

@syms t
println("t = sqrt(9r3^2 - 2r3*r4 + 9r4^2)")
var = ("R", "r1", "r2", "x3", "y3", "x4")
for i = 1:6
   @printf("%s = %s\n", var[i], res[i](sqrt(9r3^2 - 2r3*r4 + 9r4^2) => t))
end

   t = sqrt(9r3^2 - 2r3*r4 + 9r4^2)
   R = 3*r3/2 + 3*r4/2 + t/2
   r1 = (3*r3 + 3*r4 + t)*(3*r3 + 5*r4 + t)/(4*(3*r3 + r4 + t))
   r2 = 3*r3/4 + 3*r4/4 + t/4
   x3 = sqrt(2)*sqrt(r3)*sqrt(-r3 + 3*r4 + t)
   y3 = 3*r3/2 - 3*r4/2 - t/2
   x4 = sqrt(2)*sqrt(r4)*sqrt(3*r3 + 5*r4 + t)/2

丙円と丁円の直径が等しい 2r3 = 2r4 = 1 寸のときにのみ，算額の「答」通り，甲円の直径 2r1 = 3.75 寸になる。

r3 = r4 = 1/2
t = sqrt(9r3^2 - 2r3*r4 + 9r4^2)
r1 = (3*r3 + 3*r4 + t)*(3*r3 + 5*r4 + t)/(4*(3*r3 + r4 + t))
2r1

   3.75

function draw(r4, r3, more, plot_flag)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   t = sqrt(9r3^2 - 2r3*r4 + 9r4^2)
   R = 3*r3/2 + 3*r4/2 + t/2
   r1 = (3*r3 + 3*r4 + t)*(3*r3 + 5*r4 + t)/(4*(3*r3 + r4 + t))
   r2 = 3*r3/4 + 3*r4/4 + t/4
   x3 = sqrt(2)*sqrt(r3)*sqrt(-r3 + 3*r4 + t)
   y3 = 3*r3/2 - 3*r4/2 - t/2
   x4 = sqrt(2)*sqrt(r4)*sqrt(3*r3 + 5*r4 + t)/2
   if !plot_flag
       @printf("丁円の直径が %g のとき，甲円の直径は %g である。\n", 2r4, 2r1)
       @printf("R = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  x4 = %g\n", R, r1, r2, r3, x3, y3, r4, x4)
   end
   p1 = plot(showaxis=!plot_flag)
   circle(0, 0, R, :magenta)
   circle(0, R - r1, r1)
   circle22(0, r2, r2, :blue)
   circle2(x3, y3, r3, :orange)
   circle2(x4, 0, r4, :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(R, 0, " R", :magenta, :left, :bottom, delta=delta/2)
       point(0, R, " R", :magenta, :left, :bottom, delta=delta/2)
       !plot_flag && point(0, R - r1, @sprintf("甲円の直径 = %g", 2r1), :red, :center, :bottom, delta=delta/2)
       point(0, R - r1, "r1,(0,R-r1)", :red, :center, delta=-delta/2)
       point(0, r2, "乙円:r2,(0,r2)", :blue, :center, delta=-delta/2)
       point(0, -r2, "", :blue, :center, delta=-delta/2)
       point(x3, y3, "丙円:r3\n(x3,y3)", :orange, :center, delta=-delta/2)
       point(x4, 0, "丁円:r4\n(x4,0)", :green, :center, delta=-delta/2)
   end
   return p1
end;

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