長野県小諸市諸 金毘羅社 文政7年(1824)
中村信弥「改訂増補 長野県の算額」県内の算額(P.106)
http://www.wasan.jp/zoho/zoho.html
キーワード：外円，楕円
#Julia #SymPy #算額 #和算 #数学
円内に等楕円が 3 個入っている。円の直径が 99 寸のとき，楕円の面積が最大になるときの楕円の長径を求めよ。
円の半径と中心座標を \(R,\ (0,\ 0)\)
楕円の長径，短径と中心座標を \(a,\ b,\ (0,\ y)\)
円と楕円の共通接点を \( (x_1,\ y_1)\)
上と右の楕円の共通接点を \( (x_2,\ y_2)\)
とおき，以下の連立方程式を解く。
eq1, eq2, eq3 は円と楕円が共通接点 \( (x_1,\ y_1)\) を持つこと，
eq4, eq5 は 2 つの楕円が共通接点 \( (x_2,\ y_2)\) を持つことを意味している。
これらの方程式群から，楕円の中心の \(y\) 座標を求める。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms x1::positive, y1::positive, x2::positive, y2::positive, y::positive, a::positive, b::positive, R::positive
y2 = x2/sqrt(Sym(3))
eq1 = x1^2/a^2 + (y1 - y)^2/b^2 - 1
eq2 = b^2*x1/(a^2*(y - y1)) + x1/y1
eq3 = x1^2 + y1^2 - R^2
eq4 = x2^2/a^2 + (y2 - y)^2/b^2 - 1
eq5 = b^2*x2/(a^2*(y - y2))- 1/sqrt(Sym(3));
res1 = solve([eq1, eq2, eq3], (x1, y1, y))[1] # 1 of 2
(-sqrt(-(R^2*b^2 - a^4)/(a - b))/sqrt(a + b), a*sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/(a^2 - b^2), sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/a)
res2 = solve([eq4, eq5], (x2, y))[1]
(a^2/sqrt(a^2 + 3*b^2), sqrt(3*a^2 + 9*b^2)/3)
両方の解 \(y\) が等しいとすると以下のようになる。
eq = res1[3] - res2[2]
eq |> println
-sqrt(3*a^2 + 9*b^2)/3 + sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/a
方程式を \(b\) について解く。
res3 = solve(eq, b)[1]
res3 |> println
a*sqrt(9*R^2 - 12*a^2)/(3*R)
楕円の面積 \(S\) は長軸の長さ \(a\) の関数になり，以下のような外形である。
S = PI * a * a*sqrt(9*R^2 - 12*a^2)/(3*R);
円の半径を \(R = 99/2\) として，図を描いてみる。
\(a\) が 30〜40 の間で面積が最大になることがわかる。
pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(S(R => 99/2), xlims=(0, 43), xlabel="長径 a", ylabel="楕円の面積 S")
\(S\) の導関数を求め，接線の傾きが \(0\) になるときの \(a\) を求める。
max_at_a = solve(diff(S, a), a)
max_at_a[1] |> println
sqrt(2)*R/2
\(a\) の値が 35.0017856687341 のとき，面積が最大値 8888.52378442978 になる。
max_at_a[1](R => 99/2).evalf() |> println
35.0017856687341
S(R => 99/2, a => 35.0017856687341).evalf() |> println
2222.13094610745
そのときの短径 \(b\) は 20.2082903779612 である。
(a*sqrt(9*R^2 - 12*a^2)/(3*R))(R => 99/2, a => 35.0017856687341).evalf() |> println
20.2082903779612
算額において，長径，短径は差渡し径なので，現代の用語で表すものの 2 倍である。
よって，算額の答えとしては「長径は 70.0035713374682，短径は 40.4165807559224」である。
その他のパラメータは以下の通り。
\(x_1 = 24.75;\ y_1 = 42.8683;\ x_2 = 24.75;\ y_2 = 14.2894\)
\(y = 28.5788;\ a = 35.0018;\ b = 20.2083\)
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
R = 99/2
a = √2R/2
b = a*sqrt(9*R^2 - 12*a^2)/(3*R)
(x1, y1, y) = (sqrt( (-R^2*b^2 + a^4)/(a - b))/sqrt(a + b),
a*sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/(a^2 - b^2),
sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/a)
x2 = a^2/sqrt(a^2 + 3*b^2)
y2 = x2/sqrt(3)
@printf("x1 = %g; y1 = %g; x2 = %g; y2 = %g; y = %g; a = %g; b = %g\n", x1, y1, x2, y2, y, a, b)
plot()
circle(0, 0, R)
ellipse(0, y, a, b, color=:blue)
ellipse(√3y/2, -y/2, a, b, φ=240, color=:blue)
ellipse(-√3y/2, -y/2, a, b, φ=300, color=:blue)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
point(x1, y1, " (x1,y1)", :blue, :left, :bottom)
point(x2, y2, " (x2,y2)", :blue, :left)
point(0, y, " y", :blue, :left, :vcenter)
end
end;
以下のアイコンをクリックして応援してください
