東都本所亀戸村 天満宮 寛政9年(1797)
内藤忠辰校：賽祠神算
Kano Knowledge Portal
https://kanokp.lib.u-tokyo.ac.jp/item/tohoku-10010000022299
キーワード：円5個，外円，弦
#Julia #SymPy #算額 #和算 #数学
外円の中に甲円，乙円，丙円，丁円を容れる。甲円の直径が 135 寸，丙円の直径が 60 寸，丁円の直径が 54 寸のとき，外円の直径はいかほどか。
外円の半径と中心座標を \(R,\ (0,\ 0)\)
甲円の半径と中心座標を \(r_1,\ (r_1,\ y + r_1)\)
乙円の半径と中心座標を \(r_2,\ (x_2,\ y_2)\)
丙円の半径と中心座標を \(r_3,\ (x_3,\ y + r_3)\)
丁円の半径と中心座標を \(r_4,\ (x_4,\ y + r_4)\)
とおき，以下の連立方程式の解を求める。
include("julia-source.txt"); # julia-source.txt ソース
function driver(r1, r3, r4)
function H(u)
function parameters()
eq1 = x1^2 + (y + r1)^2 - (R - r1)^2
eq2 = x2^2 + y2^2 - (R - r2)^2
eq3 = x3^2 + (y + r3)^2 - (R - r3)^2
eq4 = (x1 - x2)^2 + ( (y + r1) - y2)^2 - (r1 + r2)^2
eq5 = (x1 - x4)^2 + (r1 - r4)^2 - (r1 + r4)^2
eq6 = (x2 - x3)^2 + (y2 - (y + r3))^2 - (r2 + r3)^2
eq7 = (x2 - x4)^2 + (y2 - (y + r4))^2 - (r2 + r4)^2
eq8 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2
return [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8]
end;
(R, y, x1, r2, x2, y2, x3, x4) = u
return parameters()
end;
iniv = BigFloat[65, -36, 35, 30, -24, 26, -44, -10]
iniv = BigFloat[2.8, -2.3, 1.2, 1.8, -0.84, 0.56, -1.4, -0.27].*r1
res = nls(H, ini=iniv)
res[2] || println("収束していない")
return res[1]
end;
(r1, r3, r4) = (135/2, 60/2, 54/2)
(R, y, x1, r2, x2, y2, x3, x4) = driver(r1, r3, r4)
(R, y, x1, r2, x2, y2, x3, x4) |> println
(125.0, -35.0, 47.43416490252569, 42.1875, -59.292706128157114, 57.8125, -94.86832980505137, -37.94733192202055)
# 2R: 外円の直径
2R
250.0
甲円の直径が 135 寸，丙円の直径が 60 寸，丁円の直径が 54 寸のとき，外円の直径は 250 寸である。
描画関数プログラムのソースを見る
function draw(r1, r3, r4, more=true)
pyplot(size=(500, 500), grid=false, aspectratio=1, showaxis=false, label="", fontfamily="IPAexMincho")
(R, y, x1, r2, x2, y2, x3, x4) = [65, -36, 35, 30, -24, 26, -44, -10]
(R, y, x1, r2, x2, y2, x3, x4) = driver(r1, r3, r4)
@printf("甲円の直径が %.15g，丙円の直径が %.15g，丁円の直径が %.15g のとき，外円の直径は %.15g である。\n", 2r1, 2r3, 2r4, 2R)
x = sqrt(R^2 - y^2)
plot()
circle(0, 0, R, :green)
circle(x1, y + r1, r1)
circle(x2, y2, r2, :blue)
circle(x3, y + r3, r3, :magenta)
circle(x4, y + r4, r4, :orange)
segment(-x, y, x, y)
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, 0, " 外円:R,(0,0)", :green, :left, :vcenter)
point(x1, y + r1, "甲円:r1,(x1,y_r1)", :red, :center, delta=-delta)
point(x2, y2, "乙円:r2,(x2,y2)", :blue, :center, delta=-delta)
point(x3, y + r3, "丙円:r3\n(x3,y+r3", :magenta, :center, delta=-delta)
point(x4, y + r4, "丁円:r4\n(x4,y+r4)", :orange, :center, delta=-delta)
end
end;
(r1, r3, r4) = (135/2, 60/2, 54/2)
draw(r1, r3, r4, true)
甲円の直径が 135，丙円の直径が 60，丁円の直径が 54 のとき，外円の直径は 250 である。
以下のアイコンをクリックして応援してください
