算額あれこれ

算額問題をコンピュータで解きます

算額(その623)

和算図形問題あれこれ

https://gunmawasan.web.fc2.com/kongetu-no-mondai.html
キーワード:円5個,正方形,斜線
#Julia #SymPy #算額 #和算 #数学


正方形内に斜線を引き,5 個の等円を並べる。等円の直径が 1 寸のとき,正方形の辺の長さ,斜線の長さはいかほどか。

正方形の辺の長さを \(a\), 斜線と正方形の辺の交点座標を \( (b, 0), (0, c)\) とおく。
等円の半径と中心座標を \(r, (r, a - r), (x_1, y_1), (x_2, y_2), (a - 3r, r), (a - r, r)\)
とおき,以下の連立方程式を解く。

なお,中心座標が \( (r, a - r), (a - 3r, r)\) の 2 個の等円の中心間の距離についての方程式 \( (a - 4r)^2 + (a - 2r)^2 = 36r^2\) を解けば等円の半径が \(r\) のとき,正方形の辺の長さは \(a = r(3 + \sqrt{17})\) であることは容易にわかる。つまり,\(r = 1/2\)  なら \(a = 3.5615528128088303\) である。
ただし,このようにして \(a\) を求めても \(b, c\) は簡単には求まらないし,図を描くときのために \( (x_1, y_1), (x_2, y_2)\) も求めるとすれば,以下の 7 元連立方程式を解くほうが簡単である。

r = 1/2
r*(3 + sqrt(17))

   3.5615528128088303

include("julia-source.txt");  # julia-source.txt ソース

using SymPy

@syms a::positive, b::positive, c::positive, r::positive,
     x1::positive, y1::positive, x2::positive, y2::positive
r = 1//2
eq1 = (x1 - r)^2 + (a - r - y1)^2 - 4r^2
eq2 = (x2 - x1)^2 + (y1 - y2)^2 - 4r^2
eq3 = (a - 3r - x2)^2 + (y2 - r)^2 - 4r^2
eq4 = distance(b, 0, 0, c, r, a - r) - r^2
eq5 = distance(b, 0, 0, c, x1, y1) - r^2
eq6 = distance(b, 0, 0, c, x2, y2) - r^2
eq7 = distance(b, 0, 0, c, a - 3r, r) - r^2;

# res = solve([eq2, eq3, eq4, eq5, eq6, eq7], (a, b, c, x1, y1, x2, y2));

function H(u)
   (a, b, c, x1, y1, x2, y2) = u
   return [
       (x1 - 1/2)^2 + (a - y1 - 1/2)^2 - 1,  # eq1
       (-x1 + x2)^2 + (y1 - y2)^2 - 1,  # eq2
       (y2 - 1/2)^2 + (a - x2 - 3/2)^2 - 1,  # eq3
       (-b*(-2*a*c + b + 2*c^2 + c)/(2*(b^2 + c^2)) + 1/2)^2 + (a - c*(2*a*c + 2*b^2 - b - c)/(2*(b^2 + c^2)) - 1/2)^2 - 1/4,  # eq4
       (-b*(b*x1 + c^2 - c*y1)/(b^2 + c^2) + x1)^2 + (-c*(b^2 - b*x1 + c*y1)/(b^2 + c^2) + y1)^2 - 1/4,  # eq5
       (-b*(b*x2 + c^2 - c*y2)/(b^2 + c^2) + x2)^2 + (-c*(b^2 - b*x2 + c*y2)/(b^2 + c^2) + y2)^2 - 1/4,  # eq6
       (-c*(-2*a*b + 2*b^2 + 3*b + c)/(2*(b^2 + c^2)) + 1/2)^2 + (a - b*(2*a*b - 3*b + 2*c^2 - c)/(2*(b^2 + c^2)) - 3/2)^2 - 1/4,  # eq7
   ]
end;
r = 1/2
iniv = BigFloat[3.54, 1.71, 2.79, 0.96, 2.14, 1.54, 1.29]
res = nls(H, ini=iniv)

   ([3.5615528128088303, 1.7807764064044151, 2.9211646096066226, 1.02051760426961, 2.207701875205887, 1.5410352085392203, 1.3538509376029435], true)

正方形の一辺の長さは 3.56155 寸,斜線の長さは = sqrt(b^2 + c^2) = 3.42116 寸である。

その他のパラメータは以下の通り。
\(a = 3.56155;  b = 1.78078;  c = 2.92116;  x_1 = 1.02052;  y_1 = 2.2077\)
\(x_2 = 1.54104;  y_2 = 1.35385\)

描画関数プログラムのソースを見る

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   (a, b, c, x1, y1, x2, y2) = res[1]
   斜線 = sqrt(b^2 + c^2)
   @printf("正方形の一辺の長さ = %g;  斜線の長さ = %g;  a = %g;  b = %g;  c = %g;  x1 = %g;  y1 = %g;  x2 = %g;  y2 = %g\n",
       a, 斜線, a, b, c, x1, y1, x2, y2)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5, xlims=(-0.2, 3.7), ylims=(-0.2, 3.7))
   circle(r, a - r, r, :blue)
   circle(x1, y1, r, :blue)
   circle(x2, y2, r, :blue)
   circle(a - 3r, r, r, :blue)
   circle(a - r, r, r, :blue)
   segment(0, c, b, 0, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(r, a - r, "(r,a-r)")
       point(x1, y1, "(x1,y1)")
       point(x2, y2, "(x2,y2)")
       point(a - 3r, r, "(a-3r,r)")
       point(a - r, r, "(a-r,r)")
       point(a, 0, " a", :black, :left, :top, delta=-delta/2)
       point(b, 0, " b", :green, :left, :top, delta=-delta/2)
       point(0, a, "a ", :black, :right, :vcenter)
       point(0, c, "c ", :green, :right, :vcenter)
   end
end;


以下のアイコンをクリックして応援してください