算額(その565)
群馬の算額 51-3 稲荷神社 文政11年
http://takasakiwasan.web.fc2.com/gunnsann/g051-3.html
直角三角形内に全円と正方形が内接している。この直角三角形の弦を一辺とし,それと直角をなす辺の長さを前述の正方形の一辺の長さとする直角三角形を作る。この直角三角形に内接する円の直径が 5 寸,全円の直径が 7 寸のとき,もとの直角三角形の三辺(鈎,股,弦)の和を求めよ。
元の直角三角形の三辺を「鈎」,「股」,「弦」とする。
全円の半径と中心座標を r1, (r1, r1)
正方形の一辺の長さを a
上の直角三角形に内接する円の半径と中心座標を r2, (x2, y2),上の頂点を (x, y)
として,以下の連立方程式の数値解を求める。
include("julia-source.txt");
using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
鈎2::positive, 股2::positive, 弦2::positive,
r1::positive, r2::positive, a::positive,
x::positive, y::positive,
x2::positive, y2::positive
鈎2 = a # sqrt(x^2 +(y - 鈎)^2)
x = sqrt(a^2 - (y - 鈎)^2)
弦2 = sqrt((股 - x)^2 + y^2)
股2 = 弦
eq1 = 鈎^2 + 股^2 - 弦^2
eq2 = 鈎 + 股 - 弦 - 2r1
eq3 = (鈎 - a)/a - 鈎/股
eq4 = 鈎2^2 + 股2^2 - 弦2^2
eq5 = 鈎2 + 股2 - 弦2 - 2r2
eq6 = distance(0, 鈎, 股, 0, x2, y2) - r2^2 |> numerator;
eq7 = distance(x, y, 股, 0, x2, y2) - r2^2 |> numerator;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (鈎, 股, 弦, a, y, x2, y2))
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=big"1e-40")
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
v = r.zero
end
return v, r.f_converged
end;
function H(u)
(鈎, 股, 弦, a, y, x2, y2) = u
return [
-弦^2 + 股^2 + 鈎^2, # eq1
-2*r1 - 弦 + 股 + 鈎, # eq2
-鈎/股 + (-a + 鈎)/a, # eq3
a^2 - y^2 + 弦^2 - (股 - sqrt(a^2 - (y - 鈎)^2))^2, # eq4
a - 2*r2 + 弦 - sqrt(y^2 + (股 - sqrt(a^2 - (y - 鈎)^2))^2), # eq5
-r2^2 + (x2 - 股*(x2*股 - y2*鈎 + 鈎^2)/(股^2 + 鈎^2))^2 + (y2 - 鈎*(-x2*股 + y2*鈎 + 股^2)/(股^2 + 鈎^2))^2, # eq6
-r2^2 + (x2 - (x2*(a^2 + 2*y*鈎 + 股^2 - 2*股*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2) - 鈎^2) - y*(x2*y - y*股 + y2*股 - y2*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2)))/(a^2 + 2*y*鈎 + 股^2 - 2*股*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2) - 鈎^2))^2 + (-y*(-x2*股 + x2*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2) + y*y2 + 股^2 - 股*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2))/(a^2 + 2*y*鈎 + 股^2 - 2*股*sqrt(a^2 - y^2 + 2*y*鈎 - 鈎^2) - 鈎^2) + y2)^2, # eq7
]
end;
(r1, r2) = (7, 5) .// 2
iniv = BigFloat[9.41379, 16.89655, 19.7931, 6.27586, 15.2069, 3.45, 10.55]
res = nls(H, ini=iniv)
(BigFloat[10.49999999999999999999999999999999999999999129919479365547439836695548116931581, 14.00000000000000000000000000000000000000001100960233399594221790231656113690163, 17.5000000000000000000000000000000000000000023087971276514166162692720423047526, 5.99999999999999999999999999999999999999999995377737351386752851972765331489346, 15.29999999999999999999999999999999999999999475386066794051341170073747364389192, 3.500000000000000000000000000000000000000000057743557772276350570950054497150227, 10.99999999999999999999999999999999999999999392286748897320175788883023235161591], true)
元の直角三角形の鈎,股,弦の和は 42 寸である。
その他のパラメータは以下の通り。
鈎 = 10.5; 股 = 14; 弦 = 17.5; a = 6; x = 3.6; y = 15.3; , x2 = 3.5; y2 =11
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2) = (7, 5) .// 2
(鈎, 股, 弦, a, y, x2, y2) = Float64.(res[1])
x = sqrt(a^2 - (y -鈎)^2)
@printf("鈎 = %g; 股 = %g; 弦 = %g; a = %g; x = %g; y = %g; , x2 = %g; y2 =%g\n", 鈎, 股, 弦, a, x, y, x2, y2)
@printf("元の直角三角形の鈎,股,弦の和 = %g\n", 鈎 + 股 + 弦)
plot([0, 0, 股, x, 0, 股], [鈎, 0, 0, y, 鈎, 0], color=:blue, lw=0.5)
plot!([0, a, a, 0, 0], [0, 0, a, a, 0], color=:magenta, lw=0.5)
circle(r1, r1, r1)
circle(x2, y2, r2, :green)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
# hline!([0], color=:black, lw=0.5)
# vline!([0], color=:black, lw=0.5)
point(x, y, " (x,y)", :blue, :left, :vcenter)
point(0, 鈎, " 鈎", :blue, :left, :vcenter, delta=delta/2)
point(股, 0, "股", :blue, :left, :bottom, delta=delta/2)
point(a, 0, " a", :magenta, :left, :bottom, delta=delta/2)
point(0, a, " a", :magenta, :left, :bottom, delta=delta/2)
point(r1, r1, "全円:r1,(r1,r1)", :red, :center, :top, delta=-delta/2)
point(x2, y2, "円:r2,(x2,y2)", :green, :center, :top, delta=-delta/2)
end
end;