算額(その99)
東京都渋谷区渋谷 金王八幡神社 安政6年(1859)4月
http://www.wasan.jp/tokyo/konnou2.html
参照先の画像では問題文が全く読めない
下図のようjに記号を定め,方程式を解く。
しかし,解は複素数になる。もっとも,虚部はほとんど 0 なのだが,厳密に言えば「@syms R::positive」などとした場合には,「解なし」ということになる。
解の式は恐ろしく複雑になるが,数値にすると以下のようになる。最初の組が適切な解である。
using SymPy
@syms R::positive, r::positive, x::positive, y::positive;
@syms R, r, x, y;
eq1 = r^2 + y^2 - (r + 593)^2;
eq2 = r^2 + (y - R)^2 - (R - r)^2;
eq3 = (593 - r)^2 + R^2 - (R + r)^2;
res = solve([eq1, eq2, eq3], (R, r, y))
3-element Vector{Tuple{Sym, Sym, Sym}}:
((-1186*sqrt(3) - 2965*2^(2/3)*I*(7 + 21*sqrt(3)*I)^(1/3)/6 - 593*2^(1/3)*I*(7 + 21*sqrt(3)*I)^(2/3)/6 + 1186*I/3 + 593*2^(2/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(1/3)/6 + 593*2^(1/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(2/3)/6)/(3*sqrt(3) - I), (593*sqrt(3) - 593*2^(2/3)*I*(7 + 21*sqrt(3)*I)^(1/3)/3 - 593*2^(2/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(1/3)/6 - 593*I/3 + 593*2^(1/3)*I*(7 + 21*sqrt(3)*I)^(2/3)/12 + 593*2^(1/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(2/3)/12)/(3*sqrt(3) - I), -593/3 + 4151/(9*(7/54 + 7*sqrt(3)*I/18)^(1/3)) + 593*(7/54 + 7*sqrt(3)*I/18)^(1/3))
((-1186*sqrt(3) - 4151*sqrt(3)*(28 + 84*sqrt(3)*I)^(1/3)/18 - 593*2^(1/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(2/3)/18 - 1186*2^(2/3)*I*(7 + 21*sqrt(3)*I)^(1/3)/9 - 593*2^(2/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(1/3)/9 + 593*98^(1/3)*I*(1 + 3*sqrt(3)*I)^(2/3)/18 + 1186*I/3 + 593*sqrt(3)*98^(1/3)*(1 + 3*sqrt(3)*I)^(2/3)/18 + 4151*I*(28 + 84*sqrt(3)*I)^(1/3)/18 + 2965*2^(1/3)*I*(7 + 21*sqrt(3)*I)^(2/3)/18)/(3*sqrt(3) - I), (593*sqrt(3) - 593*sqrt(3)*98^(1/3)*(1 + 3*sqrt(3)*I)^(2/3)/18 - 4151*sqrt(3)*(28 + 84*sqrt(3)*I)^(1/3)/36 - 593*2^(1/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(2/3)/36 - 593*98^(1/3)*I*(1 + 3*sqrt(3)*I)^(2/3)/18 - 593*I/3 + 593*2^(2/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(1/3)/9 + 4151*I*(28 + 84*sqrt(3)*I)^(1/3)/36 + 1186*2^(2/3)*I*(7 + 21*sqrt(3)*I)^(1/3)/9 + 2965*2^(1/3)*I*(7 + 21*sqrt(3)*I)^(2/3)/36)/(3*sqrt(3) - I), -593/3 + 593*(-1/2 - sqrt(3)*I/2)*(7/54 + 7*sqrt(3)*I/18)^(1/3) + 4151/(9*(-1/2 - sqrt(3)*I/2)*(7/54 + 7*sqrt(3)*I/18)^(1/3)))
((-1186*sqrt(3) - 593*2^(1/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(2/3)/9 - 1186*2^(1/3)*I*(7 + 21*sqrt(3)*I)^(2/3)/9 - 593*sqrt(3)*98^(1/3)*(1 + 3*sqrt(3)*I)^(2/3)/18 - 593*2^(2/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(1/3)/18 + 593*98^(1/3)*I*(1 + 3*sqrt(3)*I)^(2/3)/18 + 1186*I/3 + 2965*2^(2/3)*I*(7 + 21*sqrt(3)*I)^(1/3)/18 + 4151*sqrt(3)*(28 + 84*sqrt(3)*I)^(1/3)/18 + 4151*I*(28 + 84*sqrt(3)*I)^(1/3)/18)/(3*sqrt(3) - I), (593*sqrt(3) - 2965*2^(2/3)*I*(7 + 21*sqrt(3)*I)^(1/3)/18 - 593*2^(1/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(2/3)/18 - 593*2^(1/3)*I*(7 + 21*sqrt(3)*I)^(2/3)/9 - 593*98^(1/3)*I*(1 + 3*sqrt(3)*I)^(2/3)/18 - 593*I/3 + 593*2^(2/3)*sqrt(3)*(7 + 21*sqrt(3)*I)^(1/3)/18 + 4151*sqrt(3)*(28 + 84*sqrt(3)*I)^(1/3)/36 + 4151*I*(28 + 84*sqrt(3)*I)^(1/3)/36 + 593*sqrt(3)*98^(1/3)*(1 + 3*sqrt(3)*I)^(2/3)/18)/(3*sqrt(3) - I), -593/3 + 4151/(9*(-1/2 + sqrt(3)*I/2)*(7/54 + 7*sqrt(3)*I/18)^(1/3)) + 593*(-1/2 + sqrt(3)*I/2)*(7/54 + 7*sqrt(3)*I/18)^(1/3))
for j = 1:length(res)
println("解 $j")
println("R = $(res[j][1].evalf())")
println("r = $(res[j][2].evalf())")
println("y = $(res[j][3].evalf())")
println()
end
解 1
R = 475.549077332269 - 6.93889390390723e-18*I
r = 164.545086163906
y = 739.458905004458 + 0.e-20*I
解 2
R = -1332.45890500446 + 1.38777878078145e-17*I
r = -237.774538666135 + 3.46944695195361e-18*I
y = -263.909827672189 + 0.e-17*I
解 3
R = -329.090172327811
r = 666.229452502229 - 6.93889390390723e-18*I
y = -1068.54907733227 + 0.e-18*I
解説したページ https://sites.google.com/site/seijiimura/home/09-san-gaku-no-mondai-ni-chousen/7-jin-wang-ba-fan-shen-sheno-suan-e は
凡ミスがあるが,解としては正しい。ただし,右上の小円の中心座標はさらに求める必要がある。
ここで示した SymPy による方法では,必要な数値は(一度で)全て得られる。
算額では 463 寸ばかりとなっているが,中円の径は 475.549077332269 寸に間違いない。
using Plots
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=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="")#, fontfamily="IPAMincho")
R = 475.549077332269
r = 164.545086163906
y = 739.458905004458
plot()
circle(0, 0, 593, :green)
circle(0, R, R, :blue)
circle(0, -R, R, :blue)
circle(593 - r, 0, r)
circle(-593 + r, 0, r)
circle(r, y, r)
circle(-r, y, r)
circle(r, -y, r)
circle(-r, -y, r)
if more
point(0, R, "R ", :blue, :right)
point(0, 593, "593 ", :green, :right)
point(593-r, 0, "593-r ", :red, :top)
point(r, y, "(r,y)", :red, :top)
point(0, 0, "0 ", :brown, :top, :right)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
end
end;