https://sangaku0418.hatenablog.com/entry/2024/07/31/180947

三八埼玉県熊谷市新堀諏訪大神社弘化4年(1847)

埼玉県立図書館：埼玉県史料集第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円5個，外円，斜線
#Julia #SymPy #算額 #和算 #数学

外円の中に原点の中心を原点として，\(x\) 軸と \(y\) 軸を線対称とする 4 本の斜線を引き，区画された領域に等円 4 個を容れる。外円の直径を 1.6 寸としたとき，等円が最も大きくなるときの等円の直径を求めよ。

外円の半径と中心座標を \(R, (0, 0)\)
等円の半径と中心座標を \(r, (x, y)\)
とおき，以下の連立方程式を立てる。

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

using SymPy

@syms R::positive, r::positive, x::positive, y::positive
eq1 = x^2 + y^2 - (R - r)^2
eq2 = r^2*(R^2 + y^2) - x^2*y^2;

eq1 は \(x^2 = (R - r)^2 - y^2\) なので，eq2 の \(x^2\) に代入して整理する。

eq2 = eq2(x^2 => (R - r)^2 - y^2) |> expand
eq2 |> println

R^2*r^2 - R^2*y^2 + 2*R*r*y^2 + y^4

eq2 は，\(R, y\) で \(r\) が決まるという式なので，\(r\) について解けば \(R, y\) の関数になる。

f = solve(eq2, r)[1]
f |> println

y*(R - y)/R

これは単純に \(R\) が定数のときは放物線の式であり，\(y\) がどのような値を取ると \(f(y)\) が最大になるかは計算しなくてもわかる。
たとえば \(R = 1.6/2\) のとき \(y = 0.4\) のとき \(r = 0.2\) が最大値になる。

pyplot(size=(300, 200), grid=true, aspectratio=:none, label="")
plot(f(R => 1.6/2), xlims=(0, 0.8), xlabel="y", ylabel="r")

r が最大になるときの y，そして y がその値を取るときの r は以下のようにして求める。

g = diff(f, y)
y0 = solve(g, y)[1]
y0 |> println

R/2

r0 = f(y => y0)
r0 |> println

R/4

\(R = 1.6/2\) 寸のとき，\(y_0 = 0.4\) 寸, \(r_0 = 0.2\) 寸ゆえ，等円の直径は 0.4 寸である。

@syms x0, y0, x, y, r, R
y0 = sqrt(R^2 - x0^2)
x = sqrt( (R - r)^2 - y^2)
eq = dist2(-R, 0, x0, y0, x, y, r)
res = solve(eq, x0)[1] |> println

(R*(8*R^4 - 16*R^3*r + 8*R^3*sqrt(R^2 - 2*R*r + r^2 - y^2) + 8*R^2*r^2 - 8*R^2*r*sqrt(R^2 - 2*R*r + r^2 - y^2) - 8*R^2*y^2 + 4*R*r*y^2 - 4*R*y^2*sqrt(R^2 - 2*R*r + r^2 - y^2) - r^4 + 2*r^2*y^2) - 4*sqrt(2)*r*y*sqrt(R^3*(4*R^3 - 8*R^2*r + 4*R^2*sqrt(R^2 - 2*R*r + r^2 - y^2) + 5*R*r^2 - 4*R*r*sqrt(R^2 - 2*R*r + r^2 - y^2) - 3*R*y^2 - r^3 + r^2*sqrt(R^2 - 2*R*r + r^2 - y^2) + r*y^2 - y^2*sqrt(R^2 - 2*R*r + r^2 - y^2))))/(8*R^4 - 16*R^3*r + 8*R^3*sqrt(R^2 - 2*R*r + r^2 - y^2) + 12*R^2*r^2 - 8*R^2*r*sqrt(R^2 - 2*R*r + r^2 - y^2) - 4*R^2*y^2 - 4*R*r^3 + 4*R*r^2*sqrt(R^2 - 2*R*r + r^2 - y^2) + r^4)

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

function draw(R, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(y, r) = (R/2, R/4)
x = sqrt( (R - r)^2 - y^2)
t = sqrt(R^2 - 2*R*r + r^2 - y^2)
u1 = R*(8R^4 - 16R^3*r + 8R^3*t + 8R^2*r^2 - 8R^2*r*t - 8R^2*y^2 + 4R*r*y^2 - 4R*y^2*t - r^4 + 2r^2*y^2)
u2 = 4√2r*y*sqrt(R^3*(4R^3 - 8R^2*r + 4R^2*t + 5R*r^2 - 4R*r*t - 3R*y^2 - r^3 + r^2*t + r*y^2 - y^2*t))
v = 8R^4 - 16R^3*r + 8R^3*t + 12R^2*r^2 - 8R^2*r*t - 4R^2*y^2 - 4R*r^3 + 4R*r^2*t + r^4
x0 = (u1 - u2)/v
y0 = sqrt(R^2 - x0^2)
plot()
circle(0, 0, R, :red)
circle4(x, y, r, :blue)
segment(-R, 0, x0, y0, :green)
segment(-R, 0, x0, -y0, :green)
segment(R, 0, -x0, y0, :green)
segment(R, 0, -x0, -y0, :green)
#=
abline(0, y, -y/R, -R, R, :green)
abline(0, y, y/R, -R, R, :green)
abline(0, -y, -y/R, -R, R, :green)
abline(0, -y, y/R, -R, R, :green)
=#
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(x, y, "等円:r,(x,y)", :blue, :center, delta=-delta/2)
point(0, y, "y", :blue, :center, delta=-delta/2)
point(0, R, "R", :blue, :center, :bottom, delta=delta/2)
point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
end
end;