一八 埼玉県大里郡岡部村岡 稲荷社 文化14年(1817)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円5個，楕円，斜線
#Julia #SymPy #算額 #和算 #数学
楕円の中に2本の斜線を引き，区分された領域に天円 2 個，地円 1 個，人円 2 個を容れる。楕円の長径，短径がそれぞれ 69 寸，23 寸のとき，人円の直径はいかほどか。
楕円の長半径，短半径を \(a,\ b\)
天円の半径と中心座標を \(r_1,\ (r_2 + r_1,\ 0)\)
地円の半径と中心座標を \(r_2,\ (0,\ 0)\)
人円の半径と中心座標を \(r_3,\ (r_2 - r_3,\ 0)\)
原点を通る，地円と人円の共通接線と楕円の交点座標を \( (x_0,\ y_0)\)
とおき，以下の連立方程式を解く。
1. 天円の半径
天円の半径は，「算法助術」の公式84により求めることができる。
楕円に内接する同じ大きさの 2 円の半径 \(r\) は，楕円の長半径 \(a\)，短半径 \(b\)，原点から円の中心までの距離 \(d\)で表すことができる。
\(\displaystyle d^2 = \frac{(a^2 - b^2)(b^2 - r^2)}{b^2}\)
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a::positive, b::positive,
x0::positive, y0::positive,
r1::positive, r2::positive, r3::positive,
x01::positive, y01::positive, d
r2 = b
d = r2 + r1
eq1 = d^2 - (a^2 - b^2)*(b^2 - r1^2)/b^2
r1 = solve(eq1, r1)[1]
r1 |> println
b - 2*b^3/a^2
\(\displaystyle b - \frac{2b^3}{a^2} \)
2. 人円の半径
次いで，\(r_1\) が確定したとして，人円 \(r_3\) を求める。これは，人円と天円の中心から斜線までの距離と，原点からそれぞれの中心までの距離の関係から導かれる 3元連立方程式を解くことで求めることができる。
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = r3/(r2 - r3) - y0/sqrt(x0^2 + y0^2)
eq3 = r1/(r2 + r1) - y0/sqrt(x0^2 + y0^2)
res = solve([eq1, eq2, eq3], (x0, y0, r3))[4]; # 4 of 4
# x0
res[1]
\(\displaystyle \frac{\sqrt{a \sqrt{a + b} - \left(a^{2} - 2 b^{2}\right) \sqrt{\frac{1}{a - b}}} \sqrt{a \sqrt{a + b} + \left(a^{2} - 2 b^{2}\right) \sqrt{\frac{1}{a - b}}}}{\sqrt{a + b}}\)
# y0
res[2]
\(\displaystyle \frac{b \left(a^{2} - 2 b^{2}\right) \sqrt{\frac{1}{a - b}}}{a \sqrt{a + b}}\)
# r3
res[3]
\(\displaystyle \frac{b \left(a^{2} - 2 b^{2}\right) \sqrt{\frac{1}{a - b}}}{a^{2} \sqrt{\frac{1}{a - b}} - 2 b^{2} \sqrt{\frac{1}{a - b}} + 2 \sqrt{\frac{a^{4} - 2 a^{2} b^{2} + b^{4}}{a^{2} - b^{2}}} \sqrt{a + b}}\)
res[3](a => 69/2, b => 23/2)
\(\displaystyle 3.5\)
楕円の長径，短径が それぞれ 69 寸，23 寸のとき，人円の直径は 7 寸である。
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, r2) = (69, 23) ./ 2
b = r2
r1 = b - 2*b^3/a^2
x0 = sqrt(a*sqrt(a + b) - (a^2 - 2*b^2)*sqrt(1/(a - b)))*sqrt(a*sqrt(a + b) + (a^2 - 2*b^2)*sqrt(1/(a - b)))/sqrt(a + b)
y0 = b*(a^2 - 2*b^2)*sqrt(1/(a - b))/(a*sqrt(a + b))
r3 = b*(a^2 - 2*b^2)*sqrt(1/(a - b))/(a^2*sqrt(1/(a - b)) - 2*b^2*sqrt(1/(a - b)) + 2*sqrt( (a^4 - 2*a^2*b^2 + b^4)/(a^2 - b^2))*sqrt(a + b))
@printf("長径が %g 寸，短径が %g 寸のとき，人円の直径は %g 寸である。\n", 2a, 2b, 2r3)
plot()
ellipse(0, 0, a, b, color=:green)
circle2(r2 + r1, 0, r1)
circle2(r2 - r3, 0, r3, :magenta)
circle(0, 0, r2, :blue)
segment(-x0, -y0, x0, y0, :orange)
segment(-x0, y0, x0, -y0, :orange)
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(r2 + r1, 0, "天円:r1\n(r2+r1,r1)", :red, :center, delta=-2delta)
point(0, 3b/4, "地円:r2\n(0,0)", :blue, :center, delta=-2delta, mark=false)
point(r2 - r3, 0, "人円:r3\n(r2-r3,r3)", :black, :center, delta=-2delta, deltax=-4delta)
point(a, 0, " a", :green, :left, :vcenter)
point(0, b, "b = r2", :blue, :center, :bottom, delta=2delta)
end
end;
以下のアイコンをクリックして応援してください
