算額(その188)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1(86)
長野県中野市田上 田上観音堂 文化6年(1809)
第4問 鉤股(直角三角形)内に大中小の3円を入れる。大円と中円の径を掛けると 72,大円と小円の径を掛けると 48 であるとき,それぞれの径を求めよ。
注:算額の図は「算額(その187)https://blog.goo.ne.jp/r-de-r/e/2053d2ea5ef9100ee13ecf80457132b9」と違い垂直線,水平線に接するという条件がない場合の解である。
大円,中円,小円の半径を r1,r2,r3 とする。中円,小円の中心座標を (x2, r2), (r3, y3) とする。
中円,小円の中心から弦への距離,鉤股弦の面積の関係式を解く。
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 X::positive, Y::positive, x2::pisitive, y3::positive,
r1::positive, r2::positive, r3::positive, l::positive;
eq1 = 2r1 * 2r2 - 72
eq2 = 2r1 * 2r3 - 48
eq3 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq4 = (r1 - r3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq5 = distance(0, Y, X, 0, r3, y3) - r3^2
eq6 = distance(0, Y, X, 0, x2, r2) - r2^2
eq7 = X^2 + Y^2 - l^2
eq8 = (X + Y + l)*r1 - X * Y;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8])
solve() では解けないので,nlsolve() を使用する。
println(eq1, ",")
println(eq2, ",")
println(eq3, ",")
println(eq4, ",")
println(eq5, ",")
println(eq6, ",")
println(eq7, ",")
println(eq8, ",")
4*r1*r2 - 72,
4*r1*r3 - 48,
(-r1 + x2)^2 + (r1 - r2)^2 - (r1 + r2)^2,
(-r1 + y3)^2 + (r1 - r3)^2 - (r1 + r3)^2,
-r3^2 + (-X*(X*r3 + Y^2 - Y*y3)/(X^2 + Y^2) + r3)^2 + (-Y*(X^2 - X*r3 + Y*y3)/(X^2 + Y^2) + y3)^2,
-r2^2 + (-X*(X*x2 + Y^2 - Y*r2)/(X^2 + Y^2) + x2)^2 + (-Y*(X^2 - X*x2 + Y*r2)/(X^2 + Y^2) + r2)^2,
X^2 + Y^2 - l^2,
-X*Y + r1*(X + Y + l),
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, r3, X, Y, l, x2, y3) = u
return [
4*r1*r2 - 72,
4*r1*r3 - 48,
(-r1 + x2)^2 + (r1 - r2)^2 - (r1 + r2)^2,
(-r1 + y3)^2 + (r1 - r3)^2 - (r1 + r3)^2,
-r3^2 + (-X*(X*r3 + Y^2 - Y*y3)/(X^2 + Y^2) + r3)^2 + (-Y*(X^2 - X*r3 + Y*y3)/(X^2 + Y^2) + y3)^2,
-r2^2 + (-X*(X*x2 + Y^2 - Y*r2)/(X^2 + Y^2) + x2)^2 + (-Y*(X^2 - X*x2 + Y*r2)/(X^2 + Y^2) + r2)^2,
X^2 + Y^2 - l^2,
-X*Y + r1*(X + Y + l),
]
end;
iniv = [big"6.1", 4.2, 3.3, 24, 17, 31, 15, 13]
res = nls(H, ini=iniv)
println(res)
(BigFloat[5.748772323345360879553386071530534002107771135435234241146154053030068382725282, 3.131103301291523318846022596310963240657188641774061152337086444494206214742888, 2.08740220086101554589734764786891092285620654321669230693296468944296425038439, 24.38365332206977707890963182337988451889403111451717974855438930320079233053448, 16.62684846651738071918895370581186617641229049594184966902009096464686858816747, 29.51295714189643603899180810391979664284976102613029457759886594916383321385045, 14.23405369758393117236351841678872247352580238769692473134187249862639602572867, 12.67697555362087005366317143755402346987899215067675675336938195307101120791379], true)
using Printf
@printf("大円径 = %6.3f; 中円径 = %6.3f; 小円径 = %6.3f\n", 2res[1][1], 2res[1][2], 2res[1][3])
大円径 = 11.498; 中円径 = 6.262; 小円径 = 4.175
大円,中円,小円の径は元の単位で,11.498 寸,6.262 寸,4.175 寸である。
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, linecolor=:black, linewidth=0.5, seriestype=:shape, fillcolor=color)
else
plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
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)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, r3, X, Y, l, x2, y3) = res[1]
@printf("X = %6.3f; Y = %6.3f; l = %6.3f; x2 = %6.3f; y3 = %6.3f\n",
X, Y, l, x2, y3)
@printf("r1 = %6.3f; r2 = %6.3f; r3 = %6.3f\n", r1, r2, r3)
plot([0, X, 0, 0], [0, 0, Y, 0], color=:black, lw=0.5)
circle(r1, r1, r1, :red)
circle(x2, r2, r2, :blue)
circle(r3, y3, r3, :green)
if more == true
point(r1, r1, "大円(r1,r1,r1)")
point(x2, r2, "中円(x2,r2,r2)")
point(r3, y3, "小円(r3,y3,r3)")
point(0, 0, " O", :green, :left, :bottom)
point(X, 0, " X", :green, :left, :bottom)
point(0, Y, " Y", :green, :left, :bottom)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
X = 24.384; Y = 16.627; l = 29.513; x2 = 14.234; y3 = 12.677
r1 = 5.749; r2 = 3.131; r3 = 2.087