90 岩手県陸前高田市小友町上の坊 常膳寺観音堂 天保13年(1842)
安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：円7個，正三角形，斜線
#Julia #SymPy #算額 #和算 #数学
正三角形の中に斜線を 2 本，全円 1 個，丙円 1 個，乙円 3 個，全円 1 個を容れる。丙円の直径が与えられたとき，甲円の直径を得る術を述べよ。
正三角形の一辺の長さを \(2a\)
甲円の半径と中心座標を \(r_1, (x_3, y_3)\)
丙円の半径と中心座標を \(r_3, (0, 2r_0 - r_3)\)
乙円の半径と中心座標を \(r_2, (x_2, y_2), (0, 2r_0 + r_2)\)
全円の半径と中心座標を \(r_0, (0, r_0)\)
斜線と正三角形の一辺の交点座標を \( (b, \sqrt{3}(a - b))\)
とおき，以下の連立方程式を解く。
なお，\(r_0 = \sqrt{3}a/3, r_2 = \sqrt{3}a/9, y_3 = y_2\) である。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a, r0, r1, r2, x2, y2, r3, x3, y3, b
# eq1 = r0 - a/√Sym(3)
# eq2 = 2r2 - (√Sym(3)a - 2r0 - r2)
# res = solve([eq1, eq2], (r0, r2))
# r0 = √Sym(3)a/3
# r2 = √Sym(3)a/9
r0 = √Sym(3)a/3
r2 = √Sym(3)a/9
y3 = y2
# dist2() は使えない
eq3 = dist(-a, 0, b, √Sym(3)*(a - b), x2, y2) - r2^2
eq4 = dist(-a, 0, b, √Sym(3)*(a - b), 0, 2r0 - r3) - r3^2
eq5 = dist(a, 0, -b, √Sym(3)*(a - b), x2, y2) - r2^2
eq6 = dist(a, 0, 0, √Sym(3)a, x2, y2) - r2^2
eq7 = dist(-a, 0, b, √Sym(3)*(a - b), x3, y3) - r1^2
eq8 = (x2 - x3) - (r2 + r1);
println(eq3, ", # eq3")
println(eq4, ", # eq4")
println(eq5, ", # eq5")
println(eq6, ", # eq6")
println(eq7, ", # eq7")
println(eq8, ", # eq8")
-a^2/27 + (y2 - sqrt(3)*(a - b)*(sqrt(3)*y2*(a - b) + (a + b)*(a + x2))/(3*(a - b)^2 + (a + b)^2))^2 + (a + x2 - (a + b)*(sqrt(3)*y2*(a - b) + (a + b)*(a + x2))/(3*(a - b)^2 + (a + b)^2))^2, # eq3
-r3^2 + (a - (a + b)*(a*(a + b) + sqrt(3)*(a - b)*(2*sqrt(3)*a/3 - r3))/(3*(a - b)^2 + (a + b)^2))^2 + (2*sqrt(3)*a/3 - r3 - sqrt(3)*(a - b)*(a*(a + b) + sqrt(3)*(a - b)*(2*sqrt(3)*a/3 - r3))/(3*(a - b)^2 + (a + b)^2))^2, # eq4
-a^2/27 + (y2 - sqrt(3)*(a - b)*(sqrt(3)*y2*(a - b) + (-a - b)*(-a + x2))/( (-a - b)^2 + 3*(a - b)^2))^2 + (-a + x2 - (-a - b)*(sqrt(3)*y2*(a - b) + (-a - b)*(-a + x2))/( (-a - b)^2 + 3*(a - b)^2))^2, # eq5
-a^2/27 + (y2 - sqrt(3)*a*(sqrt(3)*a*y2 - a*(-a + x2))/(4*a^2))^2 + (-a + x2 + a*(sqrt(3)*a*y2 - a*(-a + x2))/(4*a^2))^2, # eq6
-r1^2 + (y2 - sqrt(3)*(a - b)*(sqrt(3)*y2*(a - b) + (a + b)*(a + x3))/(3*(a - b)^2 + (a + b)^2))^2 + (a + x3 - (a + b)*(sqrt(3)*y2*(a - b) + (a + b)*(a + x3))/(3*(a - b)^2 + (a + b)^2))^2, # eq7
-sqrt(3)*a/9 - r1 + x2 - x3, # eq8
function H(u)
(x2, y2, r1, b, r3, x3) = u
return [
-a^2/27 + (y2 - sqrt(3)*(a - b)*(sqrt(3)*y2*(a - b) + (a + b)*(a + x2))/(3*(a - b)^2 + (a + b)^2))^2 + (a + x2 - (a + b)*(sqrt(3)*y2*(a - b) + (a + b)*(a + x2))/(3*(a - b)^2 + (a + b)^2))^2, # eq3
-r3^2 + (a - (a + b)*(a*(a + b) + sqrt(3)*(a - b)*(2*sqrt(3)*a/3 - r3))/(3*(a - b)^2 + (a + b)^2))^2 + (2*sqrt(3)*a/3 - r3 - sqrt(3)*(a - b)*(a*(a + b) + sqrt(3)*(a - b)*(2*sqrt(3)*a/3 - r3))/(3*(a - b)^2 + (a + b)^2))^2, # eq4
-a^2/27 + (y2 - sqrt(3)*(a - b)*(sqrt(3)*y2*(a - b) + (-a - b)*(-a + x2))/( (-a - b)^2 + 3*(a - b)^2))^2 + (-a + x2 - (-a - b)*(sqrt(3)*y2*(a - b) + (-a - b)*(-a + x2))/( (-a - b)^2 + 3*(a - b)^2))^2, # eq5
-a^2/27 + (y2 - sqrt(3)*a*(sqrt(3)*a*y2 - a*(-a + x2))/(4*a^2))^2 + (-a + x2 + a*(sqrt(3)*a*y2 - a*(-a + x2))/(4*a^2))^2, # eq6
-r1^2 + (y2 - sqrt(3)*(a - b)*(sqrt(3)*y2*(a - b) + (a + b)*(a + x3))/(3*(a - b)^2 + (a + b)^2))^2 + (a + x3 - (a + b)*(sqrt(3)*y2*(a - b) + (a + b)*(a + x3))/(3*(a - b)^2 + (a + b)^2))^2, # eq7
-sqrt(3)*a/9 - r1 + x2 - x3, # eq8
]
end;
a = 10
r0 = √3a/3 # a*0.5773502691896257
r2 = √3a/9 # a*0.19245008972987523
iniv = BigFloat[0.575411, 0.34912, 0.095137, 0.662547, 0.387466, 0.284131] .* a
res = nls(H, ini=iniv)
([5.711594854986474, 3.5787338000056996, 0.9544492311458473, 6.575250125060548, 3.8641388667743293, 2.832644726541875], true)
正三角形の一辺の長さが \(2a = 20\) のとき，
全円の半径 \(r_0 = 5.773502691896257\)
甲円の半径 \(r_1 = 0.9544492311458473\)
乙円の半径 \(r_2 = 1.9245008972987523\)
丙円の半径 \(r_3 = 3.8641388667743293\)
丙円の直径/甲円の直径 = 3.8641388667743293/0.9544492311458473 = 4.048553595810754
丙円の直径は，甲円の直径の 4.048553595810754 倍である。
術は，「15972 の平方根から 121 を引き，甲円の直径を掛ければ丙円の直径がえられる」とあるが，\(\sqrt{15972} - 121 = 5.38037822383663\) であり，題意を満たさない。
描画関数プログラムのソースを見る
function draw(a, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(x2, y2, r1, b, r3, x3) = [5.711594854986474, 3.5787338000056996, 0.9544492311458473, 6.575250125060548, 3.8641388667743293, 2.832644726541875]
y3 = y2
r0 = √3a/3
r2 = √3a/9
plot([a, 0, -a, a], [0, sqrt(3)a, 0, 0], color=:green, lw=0.5)
circle(0, a/sqrt(3), r0)
circle(0, 2r0 + r2, r2, :blue)
circle2(x2, y2, r2, :blue)
circle2(x3, y3, r1, :brown)
circle(0, 2r0 - r3, r3, :orange)
segment(-a, 0, b, (a - b)*√3)
segment(a, 0, -b, (a - b)*√3)
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, √3a, " √3a", :green, :left, :vcenter)
point(a, 0, "a", :green, :left, :bottom, delta=delta/2)
point(0, r0, "全円:r0,(0,r0)", :red, :center, delta=-delta)
point(0, 2r0 - r3, "丙円:r3,(0,2r0-r3)", :orange, :center, delta=-delta)
point(x2, y2, "乙円:r2\n(x2,y2)", :blue, :left, delta=-delta, deltax=-3delta)
point(0, 2r0 + r2, " 乙円:r2\n(x2,2r0+r2)", :blue, :center, :bottom, delta=delta)
point(x3, y3, "甲円:r1,(x3,y3)", :brown, :right, delta=-4delta, deltax=-delta)
point(b, √3(a - b), " (b,√3(a-b))", :green, :left, :vcenter)
end
end;
draw(10, true)
以下のアイコンをクリックして応援してください
