13 岩手県奥州市（旧水沢市佐倉河） 胆沢城跡鎮守府八幡宮 弘化2年(1845)
安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
一関博物館 和算に挑戦 令和5年度出題問題(3)［上級問題］（高校生・一般向き）
https://www.city.ichinoseki.iwate.jp/museum/wasan/r5/hard.html
キーワード：円2個，外円，楕円
#Julia #SymPy #算額 #和算 #数学
大円の中に楕円 3 個と小円を容れる。楕円の面積が最大になるとき，大円の直径を小円の直径で表す術を述べよ。
この算額は「算額（その678）」と求めるものが違うだけで，本質的に同じ問題である。
変数値が正の値を取るように，上下を反転させて解く。
大円の半径と中心座標を \(R, (0, 0)\)
小円の半径と中心座標を \(r, (0, 0)\)
楕円の超半径，短半径，中心座標を \(a, b, (0, y)\)
楕円と円の接点座標を \( (x_1, y_1)\)
楕円同士の接点座標を \( (x_2, y_2)\)
とおき，以下の連立方程式を解く。
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));
eq1, eq2, eq3 を解き，\(x_1, y_1, y\) を求める。
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)
eq4, eq5 を解き，\(x_2, y\) を求める。
res2 = solve([eq4, eq5], (x2, y))[1]
(a^2/sqrt(a^2 + 3*b^2), sqrt(3*a^2 + 9*b^2)/3)
res1 と res2 の \(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\) を求める。
ans_b = solve(eq, b)[1]
ans_b |> println
a*sqrt(9*R^2 - 12*a^2)/(3*R)
楕円の面積 \(S\) は \(\pi a b\) である。
S = PI * a * ans_b;
面積が最大になるときの長半径を求めるには，\(S(a)\) の導関数 \(S(a)\)'を求め，\(S(a)'= 0\) を解く。
max_at_a = solve(diff(S, a), a)[1]
max_at_a |> println
sqrt(2)*R/2
楕円の面積が最大になるときの長半径は \(\sqrt{2}R/2\)，
ans_b = ans_b(a => max_at_a)
ans_b |> println
sqrt(6)*R/6
同じく短半径は \(\sqrt{6}R/6\)，
ans_y = res2[2](b =>ans_b)(a => max_at_a)
ans_y |> println
sqrt(3)*R/3
そのときの，中心座標の \(y\) 座標値は \(\sqrt{3}R/3\) である。
小円の半径は \(r = ans_y - ans_b = (\sqrt{3}/3 - \sqrt{6}/6)R\) で，
\(R = r/(\sqrt{3}/3 - \sqrt{6}/6) = r(\sqrt{6} + 2\sqrt{3})\) である。
すなわち，大円の直径は小円の直径の \( \sqrt{6} + 2\sqrt{3} = 5.913591357920932\) 倍である。
@syms r, d
apart(r / (√Sym(3)/3 - √Sym(6)/6), d) |> simplify |> factor |> println
r*(sqrt(6) + 2*sqrt(3))
術は \(\sqrt{\sqrt{288} + 18}\) である。二重根号を外せば \(\sqrt{\sqrt{288} + 18} = \sqrt{6} + 2\sqrt{3}\) で前述の式と同じになる。
√(√Sym(288) + 18) |> sympy.sqrtdenest |> println
sqrt(6) + 2*sqrt(3)
描画関数プログラムのソースを見る
function draw(r, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
R = r*(√6 + 2√3)
a = √2R/2
b = √6R/6
y = √3R/3
x1 = sqrt( (-R^2*b^2 + a^4)/(a - b))/sqrt(a + b)
y1 = a*sqrt(R^2*a^2 - R^2*b^2 - a^4 + a^2*b^2)/(a^2 - b^2)
x2 = a^2/sqrt(a^2 + 3*b^2)
y2 = y2 = √3x2/3
println( (x1, y1))
println( (x2, y2))
plot()
circle(0, 0, R)
circle(0, 0, r, :green)
ellipse(0, y, a, b, color=:blue)
ellipse(y*cosd(30), -y*sind(30), a, b, color=:blue, φ=-120)
ellipse(-y*cosd(30), -y*sind(30), a, b, color=:blue, φ=120)
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(x1, y1, "(x1,y1)", :green, :left, :bottom, delta=delta/2)
point(x2, y2, "(x2,y2)", :green, :left, delta=-delta/2)
point(0, R, "R", :red, :center, :bottom, delta=delta/2)
point(0, r, "r", :green, :center, :bottom, delta=delta/2)
point(0, y, " y=r+b", :green, :center, :bottom, delta=delta/2)
end
end;
draw(1/2, true)
以下のアイコンをクリックして応援してください
