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

算額（その1313）

算額（その1313）

百三十八 群馬県利根郡月夜野町上津 八幡神社 明治22年(1889)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円39個，円環
#Julia, #SymPy, #算額, #和算

大円の内外にそれぞれ 39 個の等円（外等円，内等円と呼ぶ）が互いに接している。
外等円と内等円の直径の差が 8.454 寸，大円と内等円の直径の和が 171.305 寸，「外余積が 109.688535160125 歩」のとき，大円，外等円，内等円の直径はいかほどか。

算額（その672）と同じであるが，問題の提示法に難がある。

大円の半径と中心座標を R, (0, 0)
外等円の半径と中心座標を r1, ((R + r1)*cos(pi/19), (R + r1)*sin(pi/19)); (R + r1)*sin(pi/19) = r1
内等円の半径と中心座標を r2, ((R - r2)*cos(pi/19), (R - r2)*sin(pi/19)); (R - r2)*sin(pi/19) = r2
とする。
未知数は R, r1, r2 の 3 個なので，方程式も 3 本用意すれば解ける。
外等円と内等円は互いに接し合っているという 2 条件は eq1, eq2 で表されている。
残り1つの条件が必要であるが，「問」には 3 個の条件が提示されている。3 番目の条件に出てくる「外余積」は有効数字が多いので，この条件を使えば正確な解が得られると思われるが，この「外余積」が何であるかわからない。
残る 2 条件のいずれを使うかで解が異なる。それは「外等円と内等円の直径の差が 8.454 寸」も「大円と内等円の直径の和が 171.305 寸」も有理数ではなく誤差を含むのが原因である。
更に，後で分かるが，eq1, eq2 に使われている「円周率」は 3.16 である（円積率=0.79）。

1.　「外等円と内等円の直径の差が 8.454 寸」を使う場合

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, r2::positive,
     円周率::positive;
円周率 =3.16
eq1 = (R + r1)*sin(円周率/19) - r1
eq2 = (R - r2)*sin(円周率/19) - r2
eq3 = 2(r1 - r2) - 8.454;

res = solve([eq1, eq2, eq3], (r1, r2, R))
res |> println

   Dict{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}(r2 => 10.6530275137650, R => 75.0022912276876, r1 => 14.8800275137650)

大円の直径 = 2*75.0022912276876 = 150.0045824553752
外等円の直径 = 2*14.8800275137650 = 29.76005502753
内等円の直径 = 2*10.6530275137650 = 21.30605502753

2.　「大円と内等円の直径の和が 171.305 寸」を使う場合

@syms R::positive, r1::positive, r2::positive,
     円周率::positive;
円周率 =3.16
eq1 = (R + r1)*sin(円周率/19) - r1
eq2 = (R - r2)*sin(円周率/19) - r2
eq3 = 2(R + r2) - 171.305;

res = solve([eq1, eq2, eq3], (r1, r2, R))
res |> println

   Dict{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}(r2 => 10.6526769444053, R => 74.9998230555947, r1 => 14.8795378424756)

大円の直径 = 2*74.9998230555947 = 149.9996461111894
外等円の直径 = 2*14.8795378424756 = 29.7590756849512
内等円の直径 = 2*10.6526769444053 = 21.3053538888106

3.　大円の直径を与える場合

r1, r2 は 大円 R の関数なので，条件は eq1, eq2 だけで十分である。求められた r1, r2 は R を含む式なので，R として任意の値を設定すれば，その場合の外等円，内等円の半径（直径）はすぐに得られる。

大円の直径 2R は，前 2 項から 150 寸と推定される（実際に「答」で 150 寸とされている）。

@syms R::positive, r1::positive, r2::positive,
     円周率::positive;
円周率 = Sym(316)/100
eq1 = (R + r1)*sin(円周率/19) - r1
eq2 = (R - r2)*sin(円周率/19) - r2
ans_r1 = solve(eq1, r1)[1]
ans_r2 = solve(eq2, r2)[1]
(ans_r1, ans_r2) |> println
(2ans_r1(R => 150/2).evalf(), 2ans_r2(R => 150/2).evalf())

   (R*sin(79/475)/(1 - sin(79/475)), R*sin(79/475)/(sin(79/475) + 1))

   (29.7591458944761, 21.3054041537713)

外等円，内等円の半径はそれぞれ
R*sin(79/475)/(1 - sin(79/475))
R*sin(79/475)/(sin(79/475) + 1)
であり，R = 150/2 を代入し 2 倍すると，それぞれの直径は 29.7591458944761, 21.3054041537713 となり，小数点以下 3 桁までとると，29.759, 21.305 となり，「答」と一致する。

4.　eq1, eq2 を，円周率 = π で計算する場合

近似計算ではなく正確な解を求める。

@syms R::positive, r1::positive, r2::positive,
     円周率::positive;
円周率 = PI
eq1 = (R + r1)*sin(円周率/19) - r1
eq2 = (R - r2)*sin(円周率/19) - r2
ans_r1 = solve(eq1, r1)[1]
ans_r2 = solve(eq2, r2)[1]
(ans_r1, ans_r2) |> println
(2ans_r1(R => 150/2).evalf(), 2ans_r2(R => 150/2).evalf())

   (R*sin(pi/19)/(1 - sin(pi/19)), R*sin(pi/19)/(sin(pi/19) + 1))

   (29.5535416156890, 21.1998138649765)

大円の直径が 150 寸のとき，外等円の直径は 29.5535 寸, 内等円の直径は 21.1998 寸である。

function draw(R, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (R*sin(pi/19)/(1 - sin(pi/19)), R*sin(pi/19)/(sin(pi/19) + 1))
   @printf("大円の直径が %g のとき，外等円の直径は %g, 内等円の直径は %g である。\n", 2R, 2r1, 2r2)
   @printf("R = %g;  r1 = %g;  r2 = %g\n", R, r1, r2)
   plot()
   circle(0, 0, R, :magenta)
   circle(0, 0, R + r1, :gray90)
   circle(0, 0, R - r2, :gray90)
   rotate(0, -r1 - R, r1, :blue, angle=360/19) 
   rotate(0, r2 - R, r2, angle=360/19) 
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
   end
end;

draw(150/2, true)

算額（その1312）

算額（その1312）

百三十七 群馬県藤岡市藤岡 金光寺 明治21年(1888)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，円弧，最大化
#Julia, #SymPy, #算額, #和算

団扇の中に円弧，甲円，乙円を容れる。団扇の直径，円弧の直径がそれぞれ 4 寸，3 寸のとき，乙円の直径が最大になるのは甲円の直径がいかほどのときか。

団扇の直径と中心座標を R, (0, 0)
円弧の直径と中心座標を r0, (0, y0); y0 < 0
甲円の半径と中心座標を r1, (0, r1 - R)
乙円の半径と中心座標を r2, (x2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, r0::positive, y0::positive,
     r1::positive, r2::positive, x2::positive, y2::positive;
y0 = 2r1 - R - r0
eq1 = x2^2 + (y2 - y0)^2 - (r0 - r2)^2
eq2 = x2^2 + y2^2 - (R - r2)^2
eq3 = x2^2 + (y2 - r1 + R)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3], (r2, x2, y2))[2]  # 2 of 2

   (-r1*(-R + r1)*(r0 - r1)/(R*r0 - r1^2), 2*sqrt(R)*sqrt(r0)*r1*sqrt(-(-R + r1)*(r0 - r1))/(R*r0 - r1^2), (-R^2*r0 + R*r0*r1 + r0*r1^2 - r1^3)/(R*r0 - r1^2))

ans_r2 = res[1];

乙円の半径 r2 は -r1*(-R + r1)*(r0 - r1)/(R*r0 - r1^2) である。
R, r0 が与えられるので，変化しうるのは甲円の半径 r1 である。
r2 と r1 の関係を下図に示す。r1 が 0.5〜1 の間で r2 が最大値をとることがわかる。

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="")
plot(ans_r2(R => 4/2, r0 => 3/2), xlims=(0,1.5), xlabel="r1", ylabel="r2")
vline!([1.37879462408182/2])
hline!([0])

r2 が最大になるときの r1 は，r2 の式を r1 で偏微分し，その式が 0 になるときの r1 を求めればよい。

r10 = solve(diff(ans_r2, r1), r1)[2];  # 2 of 4
r10 |> println
2r10(R => 4/2, r0 => 3/2) |> println

   -sqrt(2*R^2*r0^2/(9*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) + 8*R*r0/3 + 2*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3))/2 + sqrt(-2*R^2*r0^2/(9*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) + 16*R*r0/3 + (4*R^2*r0 + 4*R*r0^2)/sqrt(2*R^2*r0^2/(9*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) + 8*R*r0/3 + 2*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) - 2*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3))/2

   1.37879462408182

乙円の直径が最大になるときの甲円の直径を表す式はとても長く，SymPy では簡約化できないが，数値は求めることができる。

団扇の直径が 4 寸, 円弧の直径が 3 寸で，甲円の直径が 1.37879 寸のときに乙円の直径は最大値 0.580181 寸になる。

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

    R = 2;  r0 = 1.5;  y0 = -2.12121;  r1 = 0.689397;  r2 = 0.290091;  x2 = 0.974955;  y2 = -1.40473

function draw(R, r0, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = -sqrt(2*R^2*r0^2/(9*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) + 8*R*r0/3 + 2*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3))/2 + sqrt(-2*R^2*r0^2/(9*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) + 16*R*r0/3 + (4*R^2*r0 + 4*R*r0^2)/sqrt(2*R^2*r0^2/(9*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) + 8*R*r0/3 + 2*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3)) - 2*(-26*R^3*r0^3/27 + (2*R^2*r0 + 2*R*r0^2)^2/16 + sqrt(-R^6*r0^6/729 + (52*R^3*r0^3/27 - (2*R^2*r0 + 2*R*r0^2)^2/8)^2/4))^(1/3))/2
   (r2, x2, y2) = (-r1*(-R + r1)*(r0 - r1)/(R*r0 - r1^2), 2*sqrt(R)*sqrt(r0)*r1*sqrt(-(-R + r1)*(r0 - r1))/(R*r0 - r1^2), (-R^2*r0 + R*r0*r1 + r0*r1^2 - r1^3)/(R*r0 - r1^2))
   y0 = 2r1 - R - r0
   @printf("団扇の直径が %g, 円弧の直径が %g で，甲円の直径が %g のときに乙円の直径は最大値 %g になる。\n", 2R, 2r0, 2r1, 2r2)
   @printf("R = %g;  r0 = %g;  y0 = %g;  r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", R, r0, y0, r1, r2, x2, y2)
   plot()
   circle(0, 0, R)
   circle(0, y0, r0, :blue)
   circle(0, r1 - R, r1, :green)
   circle2(x2, y2, r2, :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, 0, "団扇:R,(0,0)", :red, :center, delta=-delta/2)
       point(0, y0, "円弧:r0,(0,y0)", :blue, :center, delta=-delta/2)
       point(0, r1 - R, "甲円:r1\n(0,r1-R)", :green, :center, delta=-delta/2)
       point(x2, y2, "乙円:r2,(x2,y2)", :magenta, :left, :bottom, delta=4delta, deltax=-4delta)
   end
end;

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

算額（その1311）

算額（その1311）

百三十七 群馬県藤岡市藤岡 金光寺 明治21年(1888)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，楕円，最大化
#Julia, #SymPy, #算額, #和算

大円と楕円が交わり，其の隙間に等円 3 個が入っている。楕円の長径が与えられたとき，等円の直径が最大になるのは楕円の短径がいかほどのときか。

楕円の長半径，短半径，中心座標を a, b, (0, 0)
大円の半径と中心座標を r1, (0, r1 - b)
等円の半径と中心座標を r2, (x2, 0), (0, b + r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r1::positive,
     r2::positive, x2::positive;
r1 = b + r2
eq1 = x2^2 + (r1 -b)^2 - (r1 + r2)^2
eq2 = x2^2 - (a^2 - b^2)*(b^2 - r2^2)/b^2
res = solve([eq1, eq2], (r2, x2))[1]

   ((8*a^2*b^2*(a - b)*(a + b)/(a^2 + 2*b^2) - 4*a^2*b^2 + 4*b^4)/(4*b*(a - b)*(a + b)), 2*sqrt(2)*a*b*sqrt(a - b)*sqrt(a + b)/(a^2 + 2*b^2))

x2 は簡約化できる。

ans_x2 = res[2] |> simplify
ans_x2 |> println

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

等円の半径も簡約化され，a, b を含む関数である。a は与えられるので，変化できるのは b である。 0 < b < √2a/2

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

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

a = 1 のとき，楕円の短半径 b と小円の半径 r2 は以下の図のようになる。

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="")
plot(ans_r2(a => 1), xlims=(0,√2/2), xlabel="b", ylabel="r2")
vline!([0.34356074972251255])
hline!([0])

r2 が最大になるときの b は，r2 の式を b で偏微分し，その式が 0 になるときの b を求めればよい。

b0 = solve(diff(ans_r2, b), b)[1]
b0 |> println
b0(a => 1).evalf() |> println

   a*sqrt(-1 + sqrt(5)/2)
   0.343560749722512

楕円の長半径 a のとき，等円の半径が最大になるのは，楕円の短半径が b = a*sqrt(√5/2 - 1) のときである。

a = 1 のとき，b = sqrt(√5/2 - 1) = 0.343560749722512 である。

そのときの等円の半径は最大値 sqrt(√5 - 2)*(√10 - √2)/4 = 0.212332220528909 になる。

r2_max = ans_r2(b => b0) |> simplify |> factor
r2_max |> println
r2_max(a => 1).evalf() |> println

   a*sqrt(-2 + sqrt(5))*(-sqrt(2) + sqrt(10))/4
   0.212332220528909

function draw(a, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = a*sqrt(-1 + sqrt(5)/2)
   r2 = a*sqrt(√5 - 2)*(√10 - √2)/4
   r1 = b + r2
   x2 = 2a*b*sqrt(2a^2 - 2b^2)/(a^2 + 2b^2)
   plot()
   ellipse(0, 0, a, b, color=:red)
   circle(0, r1 - b, r1, :green)
   circle2(x2, 0, r2, :blue)
   circle(0, b + 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(a, 0, "a", :red, :left, :bottom, delta=delta/2, deltax=delta)
       point(0, b, "b", :red, :center, :bottom, delta=2delta)
       point(0, r1 - b, "大円:r1,(0,r1-b)", :green, :center, delta=-2delta)
       point(0, b + r2, "等円:r2\n(0,b+r2)", :blue, :center, :bottom, delta=2delta)
       point(x2, 0, "等円:r2\n(x2,0)", :blue, :center, :bottom, delta=2delta)
   end
end;

draw(30, true)