(12) 京都市伏見区御香宮門前町 御香宮神社（ごこうのみやじんじゃ）文久3年(1863)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円6個，直線上
#Julia #SymPy #算額 #和算 #数学
直線上に互いに隣同士外接する甲，乙，丙，丁，戊の 5 個の円が載っており，それらの全てに外接する己円が載っている。甲，乙，丙円の直径がそれぞれ 4 寸，1.2 寸，1 寸のとき，戊円の直径はいかほどか。
甲円の半径と中心座標を \(r_1,\ (x_1,\ r_1)\)
乙円の半径と中心座標を \(r_2,\ (x_2,\ r_2)\)
丙円の半径と中心座標を \(r_3,\ (x_3,\ r_3)\)
丁円の半径と中心座標を \(r_4,\ (x_4,\ r_4)\)
戊円の半径と中心座標を \(r_5,\ (x_5,\ y_6)\)
己円の半径と中心座標を \(r_6,\ (0,\ y_6)\)
とおき，以下の連立方程式を解く。
なお，上の図は甲，乙，丙円の直径がそれぞれ 4 寸，1.29 寸，1 寸のときものである。乙円の直径が違うだけで，図は大きく変わる。与えられた条件のときの図は，後で掲示する。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms r1::positive, x1::positive, r2::positive, x2::positive,
r3::positive, x3::positive, r4::positive, x4::positive,
r5::positive, r6::positive, x6::positive, y6::positive
#@syms r1, x1, r2, x2, r3, x3, r4, x4, r6, x6, y6, r5
#(r1, r2, r3) = (4//2, 12//20, 1//2)
eq1 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2 |> expand
eq2 = (x1 - x6)^2 + (r1 - y6)^2 - (r1 + r6)^2 |> expand
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2 |> expand
eq4 = (x2 - x6)^2 + (r2 - y6)^2 - (r2 + r6)^2 |> expand
eq5 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2 |> expand
eq6 = (x3 - x6)^2 + (r3 - y6)^2 - (r3 + r6)^2 |> expand
eq7 = (x4 - x6)^2 + (r4 - y6)^2 - (r4 + r6)^2 |> expand
eq8 = x4^2 + (r4 - r5)^2 - (r4 + r5)^2 |> expand
eq9 = x6^2 + (y6 - r5)^2 - (r6 + r5)^2 |> expand;
function H(u)
(x1, x2, x3, r4, x4, r5, r6, x6, y6) = u
return [
-4*r1*r2 + x1^2 - 2*x1*x2 + x2^2, # eq1
-2*r1*r6 - 2*r1*y6 - r6^2 + x1^2 - 2*x1*x6 + x6^2 + y6^2, # eq2
-4*r2*r3 + x2^2 - 2*x2*x3 + x3^2, # eq3
-2*r2*r6 - 2*r2*y6 - r6^2 + x2^2 - 2*x2*x6 + x6^2 + y6^2, # eq4
-4*r3*r4 + x3^2 - 2*x3*x4 + x4^2, # eq5
-2*r3*r6 - 2*r3*y6 - r6^2 + x3^2 - 2*x3*x6 + x6^2 + y6^2, # eq6
-2*r4*r6 - 2*r4*y6 - r6^2 + x4^2 - 2*x4*x6 + x6^2 + y6^2, # eq7
-4*r4*r5 + x4^2, # eq8
-2*r5*r6 - 2*r5*y6 - r6^2 + x6^2 + y6^2, # eq9
]
end;
(r1, r2, r3) = (4, 1.2, 1)./2
#(r1, r2, r3) = (4.0, 1.29, 1.0)./2
iniv = BigFloat[11.5, 9.3, 8.2, 0.94, 6.8, 12.5, 1, 8.5, 2]
res = nls(H, ini=iniv)
([11.502173707608492, 9.311283477587828, 8.215838362577495, 0.9375000000000001, 6.84653196881458, 12.500000000000012, 1.0125, 8.489699641330079, 1.9874999999999998], true)
甲，乙，丙円の直径がそれぞれ 4 寸，1.2 寸，1 寸のとき，戊円の半径は 12.500000000000012 である（直径は 25 寸）。
このときの図は以下のようになる。
「術」でも，戊円の直径は 25 寸としており，間違いではない。
パラメータは以下のとおりである。
\(r_1 = 2;\ x_1 = 11.5022\)
\( r_2 = 0.6;\ x_2 = 9.31128\)
\( r_3 = 0.5;\ x_3 = 8.21584\)
\( r_4 = 0.9375;\ x_4 = 6.84653\)
\( r_5 = 12.5\)
\( r_6 = 1.0125;\ x_6 = 8.4897;\ y_6 = 1.9875\)
描画関数プログラムのソースを見る
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(x1, x2, x3, r4, x4, r5, r6, x6, y6) = res[1]
@printf("r1 = %g; x1 = %g\n", r1, x1)
@printf("r2 = %g; x2 = %g\n", r2, x2)
@printf("r3 = %g; x3 = %g\n", r3, x3)
@printf("r4 = %g; x4 = %g\n", r4, x4)
@printf("r5 = %g\n", r5)
@printf("r6 = %g; x6 = %g; y6 = %g\n", r6, x6, y6)
plot()
circle(x1, r1, r1)
circle(x2, r2, r2, :blue)
circle(x3, r3, r3, :magenta)
circle(x4, r4, r4, :purple)
circle(0, r5, r5, :green)
circle(x6, y6, r6, :orange)
segment(-r5, 0, x1+r1, 0, :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(x1, r1, "甲:r1,(x1,r1)", :red, :center, delta=-delta/2)
point(x2, r2, "乙:r2,(x2,r2)", :blue, :left, delta=-delta/2)
point(x3, r3, "丙:r3,(x3,r3)", :magenta, :center, delta=-6delta)
point(x4, r4, "丁:r4,(x4,r4)", :purple, :right, delta=-delta/2)
point(0, r5, "戊:r5,(0,r5)", :green, :center, delta=-delta/2)
point(x6, y6, "己:r6,(x6,y6)", :orange, :center, delta=-delta/2)
end
end;
以下のアイコンをクリックして応援してください
