三角形に内接する楕円 -- 算法助術の公式97
1. 不等辺三角形の場合
不等辺三角形の場合について一般解が示されている。
和算の慣例として，楕円の長軸（長径），短軸（短径）を p, q としているが，ここでは長半径，短半径を使用することにする。
include("julia-source.txt");
@syms a, b, c, h, p, q
eq97 = -(b^2 - c^2)^2*h*q^2 + (b^2 - c^2)^2*q^3 + a^4*h*(2h - q)^2 - a^4*q*(2h - q)^2 - a^2*h*p^2*(2h - q)^2
eq97 = eq97(q => 2q, p => 2p)
eq_97 = eq97/4 |> simplify
eq_97 |> println
a^4*h*(h - q)^2 - 2*a^4*q*(h - q)^2 - 4*a^2*h*p^2*(h - q)^2 - h*q^2*(b^2 - c^2)^2 + 2*q^3*(b^2 - c^2)^2
2. 二等辺三角形の場合
b = c の場合は，簡単になる。
eq97_2 = eq97(c => b) |> factor
eq97_2 |> println
-4*a^2*(-h + q)^2*(-a^2*h + 2*a^2*q + 4*h*p^2)
更に a ≠ 0, h ≠ q なので非常に簡単な式になる。
eq97_3 = eq97_2/(4a^2*(q - h)^2)
eq97_3 |> println
a^2*h - 2*a^2*q - 4*h*p^2
2.1. a, h, q を与えて p を知るとき
res_p = solve(eq97_3, p)[2]
res_p |> println
a*sqrt( (h - 2*q)/h)/2
a = 5; h = 4; q = 1 のとき
ans_p = res_p(a => 5, h => 4, q => 1)
ans_p |> println
ans_p.evalf() |> println
5*sqrt(2)/4
1.76776695296637
2.2. a, h, p を与えて q を知るとき
res_q = solve(eq97_3, q)[1]
res_q |> println
h/2 - 2*h*p^2/a^2
a = 5; h = 4; q = 2 のとき
ans_q = res_q(a => 5, h => 4, p => 2)
ans_q |> println
ans_q.evalf() |> println
18/25
0.720000000000000
3. 関数として定義する
function formula97(a, h; p=0, q = 0)
p == 0 && q == 0 && return "error"
p != 0 && return h/2 - 2h*p^2/a^2
q != 0 && return a*sqrt( (h - 2q)/h)/2
end;
formula97(5, 4, p=1)
1.68
formula97(5, 4, q=1)
1.7677669529663689
formula97(5, 4)
"error"
4. 図を描いて確かめる
function draw(a, h, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
plot([0, a, a/2, 0], [0, 0, h, 0], color=:blue, lw=0.5)
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
q = 1
p = formula97(a, h, q=1)
@printf("p = %g; q = %g\n", p, q)
point(a/2, q, @sprintf("長半径 = %g\n短半径 = %g", p, q), :green, :center, delta=-delta/2)
ellipse(a/2, q, p, q, color=:green)
p = 1
q = formula97(a, h, p=1)
@printf("p = %g; q = %g\n", p, q)
point(a/2, q, @sprintf("長半径 = %g\n短半径 = %g", p, q), :red, :center, delta=-delta/2)
ellipse(a/2, q, p, q, 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)
end
end;
