十三 群馬県前橋市元総社町 総社神社 享和年中(1801〜1804)
三十三 群馬県佐波郡境町伊与久 雷電神社 文化13年(1816)
三十九 群馬県前橋市元総社町 総社神社 文政6年(1823)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円6個，三角形
#Julia #SymPy #算額 #和算 #数学
三角形の中に甲円 1 個，乙円 1 個，丙円 2 個，丁円 2 個を容れる。甲円の直径が 8 寸のとき，乙円，丙円，丁円の直径はいかほどか。
三角形の底辺の長さを a，頂点の座標を \( (b, c)\)
甲円の半径と中心座標を \(r_1, (x_1, r_1)\)
乙円の半径と中心座標を \(r_2, (x_2, r_2)\)
丙円の半径と中心座標を \(r_3, (x_{31}, r_3), (x_{32}, r_3)\)
丁円の半径と中心座標を \(r_4, (x_{41}, r_4), (x_{42}, y_{42})\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
function driver(r1)
function H(u)
function parameters()
eq1 = x2 - x1 - 2sqrt(r1*r2)
eq2 = x32 - x2 - 2sqrt(r2*r3)
eq3 = x41 - x32 - 2sqrt(r3*r4)
eq4 = x1 - x31 - 2sqrt(r1*r3)
eq5 = (x1 - x42)^2 + (y42 - r1)^2 - (r1 + r4)^2
eq6 = r4/(a - x41) - r1/(a - x1)
eq7 = r3/(a - x32) - r1/(a - x1)
eq8 = r2/(a - x2) - r1/(a - x1)
eq9 = r3/x31 - r1/x1
eq10 = dist3(b, c, a, 0, x2, r2, r2)
eq11 = dist3(b, c, a, 0, x32, r3, r3)
eq12 = dist3(b, c, a, 0, x42, y42, r4)
eq13 = dist3(b, c, 0, 0, x31, r3, r3)
return [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13]
end;
(a, b, c, x1, r2, x2, r3, x31, x32, r4, x41, x42, y42) = u
return parameters()
end;
iniv = BigFloat[20.818, 5.28282, 9.60552, 7.7588, 2.27159, 12.81312, 1.27967, 2.1886, 16.22124, 0.74556, 18.20508, 5.53907, 8.57916]
res = nls(H, ini=iniv)
res[2] || println("収束していない")
return res[1]
end;
driver(8/2) |> println
[20.770561253610573, 5.261860482353991, 9.62855983402497, 6.744263886800197, 2.278228309790651, 12.781784317958207, 1.2975810578828915, 2.1878072672188953, 16.220496809759116, 0.739046475079226, 18.179039846772014, 5.542828364137464, 8.5842244903413]
甲円の直径 = 8 のとき，乙円直径 = 4.56083，丙円直径 = 2.60015，丁円直径 = 1.48236 である。
描画関数プログラムのソースを見る
function draw(r1, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, b, c, x1, r2, x2, r3, x31, x32, r4, x41, x42, y42) = (20.818, 5.28282, 9.60552, 6.7588, 2.27159, 12.81312, 1.27967, 2.1886, 16.22124, 0.74556, 18.20508, 5.53907, 8.57916)
(a, b, c, x1, r2, x2, r3, x31, x32, r4, x41, x42, y42) = res[1]
@printf("甲円の直径 = %g，乙円直径 = %g，丙円直径 = %g，丁円直径 = %g\n", 2r1, 2r2, 2r3, 2r4)
plot([0, a, b, 0], [0, 0, c, 0], color=:blue, lw=0.5)
circle(x1, r1, r1)
circle(x2, r2, r2, :green)
circle(x31, r3, r3, :magenta)
circle(x32, r3, r3, :magenta)
circle(x41, r4, r4, :brown)
circle(x42, y42, r4, :brown)
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(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=-delta)
point(x2, r2, "乙円:r2\n(x2,r2)", :green, :center, delta=-delta)
point(x31, r3, " 丙円:r3,(x31,r2)", :magenta, :center, delta=-10delta, deltax=8delta)
point(x32, r3, " 丙円:r3,(x32,r2)", :magenta, :left, :bottom, delta=12delta, deltax=-5delta)
point(x41, r4, " 丁円:r4,(x41,r4)", :brown, :center, delta=-6delta, deltax=-5delta)
point(x42, y42, " 丁円:r4,(x42,y42)", :brown, :left, :bottom, delta=2delta, deltax=5delta)
point(b, c, "(b,c)", :blue, :center, :bottom, delta=delta)
point(a, 0, "a", :blue, :left, :bottom, delta=delta)
ylims!(-1, 10)
end
end;
draw(8/2, true)
以下のアイコンをクリックして応援してください
