算額（その775）
福島県二本松市 新田神社 元治元年(1864)
http://www.wasan.jp/fukusima/sinden.html
正方形内に交差する 2 個の楕円と，区画された領域に甲円 1 個，乙円 4 個が入っている。甲円の直径が 4 寸，乙円の直径が 3 寸のとき，正方形の面積はいかほどか。
計算を簡単にするために図形を 45 度回転し楕円の長軸・短軸が x-y 軸に載るようにする。
楕円の長半径と短半径と中心座標をそれぞれ a, b, (0, 0)
甲円の半径と中心座標を r1, (0, 0); 甲円の半径は楕円の短半径に等しい。r1 = b
乙円の半径と中心座標を r2, (b + r2, 0)
楕円と乙円の交点座標を (x0, y0)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms a::positive, b::positive, x0::positive, y0::positive, r1, r2, x1, y1
r2 = 3//2
r1 = b = 4//2
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*x0/(a^2*y0) + (x0 - b - r2)/y0
eq3 = (x0 - b - r2)^2 + y0^2 - r2^2
res1 = solve([eq1, eq2, eq3], (a, x0, y0));
res1[4]
(sqrt(2)*b^(3/2)*sqrt(1/(b - r2)), 2*b, sqrt(b)*sqrt(-b + 2*r2))
4 組の解が得られるが，最後のものが適解である。
楕円が内接する正方形は，楕円と正方形の一辺の接点を (x1, y1) とした以下の連立方程式を解く。
@syms x1::positive, y1::positive
r2 = 3//2
r1 = b = 4//2
a = sqrt(Sym(2))*b^(3//2)*sqrt(1/(b - r2))
eq4 = x1^2/a^2 + y1^2/b^2 - 1
eq5 = -b^2*x1/(a^2*y1) + 1
res2 = solve([eq4, eq5], (x1, y1));
res2[2]
(2*b^2*sqrt((b - r2)/(3*b - r2))/(b - r2), b*sqrt((b - r2)/(3*b - r2)))
2 組の解が得られるが，2 番目のものが適解である。
接点を通る傾き -1 の直線が x-y 軸で切り取られるものが正方形の一辺 sl である。x2 は x 切片（x 軸との交点）。
x2 = x1 + y1
sl = x2*√2
甲円の直径が 4 寸，乙円の直径が 3 寸のとき，正方形の一辺の長さは 8.48528，正方形の面積は 72 歩である。
sl = sqrt(Sym(2)) * (res2[2][1] + res2[2][2]) |> simplify
sl |> println
sl(b => 4//2, r2 => 3//2).evalf() |> println
sl(b => 4//2, r2 => 3//2)^2 |> println
sqrt(2)*b*sqrt((b - r2)/(3*b - r2))*(3*b - r2)/(b - r2)
8.48528137423857
72
その他のパラメータは以下のとおりである。
r1 = 2; r2 = 1.5; a = 5.65685; b = 2; x0 = 4; y0 = 1.41421; x1 = 5.33333; y1 = 0.666667; x2 = 6
もとの位置で図を描くのも簡単だが，わかりやすくするために回転させたままの図を描く。
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = b = 4//2
r2 = 3//2
(a, x0, y0) = (sqrt(2)*b^(3/2)*sqrt(1/(b - r2)), 2*b, sqrt(b)*sqrt(-b + 2*r2))
(x1, y1) = (2.0*b^2*sqrt((b - r2)/(3.0*b - r2))/(b - r2), b*sqrt((b - r2)/(3.0*b - r2)))
x2 = x1 + y1
sl = x2*√2
@printf("正方形の一辺の長さは %g，面積は %g\n", sl, sl^2)
@printf("r1 = %g; r2 = %g; a = %g; b = %g; x0 = %g; y0 = %g; x1 = %g; y1 = %g; x2 = %g\n",
r1, r2, a, b, x0, y0, x1, y1, x2)
plot([0, x2, 0, -x2, 0], [-x2, 0, x2, 0, -x2], color=:blue, lw=0.5)
ellipse(0, 0, a, b, color=:green)
ellipse(0, 0, b, a, color=:green)
circle(0, 0, r1, :magenta)
circle42(0, b + r2, r2)
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, 0, "甲円:r1\n(0,0)", :magenta, :center, delta=-delta)
point(b + r2, 0, "乙円:r2\n(b+r2,0)", :red, :center, delta=-delta)
point(x0, y0, "(x0,y0)", :green, :right, delta=-delta)
point(x1, y1, " (x1,y1)", :green, :left, :vcenter)
point(x2, 0, " x2", :blue, :left, delta=-delta/2)
end
end;
