算額(その247)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3(229)
長野県諏訪市下諏訪 諏訪大社下社 慶応4年(1868)
菱形の中に日円,月円,星円がそれぞれ 1 個,3 個,2 個入っている。
日円,月円の径を知って,星円の径を求めよ。
![](https://blogimg.goo.ne.jp/user_image/20/1d/1baa45ce7253f268803a1b73ad57da33.png)
菱形の中心を原点とする。x 軸上の頂点の x 座標を x,y 軸上の頂点の y 座標を y とする。
日円の半径,中心座標を r1, (0, y1) とする。
真ん中の月円の半径,中心座標を r2, (0, y2) とする。
右の月円の半径,中心座標を r3, (x3, y3) とする。r3 = r2 である。
星円の半径,中心座標を r4, (x4, r4) とする。
以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms r1::positive, y1::positive, r2::positive, y2::negative,
r3::positive, x3::positive, y3::positive,
r4::positive, x4::positive, y4::positive,
x::positive, y::positive;
r3 = r2
eq1 = x4^2 + (y1 - y4)^2 - (r1 + r4)^2
eq2 = (x3 - x4)^2 + (y3 - y4)^2 - (r3 + r4)^2
eq3 = x3^2 + (y1 - y3)^2 - (r1 + r3)^2
eq4 = x3^2 + (y2 - y3)^2 - (r2 + r3)^2
eq5 = distance(0, y, x, 0, 0, y1) - r1^2
eq6 = distance(0, y, x, 0, x4, y4) - r4^2
eq7 = distance(0, -y, x, 0, 0, y2) - r2^2
eq8 = distance(0, -y, x, 0, x3, y3) - r3^2
eq9 = (sqrt(x^2 + y^2)/x + 1) * (r1 + r2) - 2y;
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9])
println(eq1, ", # eq1")
println(eq2, ", # eq2")
println(eq3, ", # eq3")
println(eq4, ", # eq4")
println(eq5, ", # eq5")
println(eq6, ", # eq6")
println(eq7, ", # eq7")
println(eq8, ", # eq8")
println(eq9, ", # eq9")
x4^2 - (r1 + r4)^2 + (y1 - y4)^2, # eq1
-(r2 + r4)^2 + (x3 - x4)^2 + (y3 - y4)^2, # eq2
x3^2 - (r1 + r2)^2 + (y1 - y3)^2, # eq3
-4*r2^2 + x3^2 + (y2 - y3)^2, # eq4
-r1^2 + x^2*y^2*(y - y1)^2/(x^2 + y^2)^2 + (-y*(x^2 + y*y1)/(x^2 + y^2) + y1)^2, # eq5
-r4^2 + (-x*(x*x4 + y^2 - y*y4)/(x^2 + y^2) + x4)^2 + (-y*(x^2 - x*x4 + y*y4)/(x^2 + y^2) + y4)^2, # eq6
-r2^2 + x^2*y^2*(y + y2)^2/(x^2 + y^2)^2 + (-y*(-x^2 + y*y2)/(x^2 + y^2) + y2)^2, # eq7
-r2^2 + (-x*(x*x3 + y^2 + y*y3)/(x^2 + y^2) + x3)^2 + (-y*(-x^2 + x*x3 + y*y3)/(x^2 + y^2) + y3)^2, # eq8
-2*y + (1 + sqrt(x^2 + y^2)/x)*(r1 + r2), # eq9
日円,月円の半径を 5, 4 とする。
using Printf
using Plots
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)
(r1, r2) = (5, 4)
(r4, y1, y2, x3, y3, x4, y4, y, x) = u
return [
x4^2 - (r1 + r4)^2 + (y1 - y4)^2, # eq1
-(r2 + r4)^2 + (x3 - x4)^2 + (y3 - y4)^2, # eq2
x3^2 - (r1 + r2)^2 + (y1 - y3)^2, # eq3
-4*r2^2 + x3^2 + (y2 - y3)^2, # eq4
-r1^2 + x^2*y^2*(y - y1)^2/(x^2 + y^2)^2 + (-y*(x^2 + y*y1)/(x^2 + y^2) + y1)^2, # eq5
-r4^2 + (-x*(x*x4 + y^2 - y*y4)/(x^2 + y^2) + x4)^2 + (-y*(x^2 - x*x4 + y*y4)/(x^2 + y^2) + y4)^2, # eq6
-r2^2 + x^2*y^2*(y + y2)^2/(x^2 + y^2)^2 + (-y*(-x^2 + y*y2)/(x^2 + y^2) + y2)^2, # eq7
-r2^2 + (-x*(x*x3 + y^2 + y*y3)/(x^2 + y^2) + x3)^2 + (-y*(-x^2 + x*x3 + y*y3)/(x^2 + y^2) + y3)^2, # eq8
-2*y + (1 + sqrt(x^2 + y^2)/x)*(r1 + r2), # eq9
]
end;
iniv = [big"1.0", 1.9, -3.1, 3.1, -1.3, 4, 1.9, 5.4, 9.3]
res = nls(H, ini=iniv);
println(res[2] ? "収束した" : "収束していない")
収束した
names = ["r4", "y1", "y2", "x3", "y3", "x4", "y4", "y ", "x "]
for i in 1:length(names)
println("$(names[i]) = $(Float64(res[1][i]))")
end
r4 = 1.7470608602667996
y1 = 3.9418436943485617
y2 = -5.058156305651438
x3 = 7.166451331820933
y3 = -1.5026007500958825
x4 = 6.740960675000237
y4 = 4.228687607033678
y = 9.523406750862943
x = 19.195039966835868
星円の半径は 1.7470608602667996。
![](https://blogimg.goo.ne.jp/user_image/05/0f/d345497d10f553d41f208c82478d90a3.png)
using Plots
using Printf
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r4, y1, y2, x3, y3, x4, y4, y, x) = res[1]
(r1, r2) = (5, 4)
r3 = r2
# y = 5/2 + 5*sqrt(3)/3
# x = sqrt(3) * y
@printf("r1 = %.5f; r2 = %.5f; \nr4 = %.5f; y1 = %.5f; y2 = %.5f; \nx3 = %.5f; y3 = %.5f; x4 = %.5f; \ny4 = %.5f; y = %.5f; x = %.5f\n", r1, r2, r4, y1, y2, x3, y3, x4, y4, y, x)
plot([x, 0, -x, 0, x], [0, y, 0, -y, 0], color=:black, lw=0.5)
circle(0, y1, r1)
circle(0, y2, r2, :blue)
circle(x3, y3, r3, :blue)
circle(-x3, y3, r3, :blue)
circle(x4, y4, r4, :green)
circle(-x4, y4, r4, :green)
if more == true
point(0, y1, " y1\n 日円", :red, :left, :bottom)
point(0, y2, " y2\n 月円", :blue, :left, :bottom)
point(x3, y3, "(x3,y3)\n月円", :blue, :left, :center)
point(x4, y4, "(x4,y4)\n星円", :green, :left, :center)
point(0, y, " y", :black, :left, :bottom)
point(x, 0, "x", :black)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;