算額（その311）

算額（その311）

「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県菰野町 伎留太神社 寛政9年(1797)

問題文1
扇形の中に直径の等しい 3 個の円が入っている。条件として，弦（扇形を形成する左右の辺の端を結んだ線分），矢（円弧の中点から弦におろした垂線），左右斜（扇形の半径から扇形と中心が一致する内部の円の半径を引いた線分）が与えられているとき円の直径を求めよ。

扇の中心から扇の先端までの距離（扇の外円の半径）を R，等円の半径を r，右下の等円の中心座標を(r, R - 矢) として，以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, r::positive, y1::positive,
   左右斜::negative, 矢::positive, 弦::positive;

(左右斜, 矢, 弦) = (50, 40, 130)

eq1 = r^2 + (R - r - y1)^2 - 4r^2
eq2 = r^2 + y1^2 - (r + (R - 左右斜))^2
eq3 = (弦/2)^2 - (R^2 - (R - 矢)^2);

eq1 を y1 について解く。

solve(eq1, y1) |> println

   Sym[R - r + sqrt(3)*r, R - sqrt(3)*r - r]

eq2 の y1 に代入する。

eq22 = eq2(y1 => R - r - sqrt(Sym(3))r) |> expand
eq22 |> println

   -4*R*r - 2*sqrt(3)*R*r + 2*R*左右斜 + 2*sqrt(3)*r^2 + 4*r^2 + 2*r*左右斜 - 左右斜^2

eq22 を R について解く。

solve(eq22, R)[1] |> println

   (sqrt(3)*r^2 + 2*r^2 + r*左右斜 - 左右斜^2/2)/(sqrt(3)*r + 2*r - 左右斜)

eq3 に代入する。

eq32 = eq3(R => (sqrt(Sym(3))*r^2 + 2*r^2 + r*左右斜 - 左右斜^2/2)/(sqrt(Sym(Sym(3)))*r + 2*r - 左右斜))
eq32 |> println

   弦^2/4 + (-矢 + (sqrt(3)*r^2 + 2*r^2 + r*左右斜 - 左右斜^2/2)/(sqrt(3)*r + 2*r - 左右斜))^2 - (sqrt(3)*r^2 + 2*r^2 + r*左右斜 - 左右斜^2/2)^2/(sqrt(3)*r + 2*r - 左右斜)^2

左右斜, 矢, 弦 を代入すると r のみの式になる。

eq33 = eq32(左右斜 => 50, 矢 => 40, 弦 => 130)

   (-40 + (sqrt(3)*r^2 + 2*r^2 + 50*r - 1250)/(sqrt(3)*r + 2*r - 50))^2 + 4225 - (sqrt(3)*r^2 + 2*r^2 + 50*r - 1250)^2/(sqrt(3)*r + 2*r - 50)^2

二通りの解（r）が得られる。

res = solve(eq33);
length(res)

   2

いずれも虚数として表示されるが，虚部は誤差範囲で 0 である。
したがって，r は 45.2629281760725，14.1521122023713 であるが，後者が適解である。

solve(eq33)[1].evalf(), solve(eq33)[2].evalf()

   (45.2629281760725 - 0.e-22*I, 14.1521122023713 + 0.e-21*I)

このあとさらに y1, R も求めなければならない。

---

最初から連立方程式に，左右斜, 矢, 弦 を代入して解を求めれば，R, r, y1 が一度に求まる。

using SymPy

@syms R::positive, r::positive, y1::positive,
   左右斜::negative, 矢::positive, 弦::positive;

(左右斜, 矢, 弦) = (50, 40, 130)

eq1 = r^2 + (R - r - y1)^2 - 4r^2
eq2 = r^2 + y1^2 - (r + (R - 左右斜))^2
eq3 = (弦/2)^2 - (R^2 - (R - 矢)^2);
res = solve([eq1, eq2, eq3], (R, r, y1))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (72.8125000000000, 14.1521122023713, 34.1482104287061)

   左右斜= 50;  矢 = 40;  弦 = 130
   R = 72.812500;  r = 14.152112;  y1 = 34.148210
   等円の直径は 2r = 28.304224

等円の直径は 2×14.1521122023713 = 28.3042244047426 である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r, y1) =  (72.8125000000000, 14.1521122023713, 34.1482104287061)
   @printf("左右斜= %d;  矢 = %d;  弦 = %d\n", 左右斜, 矢, 弦)
   @printf("R = %.6f;  r = %.6f;  y1 = %.6f\n", R, r, y1)
   @printf("等円の直径は 2r = %.6f\n",  2r)
   y = R - 矢
   x = sqrt(R^2 - y^2)
   θ = round(Int, atand(y/x))
   plot()
   circle(0, 0, R, beginangle=θ, endangle=180-θ)
   circle(0, 0, R - 左右斜, :blue, beginangle=θ, endangle=180-θ)
   circle(0, R - r, r, :green)
   circle(r, y1, r, :green)
   circle(-r, y1, r, :green)
   plot!([-x, 0, x, -x], [y, 0, y, y], color=:black, lw=0.5)
   if more
       point(0, R - r, " R-r", :green)
       point(0, y, " R-矢 ", :black, :right, :bottom)
       point(r, y1, " (r,y1)", :green, :left, :bottom)
       point(sqrt(R^2 - y^2), y, "(√(R^2-(R-矢)^2),R-矢)　", :black, :right)
       point(0, R, " R", :red, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;