28 岩手県一関市萩荘字八幡 八幡神社 弘化3年(1846)
安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：円5個，三角形，長方形
#Julia #SymPy #算額 #和算 #数学
長方形内に大円が内接する不等辺三角形を容れ，大きさの同じ小円を左側に複数個，右側に 1 個容れる。小円の直径がわかっているとき，左側の小円の個数と大円の直径を求めよ。
なお，前もって述べるが，術では「(左側の小円の個数÷2 + 1)×小円径=大円径」とあるが，この式は一般には成り立たない。
また，左側の小円の個数が 1 個の場合は「算額（その857）」と同じになる。が，その場合も，算額（その857）で述べた通り，示された条件だけでは不足であり，たとえば長方形の長辺の長さが既知でなければ一意な解はない。さらに，左側の小円の個数も未知として方程式を解くことは適切とは思えない。小円の個数も与えたうえで大円の直径を求めよとするのが妥当であろう。
長方形の長辺と短辺の長さを \(a,\ b\)
不等辺三角形の頂点の座標を \( (c,\ b)\)
大円の半径と中心座標を \(r_1,\ (x_1,\ r_1)\)
小円の半径と中心座標を \(r_2,\ (r_2,\ b - (2n - 1)r_2),\ (a - r_2,\ b - r_2)\)
とおき，以下の連立方程式を解く。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a::positive, b::positive, c::positive, n::integer,
r1::positive, x1::positive, r2::positive, d
l1 = sqrt(b^2 + c^2)
l2 = sqrt(b^2 + (a - c)^2)
S1 = (a + l1 + l2)*r1/2
S2 = (b + c*(2n - 1) + l1)*r2/2
S3 = ( (a - c) + b + l2)*r2/2
eq1 = (S2 + S3) - S1
eq2 = numerator(apart(dist(0, 0, c, b, x1, r1) - r1^2, d))
eq3 = numerator(apart(dist(c, b, a, 0, x1, r1) - r1^2, d))
eq4 = numerator(apart(dist(c, b, a, 0, a - r2, b - r2) - r2^2, d))
eq5 = (a - c) + b - l2 - 2r2;
function H(u)
(b, c, r1, x1) = u
return [
-r1*(a + sqrt(b^2 + c^2) + sqrt(b^2 + (a - c)^2))/2 + r2*(b + c*(2*n - 1) + sqrt(b^2 + c^2))/2 + r2*(a + b - c + sqrt(b^2 + (a - c)^2))/2, # eq1
-b^2*r1^2 + b^2*x1^2 - 2*b*c*r1*x1, # eq2
a^2*b^2 - 2*a^2*b*r1 - 2*a*b^2*x1 + 2*a*b*c*r1 + 2*a*b*r1*x1 - b^2*r1^2 + b^2*x1^2 - 2*b*c*r1*x1, # eq3
a + b - c - 2*r2 - sqrt(b^2 + (a - c)^2), # eq5
]
end;
a = 2
n = 1
r2 = 1/3
iniv = BigFloat[2, 1.5, 0.6, 1.2]
res = nls(H, ini=iniv)
([1.3333333333333333, 1.0, 0.49999999999999994, 1.0], true)
長方形の長辺が 2，小円の直径が 2/3，個数が 1 のとき，大円の直径は 1 である。
「術」の「(左側の小円の個数÷2 + 1)×小円径=大円径」は一般には成り立たない。
\(a = 4;\ b = 1.5;\ c = 2;\ n = 1;\ r_1 = 0.666667;\ x_1 = 2\)
小円の直径が 1，個数が 1，長方形の長辺が 4 のとき，大円の直径は 1.33333 である。
術による大円の直径 = 1.5
\(a = 4;\ b = 2.70711;\ c = 2.70711;\ n = 2;\ r_1 = 1;\ x_1 = 2.41421\)
小円の直径が 1，個数が 2，長方形の長辺が 4 のとき，大円の直径は 2 である。
術による大円の直径 = 2
\(a = 4;\ b = 4.10074;\ c = 2.83875;\ n = 3;\ r_1 = 1.23801;\ x_1 = 2.36272\)
小円の直径が 1，個数が 3，長方形の長辺が 4 のとき，大円の直径は 2.47602 である。
術による大円の直径 = 2.5
\(a = 4;\ b = 5.53944;\ c = 2.88985;\ n = 4;\ r_1 = 1.39379;\ x_1 = 2.29917\)
小円の直径が 1，個数が 4，長方形の長辺が 4 のとき，大円の直径は 2.78757 である。
術による大円の直径 = 3
\(a = 4;\ b = 7;\ c = 2.91667;\ n = 5;\ r_1 = 1.5;\ x_1 = 2.25\)
小円の直径が 1，個数が 5，長方形の長辺が 4 のとき，大円の直径は 3 である。
術による大円の直径 = 3.5
\(a = 20;\ b = 15.5568;\ c = 18.9657;\ n = 15;\ r_1 = 5.17517;\ x_1 = 14.4693\)
小円の直径が 1，個数が 15，長方形の長辺が 20 のとき，大円の直径は 10.3503 である。
術による大円の直径 = 8.5
描画関数プログラムのソースを見る
a = 20
n = 15
r2 = 1/2
iniv = BigFloat[2, 1.5, 0.6, 1.2].*10
res = nls(H, ini=iniv)
function draw(more=false)
pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(b, c, r1, x1) = res[1]
str = @sprintf("2r2=%g, n=%g, a=%g\n --> 2r1=%g", 2r2, n, a, 2r1)
@printf("a = %g; b = %g; c = %g; n = %g; r1 = %g; x1 = %g\n", a, b, c, n, r1, x1)
@printf("小円の直径が %g，個数が %g，長方形の長辺が %g のとき，大円の直径は %g である。\n", 2r2, n, a, 2r1)
@printf("術による大円の直径 = %g\n", (n/2 + 1)*2r2)
plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:blue, lw=0.5)
circle(x1, r1, r1)
circle(a - r2, b - r2, r2, :green)
for i = 1:n
circle(r2, b - (2i - 1)*r2, r2, :orange)
end
plot!([0, c, a], [0, b, 0], color=:magenta, 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(x1, 0.1r1, str, :black, :center, :bottom, delta=delta/2, mark=false)
point(x1, r1, "大円:r1,(x1,r1)", :red, :center, delta=-delta/2)
point(r2, b - (2n - 1)r2, "小円:r2,(r2,b-(2n-1)r2)", :black, :left, :bottom, delta=delta/2)
point(a - r2, b - r2,"小円:r2,(a-r2,b-r2)", :black, :right, delta=-delta/2)
end
end;
以下のアイコンをクリックして応援してください
