新潟県長岡市 蒼柴神社 享和元年(1801)
和算の館
http://www.wasan.jp/niigata/aoshi.html
涌田和芳，外川一仁: 長岡蒼柴神社の算額，長岡工業高等専門学校研究紀要，第42巻，第2号(2006)
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_41-45/vol_42_2/42_2_1wakuta.pdf
キーワード：円5個，外円，弦，斜線
#Julia #SymPy #算額 #和算 #数学
外円内に，弦と 2 本の斜線を引き，区分された領域に甲円と乙円を 2 個ずつ入れる。大円の直径が 521 寸で，弦と矢の長さの差を最大にするとき，乙円の直径を求めよ。
注：矢（し）とは，弓形の孤の中点から弦におろした垂線
外円の半径と中心座標を \(R,\ (0,\ 0)\)
弦，矢の長さおよびその差を（そのまま）弦，矢，弦矢差
甲円の半径と中心座標を \(r_1,\ (0,\ R - r_1), (0,\ R - 矢 + r_1)\)
乙円の半径と中心座標を \(r_2,\ (x_2,\ y_2)\)
とおく。
まず，弦と矢の長さの差が最大になるときの矢と弦を決定する。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R::positive, 弦::positive, 矢::positive, x::positive,
r1::positive, r2::positive, x2::positive, y2::positive, d
# R = 521//2
弦 = 2sqrt(R^2 - (R - 矢)^2)
弦矢差 = 弦 - 矢
\(- 矢 + 2 \sqrt{R^{2} - \left(R - 矢\right)^{2}}\)
たとえば，\(R = 521/2\) のとき，図を描いてみると，矢が 150 前後のときに弦矢差が最大になることがわかる。
弦矢差(R => 521/2)
\(- 矢 + 521.0 \sqrt{1 - \left(1 - 0.00383877159309021 矢\right)^{2}}\)
pyplot(size=(300, 150), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(弦矢差(R => 521/2), xlims=(0, 400), xlabel="矢", ylabel="弦 - 矢")
最大値を取るときの矢の値を求めるには，導関数を求め，導関数が 0 になるときの値を求めればよい。
g = diff(弦矢差, 矢)
\(\displaystyle \frac{2 \left(R - 矢\right)}{\sqrt{R^{2} - \left(R - 矢\right)^{2}}} - 1\)
# 矢弦差が最大になるときの矢
mx_a = solve(g, 矢)[1]
@show(mx_a)
mx_a = R*(5 - sqrt(5))/5
\(\displaystyle \frac{R \left(5 - \sqrt{5}\right)}{5}\)
mx_a(R => 521/2).evalf()
\(144.000858372261\)
# 弦と矢の差
弦矢差(矢 => mx_a) |> simplify
\(R \left(-1 + \sqrt{5}\right)\)
弦矢差(矢 => mx_a)(R => 521/2).evalf()
\(321.995708138695\)
# 弦
弦(矢 => mx_a) |> simplify
\(\displaystyle \frac{4 \sqrt{5} R}{5}\)
弦(矢 => mx_a)(R => 521/2).evalf()
\(465.996566510956\)
外円の半径が，\(R = 521/2\) のとき，
矢が 144.000858372261 のとき，弦と矢の差が最大値 321.995708138695 になる。ちなみにそのときの弦は 465.996566510956 である。
弦と矢の差が最大となるときの矢，弦が決まったので，その状況における乙円の半径を求める。
算額の問では乙円の半径を要求しているが，甲円の半径は無関係である。
using SymPy
@syms 矢, 弦, R::positive, r2::positive
s5 = √Sym(5)
矢 = R*(5 - s5)/5
弦 = 2sqrt(R^2 - (R - 矢)^2)
AB = sqrt( (弦/2)^2 + (矢/2)^2) |> simplify |> sympy.sqrtdenest
\(\displaystyle \frac{R \sqrt{110 - 10 \sqrt{5}}}{10}\)
BD = sqrt( (R - r2)^2 - (R - 矢/2)^2) |> simplify # OD^2 - OB^2
@show(BD)
BD = sqrt(-R^2*(sqrt(5) + 5)^2 + 100*(R - r2)^2)/10
\(\displaystyle \frac{\sqrt{- R^{2} \left(\sqrt{5} + 5\right)^{2} + 100 \left(R - r_{2}\right)^{2}}}{10}\)
eq2 = (矢/2)/AB ⩵ r2/BD
\(\displaystyle \frac{5 - \sqrt{5}}{\sqrt{110 - 10 \sqrt{5}}} = \frac{10 r_{2}}{\sqrt{- R^{2} \left(\sqrt{5} + 5\right)^{2} + 100 \left(R - r_{2}\right)^{2}}}\)
ans_r2 = solve(eq2, r2)[1] |> sympy.sqrtdenest |> simplify
@show(ans_r2)
ans_r2 = R*(5 - sqrt(5))/20
\(\displaystyle \frac{R \left(5 - \sqrt{5}\right)}{20}\)
ans_r2(R => 521/2).evalf() *2 |> println
72.0004291861305
乙円の半径 \(r_2\) は，外円の半径 \(R\) の \(\displaystyle \frac{5 - \sqrt{5}}{20} = 0.138196601125011\) 倍である。
外円の直径が 521 寸のとき，乙円の直径は 72.0004291861305 寸である。
矢と弦がわかれば，甲円の半径は以下の eq1, eq2 の連立方程式を解くことで得られる。
しかし，式が複雑すぎるので，数値解を求めることにする。
@syms R, x
矢 = (5 - √Sym(5))R/5
弦 = 2sqrt(R^2 - (R - 矢)^2)
y = sqrt(R^2 - x^2)
eq1 = dist(x, sqrt(R^2 - x^2), -弦/2, R - 矢, 0, R - r1) - r1^2;
eq2 = dist(x, sqrt(R^2 - x^2), -弦/2, R - 矢, 0, R - 矢 + r1) - r1^2;
function driver(R)
# r1, x は R の倍数にすぎない
# (r1, x) = R.*(0.13504615662994296, 0.4861961875171891)
function H(u)
(r1, x) = u
a = sqrt(R^2 - x^2)
b = 1 - 1/√5
c = -x - sqrt(R^2 - (R - R*b)^2)
d = R - R*b - a
e = c^2 + d^2
f = d + r1
g = R - r1 - a
h = g*d - x*c
i = f*d - x*c
return [
(x + c*h/e)^2 + (g - h*d/e)^2 - r1^2, # eq1
(x + c*i/e)^2 + (f - i*d/e)^2 - r1^2 # eq2
]
end
iniv = R.*BigFloat[0.135, 0.486]
res = nls(H, ini=iniv)
return res[1]
end;
# driver(1) # (0.13504615662994296, 0.4861961875171892)
甲円の半径 \(r_1\) は 外円の半径 \(R\) の 0.13504615662994296 倍である。
外円の直径が 521 寸のとき，甲円の直径は 70.35904760420028 寸である。
描画関数プログラムのソースを見る
function draw(R, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
# (r1, x) = driver(R)
(r1, x) = R.*(0.13504615662994296, 0.4861961875171892)
矢 = (5 - √5)R/5
弦 = 2sqrt(R^2 - (R - 矢)^2)
y = sqrt(R^2 - x^2)
r2 = R*(5 - √5)/20
x2 = sqrt(-R^2*(√5 + 5)^2 + 100*(R - r2)^2)/10
y2 = R - 矢/2
@printf("弦 = %g; 矢 = %g; 弦と矢の差 = %g\n", 弦, 矢, 弦 - 矢)
@printf("乙円の直径 = %.15g\n", 2r2)
@printf("r1 = %g; r2 = %g; x2 = %g; y2 = %g; x = %g; y = %g\n", r1, r2, x2, y2, x, y)
plot()
circle(0, 0, R)
segment(-弦/2, R - 矢, 弦/2, R - 矢, :blue)
circle2(x2, y2, r2, :green)
circle(0, R - r1, r1, :magenta)
circle(0, R - 矢 + r1, r1, :magenta)
segment(x, y, -弦/2, R - 矢)
segment(-x, y, 弦/2, R - 矢)
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(0, R - 矢, "R-矢", :blue, :center, delta=-delta/2)
point(x, y, "(x,y)", :green, :left, :bottom, delta=delta/2)
point(弦/2, R - 矢, "弦/2 ", :blue, :right, delta=-delta/2)
point(0, R - r1, "甲円:r1,(0,R-r1)", :black, :center, :bottom, delta=delta/2)
point(0, R - 矢/2, "R-矢/2 ", :black, :right, :vcenter)
point(0, R - 矢 + r1, "甲円:r1,(0,R-矢+r1)", :black, :center, delta=-delta/2)
point(x2, y2, "乙円:r2,(x2,y2)", :black, :center, delta=-delta/2)
point(0, R, "R", :red, :center, :bottom, delta=delta/2)
end
end;
#draw(521/2, true)
draw(1, true)
弦 = 1.78885; 矢 = 0.552786; 弦と矢の差 = 1.23607
乙円の直径 = 0.276393202250021
r1 = 0.135046; r2 = 0.138197; x2 = 0.468079; y2 = 0.723607; x = 0.486196; y = 0.87385
以下のアイコンをクリックして応援してください
