算額(その613)
高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/6
円弧内に斜線と等円 2 個を入れる。等円の直径が 4 寸で,弦の長さから「子 = 14 寸」を引いた丑の長さ」(図での xa - x)を求めよ。
円弧をなす円の半径と中心座標を R, (0, 0)
弦と y 軸の交点座標を (0, y) とする。
等円の半径と中心座標を r, (x, y + r), (x2, y2); x = 14 - sqrt(R^2 - y^2)
include("julia-source.txt");
using SymPy
@syms r::positive, R::positive, x::positive, y::positive,
x2::negative, y2::positive, y3::positive, xa::positive, xb::positive
r = 4//2
xa = sqrt(R^2 - y^2)
xb = sqrt(R^2 - y3^2)
eq1 = (14 - sqrt(R^2 - y^2))^2 + (y + r)^2 - (R - r)^2
eq2 = x2^2 + y2^2 - (R - r)^2
eq3 = (y3 - y)/(xa + xb) * y2/x2 + 1
eq4 = distance(-xa, y, xb, y3, 14 - sqrt(R^2 - y^2), y + r) - r^2
eq5 = (xb - xa)^2/4 + ((y3 - y)/2+y)^2 - (R - 2r)^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=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, y, x2, y2, y3) = u
return [
(14 - sqrt(R^2 - y^2))^2 - (R - 2)^2 + (y + 2)^2, # eq1
x2^2 + y2^2 - (R - 2)^2, # eq2
1 + y2*(-y + y3)/(x2*(sqrt(R^2 - y^2) + sqrt(R^2 - y3^2))), # eq3
(y - ((y + 2)*(R^2 - y*y3 + sqrt(R^2 - y^2)*sqrt(R^2 - y3^2)) - (sqrt(R^2 - y^2) + sqrt(R^2 - y3^2))*(7*y - 7*y3 + sqrt(R^2 - y^2) + sqrt(R^2 - y3^2)))/(R^2 - y*y3 + sqrt(R^2 - y^2)*sqrt(R^2 - y3^2)) + 2)^2 + (-sqrt(R^2 - y^2) - ((14 - sqrt(R^2 - y^2))*(R^2 - y*y3 + sqrt(R^2 - y^2)*sqrt(R^2 - y3^2)) - (y - y3)*(7*y - 7*y3 + sqrt(R^2 - y^2) + sqrt(R^2 - y3^2)))/(R^2 - y*y3 + sqrt(R^2 - y^2)*sqrt(R^2 - y3^2)) + 14)^2 - 4, # eq4
-(R - 4)^2 + (y/2 + y3/2)^2 + (-sqrt(R^2 - y^2) + sqrt(R^2 - y3^2))^2/4, # eq5
]
end;
iniv = BigFloat[8.8, 2.5, -1.67, 6.6, 6.7]
res = nls(H, ini=iniv)
(BigFloat[9.165816326530612244897959183673469387755102040816326530612244897959183673469382, 2.839183673469387755102040816326530612244897959183673469387755102040816326530621, -2.00642857142857142857142857142857142857142857142857137843826725205845906344671, 6.879183673469387755102040816326530612244897959183673492933533457683883923317655, 7.079183673469387755102040816326530612244897959183673469387755102040816326530546], true)
丑の長さは xa - x = 8.715 - 5.285 = 3.43 寸である。
R = 9.16582; x = 5.285; y = 2.83918; x2 = -2.00643; y2 = 6.87918; y3 = 7.07918; xa = 8.715; xb = 5.82214
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r = 2
(R, y, x2, y2, y3) = res[1]
xa = sqrt(R^2 - y^2)
xb = sqrt(R^2 - y3^2)
x = 14 - xa
@printf("R = %g; x = %g; y = %g; x2 = %g; y2 = %g; y3 = %g; xa = %g; xb = %g\n", R, x, y, x2, y2, y3, xa, xb)
@printf("丑 = %g\n", xa - x)
plot()
circle(0, 0, R, beginangle=0, endangle=180)
circle(x, y + r, r)
circle(x2, y2, r)
plot!([sqrt(R^2 - y3^2), -sqrt(R^2 - y^2), sqrt(R^2 - y^2)], [y3, y, y], color=:blue, lw=0.5)
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(0, y, " y", :blue, :left, delta=-delta/2)
point(-xa, y, "(-xa,y)", :blue)
point(xa, y, "(xa,y)", :blue, :right, delta=-delta/2)
point(x, y, "(x,y)", :blue, :center, delta=-delta/2)
point(sqrt(R^2 - y3^2), y3, " (xb,y3)", :blue, :left, :vcenter)
point(x, y + r, "(x,y+r)", :red, :center, delta=-delta/2)
point(x2, y2, "(x2,y2)", :red, :center, delta=-delta/2)
point(R, 0, "R ", :red, :right, :bottom, delta=delta/2)
segment(x, y, xa, y, lw=2)
point((x+xa)/2, y, "丑", :black, :center, :bottom, delta=delta, mark=false)
end
end;