算額（その116）
一関市博物館 > 和算に挑戦 > 令和4年度出題問題(3) [上級問題] (高校生・一般向き)
文政13年(1830)刊『算法新書』より
https://www.city.ichinoseki.iwate.jp/museum/wasan/r4/hard.html
三角形を線分（界斜）で分けた2つの三角形それぞれに内接する半径の等しい2つの円がある。三角形の辺の大きい順に 15 寸，11 寸，6 寸としたとき，界斜の長さはいかほどか。
下図のように記号を定め，方程式を解く。方針としては 2 つの円それぞれの中心から接線までの距離に関する条件と，三角形の面積 20√2 について「底辺×高さ/2」とヘロンの公式による面積に関しての方程式である。
例によって，方程式が複雑になると solve() では解けなくなってしまうので nlsolve() で数値解を得ることにする。
using SymPy
function distance(x1, y1, x2, y2, x0, y0)
p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
l = sympy.Line(p1, p2)
l.distance(sympy.Point(x0, y0))^2 # 二乗距離を返す！！
end;
@syms x1, x2, x3, x4, r
h = 8*sqrt(Sym(2))/3
eq1 = distance(0, 0, x1, h, x3, r) - r^2
eq2 = distance(x2, 0, x1, h, x3, r) - r^2
eq3 = distance(x2, 0, x1, h, x4, r) - r^2
eq4 = distance(x1, h, 15, 0, x4, r) - r^2
a = sqrt(h^2 + (x1 - x2)^2)
eq5 = r * (16 + a) - 20sqrt(Sym(2));
# solve([eq1, eq2, eq3, eq4, eq5])
println(eq1, ",")
println(eq2, ",")
println(eq3, ",")
println(eq4, ",")
println(eq5, ",")
-r^2 + (r - 8*(16*r + 3*sqrt(2)*x1*x3)/(9*x1^2 + 128))^2 + (-3*x1*(8*sqrt(2)*r + 3*x1*x3)/(9*x1^2 + 128) + x3)^2,
-r^2 + (r - sqrt(2)*(64*sqrt(2)*r - 24*x1*x2 + 24*x1*x3 + 24*x2^2 - 24*x2*x3)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2 + (x3 - (24*sqrt(2)*r*x1 - 24*sqrt(2)*r*x2 + 9*x1^2*x3 - 18*x1*x2*x3 + 9*x2^2*x3 + 128*x2)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2,
-r^2 + (r - sqrt(2)*(64*sqrt(2)*r - 24*x1*x2 + 24*x1*x4 + 24*x2^2 - 24*x2*x4)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2 + (x4 - (24*sqrt(2)*r*x1 - 24*sqrt(2)*r*x2 + 9*x1^2*x4 - 18*x1*x2*x4 + 9*x2^2*x4 + 128*x2)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2,
-r^2 + (r - sqrt(2)*(64*sqrt(2)*r + 24*x1*x4 - 360*x1 - 360*x4 + 5400)/(9*x1^2 - 270*x1 + 2153))^2 + (x4 - 3*(8*sqrt(2)*r*x1 - 120*sqrt(2)*r + 3*x1^2*x4 - 90*x1*x4 + 675*x4 + 640)/(9*x1^2 - 270*x1 + 2153))^2,
r*(sqrt((x1 - x2)^2 + 128/9) + 16) - 20*sqrt(2),
using NLsolve
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=1e-14)
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
v = r.zero
end
return v, r.f_converged
end;
function H(u)
(x1, x2, x3, x4, r) = u
return [
-r^2 + (r - 8*(16*r + 3*sqrt(2)*x1*x3)/(9*x1^2 + 128))^2 + (-3*x1*(8*sqrt(2)*r + 3*x1*x3)/(9*x1^2 + 128) + x3)^2,
-r^2 + (r - 8*(16*r - 3*sqrt(2)*x1*x2 + 3*sqrt(2)*x1*x3 + 3*sqrt(2)*x2^2 - 3*sqrt(2)*x2*x3)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2 + (x3 - (24*sqrt(2)*r*x1 - 24*sqrt(2)*r*x2 + 9*x1^2*x3 - 18*x1*x2*x3 + 9*x2^2*x3 + 128*x2)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2,
-r^2 + (r - sqrt(2)*(64*sqrt(2)*r - 24*x1*x2 + 24*x1*x4 + 24*x2^2 - 24*x2*x4)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2 + (x4 - (24*sqrt(2)*r*x1 - 24*sqrt(2)*r*x2 + 9*x1^2*x4 - 18*x1*x2*x4 + 9*x2^2*x4 + 128*x2)/(9*x1^2 - 18*x1*x2 + 9*x2^2 + 128))^2,
-r^2 + (r - sqrt(2)*(64*sqrt(2)*r + 24*x1*x4 - 360*x1 - 360*x4 + 5400)/(9*x1^2 - 270*x1 + 2153))^2 + (x4 - 3*(8*sqrt(2)*r*x1 - 120*sqrt(2)*r + 3*x1^2*x4 - 90*x1*x4 + 675*x4 + 640)/(9*x1^2 - 270*x1 + 2153))^2,
r*(sqrt((x1 - x2)^2 + 128/9) + 16) - 20*sqrt(2),
]
end;
iniv = [10.8, 9.5, 8.3, 11.7, 1.25]
res = nls(H, ini=iniv)
println(res)
([10.333333333333236, 8.99999999999995, 7.999999999999937, 10.999999999999943, 1.4142135623730963], false)
多分，x1 = 10.333333333333236 = 31/3，r = 1.4142135623651715 = √2 であろう。
界斜の長さは
h = 8*sqrt(2)/3
a = sqrt(h^2 + (x1 - x2)^2) = 4 である。
using Plots
using Printf
function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
θ = beginangle:0.1:endangle
x = r.*cosd.(θ)
y = r.*sind.(θ)
plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;
function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
mark && scatter!([x], [y], color=color, markerstrokewidth=0)
annotate!(x, y, text(string, 10, position, color, vertical))
end;
function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(x1, x2, x3, x4, r) = res[1]
h = 8*sqrt(Sym(2))/3
a = sqrt(h^2 + (x1 - x2)^2)
plot([0, 15, x1, 0], [0, 0, h, 0])
circle(x3, r, r)
circle(x4, r, r)
segment(x2, 0, x1, h)
if more == true
@printf("x1 = %.3f, h = %.3f, x2 = %.3f, x3 = %.3f, x4 = %.3f, r = %.5f\n", x1, h, x2, x3, x4, r)
@printf("界斜の長さ = %.1f", a)
point(x1, h, " (x1,h)", :blue, :left, :bottom)
point(x2, 0, "x2", :blue)
point(x3, r, "(x3,r)", :red, :top)
point(x4, r, "(x4,r)", :red, :top)
hline!([0], color=:black, lw=0.5, ylims=(-1, 5))
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end
draw(true)
savefig("fig1.png");
x1 = 10.333, h = 3.771, x2 = 9.000, x3 = 8.000, x4 = 11.000, r = 1.41421
界斜の長さ = 4.0
