長野県長野市戸隠 戸隠山中院権現堂 安永5年(1776)
中村信弥「改訂増補 長野県の算額」県内の算額(P.26)
http://www.wasan.jp/zoho/zoho.html
キーワード：円4個，外円
#Julia #SymPy #算額 #和算 #数学
外円の中に中円 1 個，小円 2 個が入っている。外円の面積から中円・小円の面積を引くと 120 歩である。また，中円と小円の直径の差は 5 寸である。小円の直径を求めよ。
解析解は「算額（その834）」を参照のこと
外円の半径と中心座標を \(R, (0, 0)\)
中円の半径と中心座標を \(r_1, (0, R - r_1)\)
小円の半径と中心座標を \(r_2, (r_2, y)\)
とおき，以下の連立方程式を解く。
eq4 の 3 は，この算額で用いる円周率の近似値である(120 はこれに依存している）。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R::positive, r1::positive, r2::positive, y::negative;
@syms R, r1, r2, y;
eq1 = 2r1 - 2r2 - 5
eq2 = r2^2 + (R - r1 - y)^2 - (r1 + r2)^2 |> expand
eq3 = r2^2 + y^2 - (R - r2)^2 |> expand
eq4 = 3*(R^2 - (r1^2 + 2r2^2)) - 120 |> expand
res = solve([eq1, eq2, eq3, eq4], (R, r1, r2, y))
2-element Vector{NTuple{4, Sym}}:
(-sqrt(185)/2, 5/2, 0, -sqrt(185)/2)
(sqrt(185)/2, 5/2, 0, sqrt(185)/2)
連立方程式として解くと，小円の半径 = 0 という不適切な解を与える。
eq1, eq3, eq4 を連立方程式として解き，\(R, r_1, y\) を求め，eq2 に代入し eq2_2 とする。
eq2_2 は \(r_2\) の 6 次式となり，\(44x^6 + 440x^5 + 1420x^4 + 250x^3 - 19400x^2 - 107875x - 185000 = 0\) の解になる。
res = solve([eq1, eq3, eq4], (R, r1, y))[3] # 3 of 4
(sqrt(12*r2^2 + 20*r2 + 185)/2, (2*r2 + 5)/2, -sqrt(3*r2^2 - r2*sqrt(12*r2^2 + 20*r2 + 185) + 5*r2 + 185/4))
eq2_2 = eq2(R => res[1], r1 => res[2], y => res[3])
\(\displaystyle 6 r_{2}^{2} - r_{2} \left(2 r_{2} + 5\right) - r_{2} \sqrt{12 r_{2}^{2} + 20 r_{2} + 185} + 10 r_{2} - \frac{\left(2 r_{2} + 5\right) \sqrt{12 r_{2}^{2} + 20 r_{2} + 185}}{2} - \left(2 r_{2} + 5\right) \sqrt{3 r_{2}^{2} - r_{2} \sqrt{12 r_{2}^{2} + 20 r_{2} + 185} + 5 r_{2} + \frac{185}{4}} + \sqrt{12 r_{2}^{2} + 20 r_{2} + 185} \sqrt{3 r_{2}^{2} - r_{2} \sqrt{12 r_{2}^{2} + 20 r_{2} + 185} + 5 r_{2} + \frac{185}{4}} + \frac{185}{2}\)
res = solve(eq2_2, r2)[1] # 1 of 3
\(\displaystyle \operatorname{CRootOf} {\left(44 x^{6} + 440 x^{5} + 1420 x^{4} + 250 x^{3} - 19400 x^{2} - 107875 x - 185000, 1\right)}\)
実数解は \(r_2 = 3.90671194097833\) である。
小円の直径は 7.81342388195666 である。
res.evalf() |> println
res.evalf()*2 |> println
3.90671194097833
7.81342388195666
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 = 3.90671194097833
(R, r1, y) = (sqrt(12*r2^2 + 20*r2 + 185)/2, (2*r2 + 5)/2, -sqrt(3*r2^2 - r2*sqrt(12*r2^2 + 20*r2 + 185) + 5*r2 + 185/4))
@printf("R = %.6f; r1 = %.6f; r2 = %.6f; y = %.6f\n", R, r1, r2, y)
@printf("小円の直径 = %g\n", 2r2)
plot()
circle(0, 0, R, :black)
circle(0, R - r1, r1)
circle(r2, y, r2, :blue)
circle(-r2, y, r2, :blue)
if more
point(0, R - r1, " R-r1")
point(R, 0, " R")
point(r2, y, "(r2,y)")
annotate!(0.8, -4, text("小円の直径 = " * @sprintf("%.5f", 2r2), :left, 9))
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
end
end;
draw(true)
以下のアイコンをクリックして応援してください
