和算の心（その010）
等円術（II）：三斜等円術（一般公式）
米山忠興：等円術（II）：三斜等円術（一般公式），東洋大学紀要 自然科学編，第 50 号，P.107-121, 2006/03
https://toyo.repo.nii.ac.jp/record/2534/files/shizenkagakuhen50_107-121_OCR.pdf
等円が 2 個の場合
三角形内に界斜を引き，2 区分された領域に等円を入れる。大斜 = 9/2，中斜= 4，小斜 = 4 が与えられたとき，界斜はいかほどか。
結論
以下の公式を適宜使用する。
公式1: 界斜 = sqrt((中斜 + 小斜)^2 - 大斜^2)/2
公式2: 大斜 = sqrt((中斜 + 小斜)^2 - 4界斜^2)
公式3: 中斜 = sqrt(大斜^2 + 4界斜^2) - 小斜
公式4: 小斜 = sqrt(大斜^2 + 4界斜^2) - 中斜
公式5: 等円半径 = sqrt((大斜 + 中斜 + 小斜)(大斜 - 中斜 + 小斜)(中斜 - 大斜 + 小斜)(大斜 + 中斜 - 小斜))/(2(中斜 + 大斜 + 小斜 + 2*界斜))
公式6-1: x = (中斜^2 + 大斜^2 - 小斜^2)/(2大斜)
公式6-2: y = sqrt((中斜 + 大斜 + 小斜)(大斜 - 中斜 + 小斜)(中斜 - 大斜 + 小斜)(大斜 + 中斜 - 小斜))/(2大斜)
公式7-1: 大斜*(小斜 + 界斜)/(中斜 + 小斜 + 2*界斜)
公式7-2: 大斜*(中斜 + 界斜)/(中斜 + 小斜 + 2*界斜)
公式8-1: r*(x + sqrt(x^2 + y^2))/y
公式8-2: (-a*r + a*y + r*x + r*sqrt(a^2 - 2*a*x + x^2 + y^2))/y
公式を導く手順
大斜，中斜，小斜，界斜を変数名として，
等円の半径と中心座標を r, (x1, r), (x2, r)
三角形の頂点座標を (x, y)
界斜と大斜の交点座標を (a, 0) として，以下の手順により界斜のほか，等円の直径，等円の中心座標などを求める。
まず，大斜，中斜，小斜，界斜の相互関係を表す恒等式を導く。
求めるパラメータ以外の未知パラメータ a, r を消去する。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r::positive, a::positive, x::positive, y::positive, x1::positive,
x2::positive, 大斜::positive, 中斜::positive, 小斜::positive, 界斜::positive;
eq1 = r*(小斜 + a + 界斜) - a*y
eq2 = r*(界斜 + (大斜 - a) + 中斜)- (大斜 - a)*y;
eq1, eq2 を辺々割って eq3 とする。
# eq3 = eq1.lhs/eq2.lhs - eq1.rhs/eq2.rhs # こうしてもよいが
eq3 = (a + 小斜 + 界斜)/(中斜 + 大斜 + 界斜 - a) - a/(大斜 - a);
ans_a = solve(eq3, a)[1];
ans_a |> println
ans_a(大斜 => 4.5, 中斜 => 4, 小斜 => 3, 界斜 => 2.68095132369090) |> println
大斜*(小斜 + 界斜)/(中斜 + 小斜 + 2*界斜)
2.06798918304242
注：ここで計算される a は大斜の左端点からの距離（左端点を原点とした）距離なので，中斜が左側にあるときには図を描くときに注意すること。
eq8 = y^2 + x^2 - 小斜^2
eq9 = y^2 + (a - x)^2 - 界斜^2
eq10 = y^2 + (大斜 - x)^2 - 中斜^2
eq10 |> println
y^2 - 中斜^2 + (-x + 大斜)^2
eq11 = eq9 - eq8 |> expand
eq11 |> println
a^2 - 2*a*x + 小斜^2 - 界斜^2
eq12 = eq10 - eq8 |> expand
eq12 |> println
-2*x*大斜 - 中斜^2 + 大斜^2 + 小斜^2
eq1112 = 2a*x/(2大斜*x) - (小斜^2 + a^2 - 界斜^2)/(大斜^2 - 中斜^2 + 小斜^2)
eq1112 = eq1112(a => ans_a) |> simplify
eq1112 |> println
(-大斜^2*(小斜 + 界斜)^2 + (小斜 + 界斜)*(中斜 + 小斜 + 2*界斜)*(-中斜^2 + 大斜^2 + 小斜^2) + (-小斜^2 + 界斜^2)*(中斜 + 小斜 + 2*界斜)^2)/((中斜 + 小斜 + 2*界斜)^2*(-中斜^2 + 大斜^2 + 小斜^2))
大斜，中斜，小斜，界斜の 4 変数を含む式は
(-大斜^2*(小斜 + 界斜)^2 + (小斜 + 界斜)*(中斜 + 小斜 + 2*界斜)*(-中斜^2 + 大斜^2 + 小斜^2) + (-小斜^2 + 界斜^2)*(中斜 + 小斜 + 2*界斜)^2)/((中斜 + 小斜 + 2*界斜)^2*(-中斜^2 + 大斜^2 + 小斜^2)) = 0
となる。
この式のうち，3 個に実際の値を与えて 4 番目の変数が取るべき値（解）を求めることができる。
1. 大斜，中斜，小斜をあたえて，界斜を求める
ans_界斜 = solve(eq1112, 界斜)[1]
ans_界斜 |> println
ans_界斜(大斜 => 4.5, 中斜 => 4, 小斜 => 3).evalf() |> println
sqrt(中斜^2 + 2*中斜*小斜 - 大斜^2 + 小斜^2)/2
2.68095132369090
公式1: 界斜 = sqrt((中斜 + 小斜)^2 - 大斜^2)/2
2. 中斜，小斜，界斜を与えて，大斜を求める
ans_大斜 = solve(eq1112, 大斜)[1]
ans_大斜 |> println
ans_大斜(中斜 => 4, 小斜 => 3, 界斜 => 2.68095132369090) |> println
sqrt(中斜^2 + 2*中斜*小斜 + 小斜^2 - 4*界斜^2)
4.50000000000000
公式2: 大斜 = sqrt((中斜 + 小斜)^2 - 4界斜^2)
3. 大斜，小斜，界斜を与えて，中斜を求める
ans_中斜 = solve(eq1112, 中斜)[1]
ans_中斜 |> println
ans_中斜(大斜 => 4.5, 小斜 => 3, 界斜 => 2.68095132369090) |> println
-小斜 + sqrt(大斜^2 + 4*界斜^2)
4.00000000000000
公式3: 中斜 = sqrt(大斜^2 + 4界斜^2) - 小斜
4. 大斜，中斜，界斜を与えて，小斜を求める
ans_小斜 = solve(eq1112, 小斜)[1]
ans_小斜 |> println
ans_小斜(大斜 => 4.5, 中斜 => 4, 界斜 => 2.68095132369090) |> println
-中斜 + sqrt(大斜^2 + 4*界斜^2)
3.00000000000000
公式4: 小斜 = sqrt(大斜^2 + 4界斜^2) - 中斜
5. 等円の半径を求める
eq1, eq2 による三角形の面積とヘロンの公式による面積が等しいとする連立方程式から，等円の半径 r を求めることができる。
@syms r
s = (大斜 + 中斜 + 小斜)/2
eq21 = (r*(小斜 + a + 界斜) + r*(界斜 + (大斜 - a) + 中斜)) - (大斜 - a)*y
eq22 = sqrt(s*(s - 大斜)*(s - 中斜)*(s - 小斜)) - (大斜 - a)*y/2;
res = solve([eq21, eq22], (r, y));
res[r] |> println
res[r](大斜 => 4.5, 中斜 => 4, 小斜 => 3, 界斜 => 2.68095132369090) |> println
sqrt(中斜 + 大斜 + 小斜)*sqrt(-中斜^3 + 中斜^2*大斜 + 中斜^2*小斜 + 中斜*大斜^2 - 2*中斜*大斜*小斜 + 中斜*小斜^2 - 大斜^3 + 大斜^2*小斜 + 大斜*小斜^2 - 小斜^3)/(2*中斜 + 2*大斜 + 2*小斜 + 4*界斜)
0.697585939193294
公式5: 等円半径 = sqrt((大斜 + 中斜 + 小斜)*(大斜 - 中斜 + 小斜)*(中斜 - 大斜 + 小斜)*(大斜 + 中斜 - 小斜))/(2*(中斜 + 大斜 + 小斜 + 2*界斜))
6. 三角形の頂点座標 (x, y) を求める
x, y は，以下の連立方程式で求めることができる。
eq21 = x^2 + y^2 - 中斜^2
eq22 = (大斜 - x)^2 + y^2 - 小斜^2
res2 = solve([eq21, eq22], (x, y))[1];
# x
res2[1] |> println
res2[1](大斜 => 4.5, 中斜 => 4, 小斜 => 3) |> println
(中斜^2 + 大斜^2 - 小斜^2)/(2*大斜)
3.02777777777778
公式6-1: x = (中斜^2 + 大斜^2 - 小斜^2)/(2大斜)
# y
res2[2] |> println
res2[2](大斜 => 4.5, 中斜 => 4, 小斜 => 3) |> println
sqrt(-(中斜 - 大斜 - 小斜)*(中斜 - 大斜 + 小斜)*(中斜 + 大斜 - 小斜))*sqrt(中斜 + 大斜 + 小斜)/(2*大斜)
2.61391693219105
公式6-2: y = sqrt((中斜 + 大斜 + 小斜)(大斜 - 中斜 + 小斜)*(中斜 - 大斜 + 小斜)*(大斜 + 中斜 - 小斜))/(2大斜)
7. a を求める 界斜と大斜の交点座標(a, 0)
小斜が左側にある場合
ans_a |> println
ans_a(大斜 => 4.5, 中斜 => 4, 小斜 => 3, 界斜 => 2.68095132369090) |> println
大斜*(小斜 + 界斜)/(中斜 + 小斜 + 2*界斜)
2.06798918304242
公式7-1: 大斜*(小斜 + 界斜)/(中斜 + 小斜 + 2*界斜)
中斜が左側にある場合には ans_a = 大斜 - 大斜*(小斜 + 界斜)/(中斜 + 小斜 + 2*界斜) = 大斜*(中斜 + 界斜)/(中斜 + 小斜 + 2*界斜)
ans_a = 大斜 - 大斜*(小斜 + 界斜)/(中斜 + 小斜 + 2*界斜) |> simplify
ans_a |> println
ans_a(大斜 => 4.5, 中斜 => 4, 小斜 => 3, 界斜 => 2.68095132369090) |> println
大斜*(中斜 + 界斜)/(中斜 + 小斜 + 2*界斜)
2.43201081695758
公式7-2: 大斜*(中斜 + 界斜)/(中斜 + 小斜 + 2*界斜)
8. 等円の中心座標 (x1, r), (x2, r)
eq41 = dist2(0, 0, x, y, x1, r, r)
eq42 = dist2(a, 0, x, y, x2, r, r)
res3 = solve([eq41, eq42], (x1, x2))[4] # 4 of 4
(r*(x/y + sqrt(x^2 + y^2)/y), r*sqrt(a^2 - 2*a*x + x^2 + y^2)/y + (-a*r + a*y + r*x)/y)
公式8-1: r*(x + sqrt(x^2 + y^2))/y
公式8-2: (-a*r + a*y + r*x + r*sqrt(a^2 - 2*a*x + x^2 + y^2))/y
res3[1](a => 2.43201081695758, x => 3.0277777777777777, y => 2.613916932191048, r => 0.6975859391932936) |> println
res3[2](a => 2.43201081695758, x => 3.0277777777777777, y => 2.613916932191048, r => 0.6975859391932936) |> println
1.87552974663334
3.30648107032424
界斜f(大斜, 中斜, 小斜) = sqrt((中斜 + 小斜)^2 - 大斜^2)/2
大斜f(中斜, 小斜, 界斜) = sqrt((中斜 + 小斜)^2 - 4界斜^2)
中斜f(大斜, 小斜, 界斜) = sqrt(大斜^2 + 4界斜^2) - 小斜
小斜f(大斜, 中斜, 界斜) = sqrt(大斜^2 + 4界斜^2) - 中斜
等円半径f(大斜, 中斜, 小斜, 界斜) = sqrt((大斜 + 中斜 + 小斜)*(大斜 - 中斜 + 小斜)*(中斜 - 大斜 + 小斜)*(大斜 + 中斜 - 小斜))/(2*(中斜 + 大斜 + 小斜 + 2*界斜))
xf(大斜, 中斜, 小斜) = (中斜^2 + 大斜^2 - 小斜^2)/(2大斜)
yf(大斜, 中斜, 小斜) = sqrt((中斜 + 大斜 + 小斜)*(大斜 - 中斜 + 小斜)*(中斜 - 大斜 + 小斜)*(大斜 + 中斜 - 小斜))/(2大斜)
af1(大斜, 中斜, 小斜, 界斜) = 大斜*(小斜 + 界斜)/(中斜 + 小斜 + 2*界斜)
af2(大斜, 中斜, 小斜, 界斜) = 大斜*(中斜 + 界斜)/(中斜 + 小斜 + 2*界斜)
x1f(a, x, y, r) = r*(x/y + sqrt(x^2 + y^2)/y)
x2f(a, x, y, r) = r*sqrt(a^2 - 2*a*x + x^2 + y^2)/y + (-a*r + a*y + r*x)/y;
界斜f(4.5, 4, 3) |> println
大斜f(4, 3, 2.68095132369090) |> println
中斜f(4.5, 3, 2.68095132369090) |> println
小斜f(4.5, 4, 2.68095132369090) |> println
等円半径f(4.5, 4, 3, 2.68095132369090) |> println
xf(4.5, 4, 3) |> println
yf(4.5, 4, 3) |> println
af1(4.5, 4, 3, 2.68095132369090) |> println
af2(4.5, 4, 3, 2.68095132369090) |> println
x1f(2.43201081695758, 3.0277777777777777, 2.613916932191048, 0.6975859391932936) |> println
x2f(2.43201081695758, 3.0277777777777777, 2.613916932191048, 0.6975859391932936) |> println
2.680951323690902
4.500000000000004
3.9999999999999973
2.9999999999999973
0.6975859391932936
3.0277777777777777
2.613916932191048
2.0679891830424224
2.432010816957577
1.875529746633338
3.3064810703242418
大斜, 中斜, 小斜 がそれぞれ 4.5, 4, 3 のとき，
界斜 = 2.68095; r = 0.697586; a = 2.43201; x = 3.02778; y = 2.61392; x1 = 1.87553; x2 = 3.30648
である。
算額の術では 界斜 = sqrt((中斜 + 小斜)^2 - 大斜^2)/2 となっており，正しい答えが求まっていることがわかる。
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(大斜, 中斜, 小斜) = (9//2, 4, 3)
界斜 = 界斜f(4.5, 4, 3)
r = 等円半径f(4.5, 4, 3, 界斜)
x = xf(大斜, 中斜, 小斜)
y = yf(大斜, 中斜, 小斜)
a = af2(大斜, 中斜, 小斜, 界斜)
x1 = x1f(a, x, y, r)
x2 = x2f(a, x, y, r)
@printf("界斜 = %g; r = %g; a = %g; x = %g; y = %g; x1 = %g; x2 = %g\n", 界斜, r, a, x, y, x1, x2)
plot([0, 大斜, x, 0], [0, 0, y, 0], color=:black, lw=0.5)
circle(x1, r, r)
circle(x2, r, r)
segment(a, 0, x, y, :blue)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(大斜, 0, " 大斜", :black, :left, :bottom, delta=delta)
point(x, y, " (x,y)", :black, :left, :vcenter, delta=-delta)
point(x1, r, "(x1,r)", :red, :center, delta=-delta)
point(x2, r, "(x2,r)", :red, :center, delta=-delta)
point(a, 0, " a", :black, :left, :bottom, delta=delta)
point(x/2, 2y/3, "中斜", :black, :left, mark=false)
point((大斜 + x)/2, 2y/3, "小斜", :black, :left, mark=false)
point((a + x)/2, 2y/3, " 界斜", :blue, :left, mark=false)
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
plot!(xlims=(0, 5))
else
plot!(showaxis=false)
end
end;
draw(true)
