算額（その1487）

算額（その1487）

四十七 岩手県一関市平沢 平沢白山神社 慶応2年(1866)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：二等辺三角形，楕円2個
#Julia, #SymPy, #算額, #和算

二等辺三角形の中に合同な楕円を 2 個容れる。短径が 1 寸のとき，二等辺三角形の高さが最も低くなるときの高さはいかほどか。

これは「算額（その699）」の特別な場合なので，同じように解ける。高さは楕円の長径の関数である。高さの導関数を取り，導関数=0 のときの高さを求めればよい。

二等辺三角形の底辺の長さを 2a，高さを h = a*s3 とする。パラメータとして s3 = h/a を使う。
上にある縦長楕円の長半径，短半径，中心座標を a1, b1, (0, 2b2 + b1)
下にある横長楕円の長半径，短半径，中心座標を a2, b2, (0, b2)
二等辺三角形の右側の斜辺と上の楕円，下の楕円との接点座標を (x1, y1), (x2, y2)
a1 = b2, b1 = a2
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a1::positive, b1::positive, x1::positive, y1::positive,
     a2::positive, b2::positive, x2::positive, y2::positive
@syms a::positive, s3::positive,  h::positive

b2 = 1//2
(a1, b1) = (b2, a2)
y1 = (a - x1)*s3
y2 = (a - x2)*s3
eq1 = x1^2/a1^2 + (y1 - (2b2 + b1))^2/b1^2 - 1
eq3 = x2^2/a2^2 + (y2 - b2)^2/b2^2 - 1
eq2 = -b1^2*x1/(a1^2*(y1 - (2b2 + b1))) + s3
eq4 = -b2^2*x2/(a2^2*(y2 - b2)) + s3;

res = solve([eq1, eq2, eq3, eq4], (s3, a, x1, x2))[2]  # 2 of 2

    (2*sqrt(4*a2^2 - 2*a2 + 1)/(2*a2 - 1), 2*a2^2/sqrt(4*a2^2 - 2*a2 + 1), sqrt(4*a2^2 - 2*a2 + 1)/(2*(2*a2^2 - a2 + 1)), 4*a2^2*sqrt(4*a2^2 - 2*a2 + 1)/(8*a2^2 - 2*a2 + 1))

# 高さ
h = res[1]*res[2]
h |> println

    4*a2^2/(2*a2 - 1)

高さは楕円の長半径 a2 の関数である。

# 導関数
g = diff(h, a2) |> factor

    8*a2*(a2 - 1)/(2*a2 - 1)^2

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="")
plot(g, xlims=(0.8, 1.2), xlabel="長半径", ylabel="高さの導関数")
hline!([0])
vline!([1])

# 導関数が 0 になるときの長半径
ans_a2 = solve(g, a2)[1]
ans_a2 |> println

    1

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="")
plot(res[1]*res[2], xlims=(0.8, 1.2), ylims=(3.9, 4.2), xlabel="長半径", ylabel="高さ")
hline!([4])
vline!([1])

# 高さの最小値
h(a2 => 1) |> println

    4

高さは楕円の長半径 a2 が 1，すなわち 楕円の長径が 2 のとき最小値 4 になる。

function draw(a2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (s3, a, x1, x2) = (2*sqrt(4*a2^2 - 2*a2 + 1)/(2*a2 - 1), 2*a2^2/sqrt(4*a2^2 - 2*a2 + 1), sqrt(4*a2^2 - 2*a2 + 1)/(2*(2*a2^2 - a2 + 1)), 4*a2^2*sqrt(4*a2^2 - 2*a2 + 1)/(8*a2^2 - 2*a2 + 1))
    b2 = 1/2
    (a1, b1) = (b2, a2)
    y1 = (a - x1)*s3
    y2 = (a - x2)*s3
    h = s3*a
    @printf("上楕円の長径 = %g;  上楕円の短径 = %g\n", 2a1, 2b1)
    @printf("下楕円の長径 = %g;  下楕円の短径 = %g\n", 2a2, 2b2)
    @printf("s3 = %g;  a = %g;  h = %g;  a1 = %g;  b1 = %g;  x1 = %g;  y1 = %g;  a2 = %g;  b2 = %g;  x2 = %g;  y2 = %g\n", s3, a, h, a1, b1, x1, y1, a2, b2, x2, y2)
    plot([a, 0, -a, a], [0, a*s3, 0, 0], color=:magenta, lw=0.5)
    ellipse(0, b2, a2, b2, color=:blue)
    ellipse(0, b1 + 2b2, a1, b1, color=:red)
    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, b1 + 2b2, "上楕円:a1,b1\n(0,2b2+b1)", :red, :center, delta=-delta)
        point(0, b2, "下楕円:a2,b2,(0,b2)", :blue, :center, delta=-delta)
        point(x1, y1, " (x1,y1)", :green, :left, :vcenter)
        point(x2, y2, " (x2,y2)", :green, :left, :vcenter)
    end
end;
draw(1, true)