算額（その797） - 算額あれこれ

藤田貞資：精要算法（下巻）天明元年(1781)

http://www.wasan.jp/seiyou/seiyou.html
キーワード：円3個，正方形，斜線
#Julia #SymPy #算額 #和算 #数学

正方形の中に2本の斜線を引き，できた区画の中に大円，中円，小円を容れる。大円，中円，小円の直径および2本の斜線の長さと正方形の一辺の長さの和が 55 寸のとき，正方形の一辺の長さはいかほどか。

大円の半径と中心座標を \(r_1,\ (r_1,\ a - r_1)\)
中円の半径と中心座標を \(r_2,\ (a - r_2,\ y_2)\)
小円の半径と中心座標を \(r_3,\ (r_3,\ r_3)\)
とおいて以下の連立方程式を解く。

SymPy の能力的に，解析解を求めることができないので，数値解を求める。

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

using SymPy

@syms r1::positive, r2::positive, y2::positive,
r3::positive, a::positive, b::positive,
c::positive, d::positive
eq1 = dist(b, 0, a, a, r1, a - r1) - r1^2
eq2 = dist(b, 0, a, a, a - r2, y2) - r2^2
eq3 = dist(b, 0, a, a, r3, r3) - r3^2
eq4 = dist(c,0, 0, d, r1, a - r1) - r1^2
eq5 = dist(c,0, 0, d, a - r2, y2) - r2^2
eq7 = c + d - sqrt(c^2 + d^2) - 2r3
eq8 = a + (a - b) - sqrt(a^2 + (a - b)^2) - 2r2
eq9 = 2r1 + 2r2 + 2r3 + sqrt(a^2 + (a - b)^2) + sqrt(c^2 + d^2) + a - 55;

function H(u)
(r1, r2, y2, r3, a, b, c, d) = u
return [
-r1^2 + (a - a*(a*(a - r1) + (a - b)*(-b + r1))/(a^2 + (a - b)^2) - r1)^2 + (-b + r1 - (a - b)*(a*(a - r1) + (a - b)*(-b + r1))/(a^2 + (a - b)^2))^2, # eq1
-r2^2 + (-a*(a*y2 + (a - b)*(a - b - r2))/(a^2 + (a - b)^2) + y2)^2 + (a - b - r2 - (a - b)*(a*y2 + (a - b)*(a - b - r2))/(a^2 + (a - b)^2))^2, # eq2
-r3^2 + (-a*(a*r3 + (a - b)*(-b + r3))/(a^2 + (a - b)^2) + r3)^2 + (-b + r3 - (a - b)*(a*r3 + (a - b)*(-b + r3))/(a^2 + (a - b)^2))^2, # eq3
-r1^2 + (a - d*(-c*(-c + r1) + d*(a - r1))/(c^2 + d^2) - r1)^2 + (-c + c*(-c*(-c + r1) + d*(a - r1))/(c^2 + d^2) + r1)^2, # eq4
-r2^2 + (-d*(-c*(a - c - r2) + d*y2)/(c^2 + d^2) + y2)^2 + (a - c + c*(-c*(a - c - r2) + d*y2)/(c^2 + d^2) - r2)^2, # eq5
c + d - 2*r3 - sqrt(c^2 + d^2), # eq7
2*a - b - 2*r2 - sqrt(a^2 + (a - b)^2), # eq8
a + 2*r1 + 2*r2 + 2*r3 + sqrt(a^2 + (a - b)^2) + sqrt(c^2 + d^2) - 55, # eq9
]
end;

iniv = BigFloat[44, 30, 30, 21, 124, 35, 85, 62]./10
res = nls(H, ini=iniv)

([4.0, 3.0, 3.0, 2.0, 12.0, 3.0, 8.0, 6.0], true)

正方形の一辺の長さは 12 寸である。

その他のパラメータは以下のとおりである。

\(r_1 = 4;\ r_2 = 3;\ y_2 = 3;\ r_3 = 2;\ a = 12;\ b = 3;\ c = 8;\ d = 6\)

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

function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, y2, r3, a, b, c, d) = res[1]
@printf("正方形の一辺の長さは %g 寸である。\n", a)
@printf("r1 = %g; r2 = %g; y2 = %g; r3 = %g; a = %g; b = %g; c = %g; d = %g\n", r1, r2, y2, r3, a, b, c, d)
plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
circle(r1, a - r1, r1)
circle(a - r2, y2, r2, :blue)
circle(r3, r3, r3, :green)
segment(c, 0, 0, d, :magenta)
segment(b, 0, a, a, :magenta)
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(r1, a - r1, "大円:r1,(r1,a-r1)", :red, :center, delta=-delta)
point(a - r2, y2, "中円:r2,(a-r2,y2)", :blue, :center, delta=-delta)
point(r3, r3, "小円:r3,(r3,r3)", :green, :center, delta=-delta)
point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
point(b, 0, " b", :magenta, :left, :bottom, delta=delta/2)
point(c, 0, "c ", :magenta, :right, :bottom, delta=delta/2)
point(0, d, "d ", :magenta, :right, :vcenter)
plot!(xlims=(-0.5, 12.5))
end
end;