算額（その650）
神壁算法 東都麹町 橘田彌曾八元克 天明九年
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
全円内に 2 本の弦と 2 個の等円が入っている。全円，等円の直径がそれぞれ 697寸，272寸のとき水平の弦の長さはいかほどか。
全円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x1, y1), (x2, y + r)
水平な弦と y 軸の交点座標を (0, y)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms R, r, x1, y1, x2::negative, x::negative, y::negative
eq1 = x1^2 + y1^2 - (R - r)^2
eq2 = x2^2 + (y + r)^2 - (R - r)^2
eq3 = y1/x1 * (sqrt(R^2 - x^2) - y)/(x -sqrt(R^2 - y^2)) + 1
eq4 = distance(x, sqrt(R^2- x^2), sqrt(R^2 - y^2), y, x1, y1) - r^2
eq5 = distance(x, sqrt(R^2- x^2), sqrt(R^2 - y^2), y, x2, y + r) - r^2;
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, y1, x2, x, y) = u
dy = sqrt(R^2 - x^2)
dx = sqrt(R^2 - y^2)
return [
x1^2 + y1^2 - (R - r)^2, # eq1
x2^2 - (R - r)^2 + (r + y)^2, # eq2
1 + y1*(dy - y)/(x1*(-dx + x)), # eq3
-r^2 + (x1 - (x1*(dx^2 - 2*dx*x + dy^2 - 2*dy*y + x^2 + y^2) + (dy - y)*(dx*dy - dx*y1 - dy*x1 - x*y + x*y1 + x1*y))/(dx^2 - 2*dx*x + dy^2 - 2*dy*y + x^2 + y^2))^2 + (y1 - (y1*(dx^2 - 2*dx*x + dy^2 - 2*dy*y + x^2 + y^2) + (dx - x)*(dx*dy - dx*y1 - dy*x1 - x*y + x*y1 + x1*y))/(dx^2 - 2*dx*x + dy^2 - 2*dy*y + x^2 + y^2))^2, # eq4
-r^2 + (x2 - (x2*(dx^2 - 2*dx*x + dy^2 - 2*dy*y + x^2 + y^2) + (dy - y)*(dx*dy - dx*r - dx*y - dy*x2 + r*x + x2*y))/(dx^2 - 2*dx*x + dy^2 - 2*dy*y + x^2 + y^2))^2 + (r + y - ((dx - x)*(dx*dy - dx*r - dx*y - dy*x2 + r*x + x2*y) + (r + y)*(dx^2 - 2*dx*x + dy^2 - 2*dy*y + x^2 + y^2))/(dx^2 - 2*dx*x + dy^2 - 2*dy*y + x^2 + y^2))^2, # eq5
]
end;
(R, r) = (697, 272) .// 2
iniv = BigFloat[110, 180, -200, -250, -90]
res = nls(H, ini=iniv)
(BigFloat[99.9999999999999999999999999999999999999999999999999999999999999999999970530364, 187.5000000000000000000000000000000000000000000000000000000000000000000016344705, -208.0000000000000000000000000000000000000000000000000000000000000000000319742255, -264.0000000000000000000000000000000000000000000000000000000000000000000687439946, -92.5000000000000000000000000000000000000000000000000000000000000000000529196718], true)
水平な弦の長さは 672 寸である。
その他のパラメータは以下の通り。
x1 = 100; y1 = 187.5; x2 = -208; x = -264; y = -92.5
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, r) = (697, 272) .// 2
(x1, y1, x2, x, y) = res[1]
弦 = 2sqrt(R^2 - y^2)
@printf("弦 = %g; x1 = %g; y1 = %g; x2 = %g; x = %g; y = %g\n",
弦, x1, y1, x2, x, y)
plot()
circle(0, 0, R, :magenta)
circle(x1, y1, r)
circle(x2, y + r, r)
segment(x, sqrt(R^2 - x^2), sqrt(R^2 - y^2), y, :green)
segment(-sqrt(R^2 - y^2), y, sqrt(R^2 - y^2), y, :blue)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
point(R, 0, " R", :magenta, :left, :bottom, delta=delta/2)
point(x1, y1, " 等円:r\n (x1,y1)", :red, :left, :vcenter)
point(x2, y + r, " (x2,y+r)", :red, :left, :vcenter)
point(0, y, " y", :blue, :left, delta=-delta/2)
point(x, sqrt(R^2-x^2), " (x,sqrt(R^2-x^2)", :green, :left, :bottom, delta=-delta/2)
point(sqrt(R^2-y^2), y, "(sqrt(R^2-y^2),y) ", :blue, :right, :top, delta=-delta/2)
end
end;
