算額（その650） - 算額あれこれ

神壁算法東都麹町橘田彌曾八元克天明9年(1789)

藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
キーワード：円3個，外円，斜線
#Julia #SymPy #算額 #和算 #数学

全円内に 2 本の弦と 2 個の等円が入っている。全円，等円の直径がそれぞれ 697寸，272寸のとき水平の弦の長さはいかほどか。

全円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x1, y1), (x2, y + r)
水平な弦と y 軸の交点座標を (0, y)
とおき，以下の連立方程式を解く。

include("julia-source.txt"); # julia-source.txt ソース

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;

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)

([100.0, 187.5, -208.0, -264.0, -92.5], 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;