算額（その637）
長野県上田市別所温泉 北向観音堂 文政11年(1828)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
2つの楕円が交差する中に等円が 4 個入っている。楕円の長径（差渡し）が 4 寸 9 分のとき，等円の直径が最大になるときの楕円の短径（差渡し）はいかほどか。
内接円の中心間の距離(eq1)，楕円と円が共通接点を持つ(eq2, eq3)，短径と円の中心の関係(eq4) があるが，eq1 〜 eq4 を一度に解くことができない。
それは，等円の半径は楕円の短径に依存するということが理由である。
そこで，例えば eq2, eq3, eq4 を連立させて r, x, y を求める。それぞれの式には b が含まれる。つまり，r, x, y は b により決定されるということである。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, b::positive, r::positive,
x0::positive, x::positive, y::positive
a = 49//20
xx = sqrt((a^2 - b^2)*(b^2 - r^2))/b
eq1 = -b^2/a^2*x/y - (-(x - xx)/y)
eq2 = x^2/a^2 + y^2/b^2 - 1
eq3 = (x - xx)^2 + y^2 - r^2
eq4 = b + r - xx;
res = solve([eq2, eq3, eq4], (r, x, y))
2-element Vector{Tuple{Sym, Sym, Sym}}:
(-800*b^3/2401 + b, -sqrt(49 - sqrt(2401 - 1600*b^2))*sqrt(sqrt(2401 - 1600*b^2) + 49)/20, b*sqrt(2401 - 1600*b^2)/49)
(-800*b^3/2401 + b, sqrt(49 - sqrt(2401 - 1600*b^2))*sqrt(sqrt(2401 - 1600*b^2) + 49)/20, b*sqrt(2401 - 1600*b^2)/49)
2 組の解が得られる。本質的にはどちらでもよいが，当初の予想と同じ符号になる 2 番目のものを適解とする。
r として得られた解をプロットしてみると，b = 1 近辺で最大値になることがわかる。
r = res[2][1]
pyplot(size=(400, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(r, xlims=(0.95, 1.05), xlab="b", ylab="r")
r が最大になるのは b = 49*sqrt(6)/120 = 1.00020831163646 のときである。
b0 = solve(diff(res[1][1], b))[1]
(b0, b0.evalf()) |> println
(49*sqrt(6)/120, 1.00020831163646)
b = 49*sqrt(6)/120 = 1.00020831163646 のときの r, x, y は以下の通りである。
r0 = res[2][1](b => b0)
x0 = res[2][2](b => b0)
y0 = res[2][3](b => b0)
(r0, x0, y0) |> println
(r0.evalf(), x0.evalf(), y0.evalf()) |> println
(49*sqrt(6)/180, sqrt(49 - 49*sqrt(3)/3)*sqrt(49*sqrt(3)/3 + 49)/20, 49*sqrt(2)/120)
(0.666805541090976, 2.00041662327293, 0.577470537969014)
その他のパラメータは以下の通り。
a = 2.45; b = 1.00021; r = 0.666806; xx = 2.00042; x = 2.00042; y = 0.577471
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
a = 4.9/2
b = 1.00020831163646
(r, x, y) = (0.666805541090976, 2.00041662327293, 0.577470537969014)
xx = sqrt((a^2 - b^2)*(b^2 - r^2))/b # 1.67
@printf("a = %g; b = %g; r = %g; xx = %g; x = %g; y = %g\n", a, b, r, x0, x, y)
plot()
ellipse(0, 0, a, b, color=:red)
ellipse(0, 0, b, a, color=:red)
circle42(0, xx, r, :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(xx, 0, "xx", :blue, :center, delta=-delta/2)
point(x, y, "(x,y)", :blue, :left, :bottom, delta=delta/2)
point(b, 0, "b ", :red, :center, delta=-delta/2)
point(a, 0, " a", :red, :center, delta=-delta/2)
end
end;
