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

68 岩手県川崎村薄衣字諏訪前（現一関市川崎町）浪分神社文政5年(1822)

安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：円6個，外円，斜線
#Julia #SymPy #算額 #和算 #数学

外円の中に弦を 4 本引き，区画された領域に等円 5 個を容れる。等円の直径が与えられたとき，外円の直径を求める術を述べよ。

外円の半径と中心座標を \(R, (0, 0)\)
等円の半径と中心座標を \(r, (0, R - r), (x_1, y_1), (x_2, y_2)\)
弦の端点座標を \( (x_{01}, \sqrt{R^2 - x_{01}^2}), (x_{02}, \sqrt{R^2 - x_{02}^2}), (0, -R)\)
とおき，以下の連立方程式を解く。

\(R, r_1\) を記号のままにして SymPy では簡単に解くことができない。数値を代入して解けば，数値解は求まる。

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

using SymPy
@syms R::positive, r::positive, x1::positive, y1::positive,
x2::positive, y2::negative, x01::positive, x02::positive
eq1 = x1^2 + y1^2 - (R - r)^2
eq2 = x2^2 + y2^2 - (R - r)^2
eq3 = dist(0, -R, x01, sqrt(R^2 - x01^2), 0, R - r) - r^2
eq4 = dist(0, -R, x01, sqrt(R^2 - x01^2), x1, y1) - r^2
eq5 = dist(0, -R, x02, sqrt(R^2 - x02^2), x1, y1) - r^2
eq6 = dist(0, -R, x02, sqrt(R^2 - x02^2), x2, y2) - r^2
eq7 = (sqrt(R^2 - x02^2) + R)/x02 * y2/x2 + 1;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (R, x1, y1, x2, y2, x01, x02))

function H(u)
(R, x1, y1, x2, y2, x01, x02) = u
return [
x1^2 + y1^2 - (R - r)^2,
x2^2 + y2^2 - (R - r)^2,
dist(0, -R, x01, sqrt(R^2 - x01^2), 0, R - r) - r^2,
dist(0, -R, x01, sqrt(R^2 - x01^2), x1, y1) - r^2,
dist(0, -R, x02, sqrt(R^2 - x02^2), x1, y1) - r^2,
dist(0, -R, x02, sqrt(R^2 - x02^2), x2, y2) - r^2,
(sqrt(R^2 - x02^2) + R)/x02 * y2/x2 + 1
]
end;
r = 2
iniv = BigFloat[10, 5, 5, 7, -3, 3, 8]
res = nls(H, ini=iniv)

([7.509090269827044, 3.8101684797908897, 3.9790315098916964, 4.870544338179145, -2.5744656631880116, 2.2798713379711115, 6.204719439604394], true)

等円の直径が 4 のとき，外円の直径は 2*7.509090269827044 = 15.018180539654088 である。

術は，「外円の直径を x，等円の直径を d として以下の d の 5 次式を解く」としている。
\(-d^5 + 4d^4 x - 14d^3 x^2 + 52d^2 x^3 - 73 d x^4 + 16x^5\)

using SymPy
@syms x, d
eq = ( ( ( (16x - 73d)*x + 52d^2)*x - 14d^3)*x + 4d^4)*x - d^5 |> expand
eq |> println

-d^5 + 4*d^4*x - 14*d^3*x^2 + 52*d^2*x^3 - 73*d*x^4 + 16*x^5

d = 4 のときの解を求める。5 個の解のうち，3 番目のものが適解 15.0181805396541 である。

res = solve(eq(d => 4));
res[3].evalf() |> println

15.0181805396541

描画関数プログラムのソースを見る

function draw(r, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, x1, y1, x2, y2, x01, x02) = res[1]
plot()
circle(0, 0, R, :green)
circle(0, R - r, r)
circle2(x1, y1, r)
circle2(x2, y2, r)
plot!([-x01, 0, x01], [sqrt(R^2 - x01^2), -R, sqrt(R^2 - x01^2)], color=:blue, lw=0.5)
plot!([-x02, 0, x02], [sqrt(R^2 - x02^2), -R, sqrt(R^2 - x02^2)], color=:blue, lw=0.5)
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(0, R, "R", :green, :center, :bottom, delta=delta)
point(x1, y1, "r,(x1,y1)", :red, :center, delta=-delta)
point(x2, y2, "r,(x2,y2)", :red, :center, delta=-delta)
point(x01, sqrt(R^2 - x01^2), "(x01,y01)", :blue, :left, :bottom, delta=delta)
point(x02, sqrt(R^2 - x02^2), "(x02,y02)", :blue, :left, :bottom, delta=delta)
point(0, R - r, " R-r", :blue, :left, :vcenter)
end
end;