算額(その108)
東京都府中市宮町 大國魂神社 明治18年
大國魂神社の算額(上段第2問)についての考察
http://www2.ttcn.ne.jp/~nagai/waseda/wasan/ok2.pdf
外円内に甲円,乙円,丙円が内接している。更に弦にも接している。甲円,乙円,丙円の径をそれぞれ 4 寸,13 寸,8 寸 4 分 5 厘のとき,外円の径を求めよ。
図のように記号を定め,方程式を立てる。
外円の半径を R とする。甲円,乙円,丙円の中心座標と半径を (0, 400 - R, 400),(x4, y4, 1300), (x, y, 845) とする。
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, x3, y3, x4, y4, x, y, R, r, x0, y0;
eq1 = x^2 + y^2 - (R - 845)^2
eq2 = distance(x1, (800 - R), x3, y3, x4, y4) - 1300^2
eq3 = (x - x0)^2 + (y - y0)^2 - 845^2
eq4 = distance(x1, (800 - R), x3, y3, x, y) - 845^2;
eq5 = R + y4 - 2100
eq6 = x4^2 + y4^2 - (R - 1300)^2
eq7 = x3^2 + y3^2 - R^2
eq8 = x1^2 + (R - 800)^2 - R^2
eq9 = x0 - (x1 + x3) / 2
eq10 = y0 - ((800 - R) + y3) / 2;
方程式は solve() では解けないので,nlsolve() で数値解を求める。
eq1 |> println
eq2 |> println
eq3 |> println
eq4 |> println
eq5 |> println
eq6 |> println
eq7 |> println
eq8 |> println
eq9 |> println
eq10 |> println
x^2 + y^2 - (R - 845)^2
(x4 - (x1*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) - (x1 - x3)*(R^2 + R*y3 + R*y4 - 1600*R + x1^2 - x1*x3 - x1*x4 + x3*x4 + y3*y4 - 800*y3 - 800*y4 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 + (y4 - ((800 - R)*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) + (R + y3 - 800)*(R^2 + R*y3 + R*y4 - 1600*R + x1^2 - x1*x3 - x1*x4 + x3*x4 + y3*y4 - 800*y3 - 800*y4 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 - 1690000
(x - x0)^2 + (y - y0)^2 - 714025
(x - (x1*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) - (x1 - x3)*(R^2 + R*y + R*y3 - 1600*R - x*x1 + x*x3 + x1^2 - x1*x3 + y*y3 - 800*y - 800*y3 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 + (y - ((800 - R)*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) + (R + y3 - 800)*(R^2 + R*y + R*y3 - 1600*R - x*x1 + x*x3 + x1^2 - x1*x3 + y*y3 - 800*y - 800*y3 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 - 714025
R + y4 - 2100
x4^2 + y4^2 - (R - 1300)^2
-R^2 + x3^2 + y3^2
-R^2 + x1^2 + (R - 800)^2
x0 - x1/2 - x3/2
R/2 + y0 - y3/2 - 400
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10])
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, x3, y3, x, y, R, x4, y4, x0, y0) = u
return [
x^2 + y^2 - (R - 845)^2,
(x4 - (x1*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) - (x1 - x3)*(R^2 + R*y3 + R*y4 - 1600*R + x1^2 - x1*x3 - x1*x4 + x3*x4 + y3*y4 - 800*y3 - 800*y4 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 + (y4 - ((800 - R)*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) + (R + y3 - 800)*(R^2 + R*y3 + R*y4 - 1600*R + x1^2 - x1*x3 - x1*x4 + x3*x4 + y3*y4 - 800*y3 - 800*y4 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 - 1690000,
(x - x0)^2 + (y - y0)^2 - 714025,
(x - (x1*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) - (x1 - x3)*(R^2 + R*y + R*y3 - 1600*R - x*x1 + x*x3 + x1^2 - x1*x3 + y*y3 - 800*y - 800*y3 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 + (y - ((800 - R)*(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000) + (R + y3 - 800)*(R^2 + R*y + R*y3 - 1600*R - x*x1 + x*x3 + x1^2 - x1*x3 + y*y3 - 800*y - 800*y3 + 640000))/(R^2 + 2*R*y3 - 1600*R + x1^2 - 2*x1*x3 + x3^2 + y3^2 - 1600*y3 + 640000))^2 - 714025,
R + y4 - 2100,
x4^2 + y4^2 - (R - 1300)^2,
-R^2 + x3^2 + y3^2,
-R^2 + x1^2 + (R - 800)^2,
x0 - x1/2 - x3/2,
R/2 + y0 - y3/2 - 400,
]
end;
iniv = [-1500.0, 900, 2000, -1000, 800, 2200, 900, -100, -400, 300]
res = nls(H, ini=iniv)
println(res)
([-1700.0, 873.9999999999999, 2025.75, -1089.0000042590011, 816.7499943213318, 2206.25, 900.0, -106.24999999999997, -413.00000000000006, 309.74999999999994], true)
外円の半径は 2206.25 である。元の単位での径(直径)は 22 寸 0 分 6 厘 2 毛 5 糸である。
using Plots
using Printf
function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
θ = beginangle:0.1:endangle
x = r.*cosd.(θ)
y = r.*sind.(θ)
if fill
plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
else
plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end
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=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(x1, x3, y3, x, y, R, x4, y4, x0, y0) = res[1]
println("R = $R")
plot()
circle(0, 0, R)
circle(x, y, 845, :blue)
circle(x4, y4, 1300, :magenta)
circle(0, 400-R, 400, :green)
segment(x1, 800 - R, -x1, 800 - R)
segment(x1, 800 - R, x3, y3)
if more
point(x4, y4, "(x4,y4)")
point(0, 400-R, "400-R")
point(x, y, "(x,y)")
point(x0, y0, "(x0,y0)")
point(x1, 800 - R, "(x1,800-R)")
point(x3, y3, "(x3,y3)")
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
end
end;
※コメント投稿者のブログIDはブログ作成者のみに通知されます