算額（その695）

算額（その695）

神壁算法 關流藤田貞資門人 東都四谷右京殿町 中村幸藏永著 寛政八年
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

外円内に三斜（不等辺三角形），と 9 個の円が入っている。乙円，丁円，己円の直径がそれぞれ，3 寸，2 寸，1 寸のとき，外円の直径はいかほどか。

円周上にある三角形の頂点座標を (x0, y0), (x, y); y0 = 2r1 - R, x0 = sqrt(R^2 - y0^2), x = -sqrt(R^2 - y^2)
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, r1 - R)
乙円の半径と中心座標を r2, (x2, y0 - r2);
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (x4, y4)
戊円の半径と中心座標を r5, (x5, y5)
己円の半径と中心座標を r6, (x6, y6)
とおき，以下の連立方程式を解く。

eq7, eq13 は「算法助術」の公式29 である。

中村信弥編著: 和算の図形公式--『算法助術』--- 
　　http://www.wasan.jp/kosiki/kosiki.html

include("julia-source.txt");

using SymPy
@syms r2, r4, r6, R, x, y, r1, x2, r3, x3, y3, x4, y4, r5, x5, y5, x6, y6

(r2, r4, r6) = (3, 2, 1) .// 2
y0 = 2r1 - R
x0 = sqrt(R^2 - y0^2)
x = -sqrt(R^2 - y^2)
len1 = sqrt((x0 - x)^2 + (y - y0)^2)
len2 = sqrt((x0 + x)^2 + (y - y0)^2)
len3 = sqrt(R^2 - (2r1 - R)^2)

eq1 = x2^2 + (y0- r2)^2 - (R - r2)^2
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2

eq3 = x3^2 + y3^2 - (R - r3)^2
eq4 = x4^2 + y4^2 - (R - r4)^2
eq5 = (x3 - x4)^2 + (y3 - y4)^2 - (r3 + r4)^2
eq6 = (R - 2r3)^2 + (len1/2)^2 - R^2
eq7 = len1^2 - 4*2R*2r4  # 公式29
eq8 = (2r1 - R - y)/(len3 - x) * y3/x3 + 1

eq9 = x5^2 + y5^2 - (R - r5)^2
eq10 = x6^2 + y6^2 - (R - r6)^2
eq11 = (x5 - x6)^2 + (y5 - y6)^2 - (r5 + r6)^2
eq12 = (R - 2r5)^2 + (len2/2)^2 - R^2
eq13 = len2^2 - 4*2R*2r6  # 公式29
eq14 = (2r1 - R - y)/(len3 + x) * y5/x5 - 1;

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)
   (R, y, r1, x2, r3, x3, y3, x4, y4, r5, x5, y5, x6, y6) = u
   return [
       x2^2 - (R - 3/2)^2 + (-R + 2*r1 - 3/2)^2,  # eq1
       x2^2 + (r1 - 3/2)^2 - (r1 + 3/2)^2,  # eq2
       x3^2 + y3^2 - (R - r3)^2,  # eq3
       x4^2 + y4^2 - (R - 1)^2,  # eq4
       -(r3 + 1)^2 + (x3 - x4)^2 + (y3 - y4)^2,  # eq5
       -R^2 + (R - 2*r3)^2 + (sqrt(R^2 - y^2) + sqrt(R^2 - (-R + 2*r1)^2))^2/4 + (R - 2*r1 + y)^2/4,  # eq6
       -16*R + (sqrt(R^2 - y^2) + sqrt(R^2 - (-R + 2*r1)^2))^2 + (R - 2*r1 + y)^2,  # eq7
       1 + y3*(-R + 2*r1 - y)/(x3*(sqrt(R^2 - y^2) + sqrt(R^2 - (-R + 2*r1)^2))),  # eq8
       x5^2 + y5^2 - (R - r5)^2,  # eq9
       x6^2 + y6^2 - (R - 1/2)^2,  # eq10
       -(r5 + 1/2)^2 + (x5 - x6)^2 + (y5 - y6)^2,  # eq11
       -R^2 + (R - 2*r5)^2 + (-sqrt(R^2 - y^2) + sqrt(R^2 - (-R + 2*r1)^2))^2/4 + (R - 2*r1 + y)^2/4,  # eq12
       -8*R + (-sqrt(R^2 - y^2) + sqrt(R^2 - (-R + 2*r1)^2))^2 + (R - 2*r1 + y)^2,  # eq13
       -1 + y5*(-R + 2*r1 - y)/(x5*(-sqrt(R^2 - y^2) + sqrt(R^2 - (-R + 2*r1)^2))),  # eq14
   ]
end;

(r2, r4, r6) = (3, 2, 1) .// 2
iniv = BigFloat[12, 5.5, 3, 4.25, 1.35, 2.75, 3.75, 4.5, 2.2, 0.5, -4.6, 3, -4, 3.8]
res = nls(H, ini=iniv)

   (BigFloat[5.999999999999999999999999999999999999999999999999999999999999999999999999999931, 5.656854249492380195206754896838792314278687501507792292706718951962929913848408, 2.999999999999999999999999999999999999999999999999999999999999999999999999999965, 4.242640687119285146405066172629094235709015626130844219530039213972197435386254, 1.267949192431122706472553658494127633057194746189619371944193020548066983091185, 2.732050807568877293527446341505872366942805253810380628055806979451933016908746, 3.863703305156273146998972798915589470535619356033618201609372305241692855919384, 4.416153642713558008419684820119938656145252405562964502449789354746723815760582, 2.344693370986443557691308802220881143166169766978954662289509621552255928898189, 0.5505102572168219018027159252941086080340525193433298715673074327490396225426752, -4.449489742783178098197284074705891391965947480656670128432692567250960377457299, 3.14626436994197234232913506571557044551247712918732870123248671744266549537083, -3.802437397408490550876329002854725019354381540512873852459501536515540189749589, 3.973848240533267382224213011697627892737753679640435536379640797303909768676809], true)

外円の直径は 12 寸である。

その他のパラメータは以下のとおり。
R = 6;  x = -2;  y = 5.65685;  r1 = 3;  x2 = 4.24264;  r3 = 1.26795;  x3 = 2.73205;  y3 = 3.8637;  x4 = 4.41615;  y4 = 2.34469;  r5 = 0.55051;  x5 = -4.44949;  y5 = 3.14626;  x6 = -3.80244;  y6 = 3.97385

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r4, r6) = (3, 2, 1) .// 2
   (R, x, y, r1, x2, r3, x3, y3, x4, y4, r5, x5, y5, x6, y6) =  [6, -2.4, 5.5, 3, 4.25, 1.35, 2.75, 3.75, 4.5, 2.2, 0.5, -4.6, 3, -4, 3.8]
   (R, y, r1, x2, r3, x3, y3, x4, y4, r5, x5, y5, x6, y6) = res[1]
   x = -sqrt(R^2 - y^2)
   @printf("R = %g;  x = %g;  y = %g;  r1 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  y3 = %g;  x4 = %g;  y4 = %g;  r5 = %g;  x5 = %g;  y5 = %g;  x6 = %g;  y6 = %g\n", R, x, y, r1, x2, r3, x3, y3, x4, y4, r5, x5, y5, x6, y6)
   y0 = 2r1 - R
   x0 = sqrt(R^2 - y0^2)
   plot()
   circle(0, 0, R, :blue)
   segment(-x0, y0, x0, y0)
   segment(x0, y0, x, y)
   segment(-x0, y0, x, y)
   circle(0, r1 - R, r1)
   circle(x2, y0 - r2, r2, :green)
   circle(-x2, y0 - r2, r2, :green)
   circle(x3, y3, r3, :magenta)
   circle(x4, y4, r4, :orange)
   θ1 = 2atand(y3/x3) - atand(y4/x4)
   circle((R - r4)*cosd(θ1), (R - r4)*sind(θ1), r4, :orange)
   circle(x5, y5, r5, :tomato)
   circle(x6, y6, r6, :brown)
   θ2 = 2atand(y5/x5) - atand(y6/x6)
   circle((r6 - R)*cosd(θ2), (r6 - R)*sind(θ2), r6, :brown)
   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)", :black, :right, :bottom, delta=delta/2)
       point(0, R, " R", :blue, :left, :bottom, delta=delta/2)
       point(x0, y0, "(x0,y0)   ", :black, :right, :bottom, delta=delta/2)
       point(0, r1 - R, " 甲円:r1\n (0,r1-R)", :red, :left, :vcenter)
       point(x2, y0 - r2, " 乙円:r2\n (x2,y0-r2)", :green, :center, delta=-delta/2)
       point(x3, y3, "丙円:r3\n(x3,y3)", :magenta, :center, delta=-delta/2)
       point(x4, y4, "丁円:r4\n(x4,y4)", :black, :center, delta=-delta/2)
       point(x5, y5, " 戊円:r5,(x5,y5)", :tomato, :left, :vcenter)
       point(x6, y6, " 己円:r6,(x6,y6)", :brown, :left, :vcenter)
   end
end;

算額（その694）

算額（その694）

神壁算法 關流藤田貞資門人 久世大和守家士 平井彌五太夫正義 寛政八年
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

鈎股弦（直角三角形）内に，斜線，全円，甲円，乙円が入っている。全円，甲円の直径がそれぞれ 19寸7分5厘，3寸1分6厘 のとき，乙円の直径はいかほどか。

直角を挟む2辺（鈎，股）の長さを「鈎」，「股」
斜線と斜辺の交点座標を (x, y)
全円の半径と中心座標を r0, (r0, r0)
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (r2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r0, r1, 鈎, 股, x1, r2, y2, x, y;

eq1 = 鈎 + 股 - sqrt(鈎^2 + 股^2) - 2r0
eq2 = (r0 - x1)^2 + (r0 - r1)^2 - (r0 + r1)^2
eq3 = (r0 - r2)^2 + (r0 - y2)^2 - (r0 + r2)^2
eq4 = dist(0, 0, x, y, x1, r1) - r1^2
eq5 = dist(0, 0, x, y, r2, y2) - r2^2
eq6 = (鈎 + 股 + sqrt(鈎^2 + 股^2))*r0 - 鈎*股
eq7 = 鈎/股 - y/(股 - x);

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)
   (鈎, 股, x1, r2, y2, x, y) = u
   return [
       -2*r0 + 股 + 鈎 - sqrt(股^2 + 鈎^2),  # eq1
       (r0 - r1)^2 - (r0 + r1)^2 + (r0 - x1)^2,  # eq2
       (r0 - r2)^2 - (r0 + r2)^2 + (r0 - y2)^2,  # eq3
       -r1^2 + (r1 - y*(r1*y + x*x1)/(x^2 + y^2))^2 + (-x*(r1*y + x*x1)/(x^2 + y^2) + x1)^2,  # eq4
       -r2^2 + (r2 - x*(r2*x + y*y2)/(x^2 + y^2))^2 + (-y*(r2*x + y*y2)/(x^2 + y^2) + y2)^2,  # eq5
       r0*(股 + 鈎 + sqrt(股^2 + 鈎^2)) - 股*鈎,  # eq6
       -y/(-x + 股) + 鈎/股,  # eq7
   ]
end;

(r0, r1) = (1975, 316) .// 200
iniv = BigFloat[26, 49, 2, 0.6, 5, 5, 24]
res = nls(H, ini=iniv)

   (BigFloat[26.50958690810937781981564610070982992662591012710315095504511369567094160353881, 48.6025397559463070232560931975192400228754697146791710191805408688433795666828, 1.974999999999999999999999999999999999999999999999999999999999999999999999999969, 0.5700002001441297307780475616969606044325776049176168848094355451524229609489859, 5.130001801297167577002428055272645439893198444258551963284919906371806648542264, 5.312669489352098474921664382153576018135850078829535519263800791306165244941906, 23.61186439712043766631850836512700452504822257257571341895022573913851219974303], true)

乙円の直径 = 1.1400004002882596

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

鈎 = 26.5096;  股 = 48.6025;  x1 = 1.975;  r2 = 0.57;  y2 = 5.13;  x = 5.31267;  y = 23.6119

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1) = (1975, 316) .// 200
   (鈎, 股, x1, r2, y2, x, y) = res[1]
   println("乙円の直径 = $(Float64(2r2))")
   @printf("鈎 = %g;  股 = %g;  x1 = %g;  r2 = %g;  y2 = %g;  x = %g;  y = %g\n", 鈎, 股, x1, r2, y2, x, y)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(r0, r0, r0, :blue)
   circle(x1, r1, r1)
   circle(r2, y2, r2, :green)
   segment(0, 0, x, y, :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(x, y, " (x,y)", :black, :left, :bottom)
       point(r0, r0, " 全円:r0\n (r0,r0)", :blue, :left, :vcenter)
       point(x1, r1, "  甲円:r1,(x1,r1)", :red, :left, :bottom, delta=delta/2)
       point(r2, y2, " 乙円:r2,(r2,y2)", :green, :left, :vcenter)
       point(股, 0, " 股", :black, :left, :bottom, delta=delta/2)
       point(0, 鈎, " 鈎", :black, :left, :bottom, delta=delta/2)
   end
end;

 

算額（その693）

算額（その693）

神壁算法 關流藤田貞資門人 毛利石見守家士 山邊平助清誠 寛政七年
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

算額（その602） に似ている（求めるものが違うが，本質的に同じ）。

外円内に長方形を n 個入れる。それぞれの長方形の長辺の 2 頂点は外円の円周上にあり，残りの 2 頂点は隣の長方形と共有している。
長方形の長辺が 3 寸のとき，外積が最小（n 個の面積の和が最大）になるのは，外円の直径ががいかほどのときか。

外円の半径を R, 内接する長方形の長辺の長さを 2a とする。
長方形の頂点座標 (x, a), (sqrt(R^2 - a^2), a) を求める。
長方形の面積の和は n*(2a * (sqrt(R^2 - a^2) - x)) である。
外積は外円の面積から長方形の面積の和を差し引いたものである。これを B とすれば，B = πR^2 - n*(2a * (sqrt(R^2 - a^2) - x)) となる。

なお，上の図は「神壁算法」にある図に近いように描いたものであるが，これは外円の直径が 13 寸のときのもので，外積は 81.4913891769524 なので，最小値からはかけ離れている。

実際に外積が最小となる場合の図は下の方に示す。

include("julia-source.txt");

using SymPy
@syms n, a, R, B;

x = a/tand(Sym(180)/n)
B = PI*R^2 - n*2a*(sqrt(R^2 - a^2) - x)
B |> println

   pi*R^2 - 2*a*n*(-a/tan(pi/n) + sqrt(R^2 - a^2))

たとえば，問のように n = 10, a = 3/2 の場合には，以下の B2 のようになる。

B2 = B(n => 10, a => 3//2)
B2 |> println

   pi*R^2 - 30*sqrt(R^2 - 9/4) + 45/sqrt(1 - 2*sqrt(5)/5)

B2 の値は R により変化するが，R = 5 近辺で最小値を取ることがわかる。

using Plots
pyplot(size=(300, 150), grid=false, aspectratio=:none, label="")
plot(B2, xlims=(0,10), xlabel="外円の半径", ylabel="黒積")

B が最小値となるのは，B の導関数が 0 になるときである。

まず，導関数 g を求める。

g = diff(B, R)
g |> println

   -2*R*a*n/sqrt(R^2 - a^2) + 2*pi*R

導関数 = 0 となる R を求めるために solve() を用いる。

3 組の解が得られるが，3 番目のものが適解である。

R0 = solve(g, R)
R0 |> println
R0[3] |> println

   Sym[0, -a*sqrt(n^2 + pi^2)/pi, a*sqrt(n^2 + pi^2)/pi]
   a*sqrt(n^2 + pi^2)/pi

a, n に具体的な値を代入すれば，黒積が最小になるときの外円の半径が得られる。

R0[3](a=>3/2, n => 10).evalf() |> println

   5.00472439995710

外円の半径が与えられると黒積が得られる。

B(R => R0[3]) |> println

   a^2*(n^2 + pi^2)/pi - 2*a*n*(-a/tan(pi/n) + sqrt(a^2*(n^2 + pi^2)/pi^2 - a^2))

a, n に具体的な値を代入すれば，黒積が最小のときの値が得られる。

B(R => R0[3], a=>3/2, n => 10).evalf() |> println

   73.9446182521105

長辺の長さが 3 の 10 個の長方形を入れるとき，外円の直径が 10.0094 のときに黒積が最小値 73.9446 をとる。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (n, a) = (10, 3//2)  # 条件として与えられる
   x = a/tand(Sym(180)/n)
   R = a*sqrt(n^2 + pi^2)/pi
   B = a^2*(n^2 + pi^2)/pi - 2*a*n*(-a/tan(pi/n) + sqrt(a^2*(n^2 + pi^2)/pi^2 - a^2))
   B0 = pi*R^2 - n*2a*(sqrt(R^2 - a^2) - x)
   @printf("長辺の長さが %g の %i 個の長方形を入れるとき，外円の直径が %g のときに黒積が最小値 %g をとる\n", 2a, n, 2R, B)
   plot()
   circle(0, 0, R, :blue)
   XY = [x sqrt(R^2 - a^2) sqrt(R^2 - a^2) x x;
           -a -a a a -a]
   for deg in 0:360/n:359
       XY2 = [cosd(deg) -sind(deg); sind(deg) cosd(deg)] * XY
       plot!(XY2[1, :], XY2[2, :])
   end
   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", :blue, :left, :bottom, delta=delta/2)
       point(x, a, "(x,a) ", :black, :right, :top)
       point(√(R^2-a^2), a, "(√(R^2-a^2),a) ", :black, :right, :bottom, delta=delta)
   end
end;