13 岩手県奥州市（旧水沢市佐倉河） 胆沢城跡鎮守府八幡宮 弘化2年(1845)
安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：円6個，外円，楕円，弦
#Julia #SymPy #算額 #和算 #数学
全円の中に，弦 1 本，楕円 1 個，大円 4 個，小円 1 個を容れる。小円の直径が与えられたとき，全円の直径を求める術を述べよ。
一般性を失うことなく，楕円の長径，短径を与えて小円と全円の直径を求め，小円と全円の関係式を求めることにしてよい。
図を時計回りに 90° 回転させて解く。
楕円の長半径，短半径，中心座標を \(a, b, (0, 0)\)
全円の半径と中心座標を \(R, (a - R, 0)\)
大円の半径と中心座標を \(r_1, (x_1, y_1), (-r_1, 0)\)
小円の半径と中心座標を \(r_2, (-a - r2, 0)\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a::positive, b::positive,
r1::positive, x1::negative, y1::positive,
r2::positive, R::positive,
x0::negarive, y0::positive
# eq0 = 4(a^2 - b^2)*(b^2 - r1^2)/b^2 - (2r1)^2 # 算法助術の公式84
r1 = b*sqrt(a^2 - b^2)/a
x1 = r1 - a
eq1 = r2 - (R - a)
eq2 = x0^2/a^2 + y0^2/b^2 - 1
eq3 = (x0 - x1)^2 + (y0 - y1)^2 - r1^2
eq4 = -b^2*x0/(a^2*y0) + (x0 - x1)/(y0 - y1)
eq5 = (x1 - (a - R))^2 + y1^2 - (R - r1)^2;
SymPy では能力的に解けないので，数値解を求める。
function driver(a, b)
function H(u)
(R, r2, y1, x0, y0) = u
return [
-R + a + r2, # eq1
-1 + y0^2/b^2 + x0^2/a^2, # eq2
(y0 - y1)^2 + (a + x0 - b*sqrt(a^2 - b^2)/a)^2 - b^2*(a^2 - b^2)/a^2, # eq3
(a + x0 - b*sqrt(a^2 - b^2)/a)/(y0 - y1) - b^2*x0/(a^2*y0), # eq4
y1^2 - (R - b*sqrt(a^2 - b^2)/a)^2 + (R - 2*a + b*sqrt(a^2 - b^2)/a)^2, # eq5
]
end;
iniv = BigFloat[6, 0.96, 3.46146, -2.65805,1.69398]
return nls(H, ini=iniv)
end;
res = driver(5, 2)
@printf("R = %g, r2 = %g, R/r2 = %g\n", res[1][1], res[1][2], res[1][1]/res[1][2])
R = 5.9428, r2 = 0.942805, R/r2 = 6.30333
res = driver(7, 3)
@printf("R = %g, r2 = %g, R/r2 = %g\n", res[1][1], res[1][2], res[1][1]/res[1][2])
R = 8.57247, r2 = 1.57247, R/r2 = 5.45159
術には「小円の直径を 5 倍すれば全円の直径が得られる」とあるが，
長半径，短半径が 5, 2 のとき，\(R/r_2\) ≒ 6.30333
長半径，短半径が 7, 3 のとき，\(R/r_2\) ≒ 5.45159
のように，\(R/r_2\) は一定ではない。
描画関数プログラムのソースを見る
function draw(a, b, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
res = driver(a, b)
(R, r2, y1, x0, y0) = res[1]
r1 = b*sqrt(a^2 - b^2)/a
x1 = r1 - a
@printf("x1= %g\n", x1)
r2 = R - a
@printf("R = %g; r2 = %g; R/r2 = %g\n", R, r2, R/r2)
@printf("a = %g; b = %g; r2 = %g; r1 = %g; x1 = %g; y1 = %g; R = %g\n", a, b, r2, r1, x1, y1, R)
p = plot()
circle(a - R, 0, R, :green)
ellipse(0, 0, a, b, color=:blue)
circle2(r1, 0, r1)
circle22(x1, y1, r1)
circle(-a - r2, 0, r2, :magenta)
y = sqrt(R^2 - (R - 2a)^2)
segment(-a, y, -a, -y, :black)
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 - R, 0, "全円:R,(a-R,0)", :green, :right, :bottom, delta=delta, deltax=5delta)
point(x1, y1, "大円:r1,(x1,y1)", :red, :center, delta=-delta)
point(-r1, 0, "大円:r1,(-r1,0)", :red, :center, delta=-delta)
point(-a - r2, 0, "小円:r2,(-a-r2,0)", :black, :center, delta=-delta, deltax=4delta)
point(a, 0, " a", :blue, :left, :vcenter)
point(0, b, "b", :blue, :center, :bottom, delta=delta)
else
point(-a, y1, @sprintf(" R/r2 = %g", R/r2), :black, :left, :vcenter, mark=false)
end
return p
end;
draw(5, 2, true)
以下のアイコンをクリックして応援してください
