算額(その407)
埼玉県神川町 光明寺 大正3年(1914)
https://toyo.repo.nii.ac.jp/?action=repository_action_common_download&item_id=2534&item_no=1&attribute_id=18&file_no=1
三角形内に界斜を引き,2 区分された領域に等円を入れる。大斜,中斜,小斜それぞれが 315 寸,257 寸,68 寸のとき,界斜はいかほどか。
数値が与えれているが,算額(その406)と同じである。
大斜,中斜,小斜,界斜を変数名として,
等円の半径と中心座標を r, (x1, r), (x2, r)
三角形の頂点座標を (x, y)
界斜と大斜の交点座標を (a, 0) として,以下の連立方程式を立て,nlsolve() により数値解を求める。
include("julia-source.txt");
using SymPy
@syms r::positive, a::positive, x::positive, y::positive, x1::positive,
x2::positive, 大斜::positive, 中斜::positive, 小斜::positive, 界斜::positive;
eq1 = (界斜 + 中斜 + a)r + (界斜 + 小斜 + 大斜 - a)r - 大斜*y
eq2 = x^2 + y^2 - 中斜^2
eq3 = (大斜 - x)^2 + y^2 - 小斜^2
eq4 = distance(a, 0, x, y, x1, r) - r^2
eq5 = distance(a, 0, x, y, x2, r) - r^2
eq6 = distance(0, 0, x, y, x1, r) - r^2
eq7 = distance(大斜, 0, x, y, x2, r) - r^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7])
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)
(界斜, r, a, x, y, x1, x2) = u
return [
r*(a + 中斜 + 界斜) + r*(-a + 大斜 + 小斜 + 界斜) - y*大斜, # eq1
x^2 + y^2 - 中斜^2, # eq2
y^2 - 小斜^2 + (-x + 大斜)^2, # eq3
-r^2 + (r - y*(a^2 - a*x - a*x1 + r*y + x*x1)/(a^2 - 2*a*x + x^2 + y^2))^2 + (x1 - (x1*(a^2 - 2*a*x + x^2 + y^2) - y*(a*r - a*y - r*x + x1*y))/(a^2 - 2*a*x + x^2 + y^2))^2, # eq4
-r^2 + (r - y*(a^2 - a*x - a*x2 + r*y + x*x2)/(a^2 - 2*a*x + x^2 + y^2))^2 + (x2 - (x2*(a^2 - 2*a*x + x^2 + y^2) - y*(a*r - a*y - r*x + x2*y))/(a^2 - 2*a*x + x^2 + y^2))^2, # eq5
-r^2 + (r - y*(r*y + x*x1)/(x^2 + y^2))^2 + (-x*(r*y + x*x1)/(x^2 + y^2) + x1)^2, # eq6
-r^2 + (r - y*(r*y + x*x2 - x*大斜 - x2*大斜 + 大斜^2)/(x^2 - 2*x*大斜 + y^2 + 大斜^2))^2 + (x2 - (x2*(x^2 - 2*x*大斜 + y^2 + 大斜^2) + y*(r*x - r*大斜 - x2*y + y*大斜))/(x^2 - 2*x*大斜 + y^2 + 大斜^2))^2, # eq7
]
end;
(大斜, 中斜, 小斜) = (315, 257, 68)
iniv = [big"35.0", 15, 220, 250, 30, 220, 250]
res = nls(H, ini=iniv);
println(res);
(BigFloat[40.00000000000000000000000000000000000000000000000000000000000000000000012868471, 13.99999999999999999999999999999999999999999999999999999999999999999999999479218, 231.0000000000000000000000000000000000000000000000000000000000000000000001855419, 255.0, 31.99999999999999999999999999999999999999999999999999999999999999999999999999945, 223.9999999999999999999999999999999999999999999999999999999999999999999999166771, 259.0000000000000000000000000000000000000000000000000000000000000000000000208086], true)
大斜, 中斜, 小斜 がそれぞれ 4315, 257, 68 寸のとき,
界斜 = 40; r = 14; a = 231; x = 255; y = 32; x1 = 224; x2 = 259である。
算額の術では 界斜 = sqrt((中斜 + 小斜)^2 - 大斜^2)/2 となっており,正しい答えが求まっていることがわかる。
(大斜, 中斜, 小斜) = (315, 257, 68)
sqrt((中斜 + 小斜)^2 - 大斜^2)/2
40.0
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(大斜, 中斜, 小斜) = (315, 257, 68)
(界斜, r, a, x, y, x1, x2) = res[1]
#(界斜, r, a, x, y, x1, x2) = [40, 50, 150, 150, 200, 103, 210]
@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;
※コメント投稿者のブログIDはブログ作成者のみに通知されます