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