算額（その767）
宮城県大和町落合舞野 舞野正観音堂 慶応4年(1868)
http://www.wasan.jp/miyagi/mainokannondo.html
正方形の中に交差する楕円が入っており，区画された領域に4個の等円が入っている。正方形の一辺の長さが与えられたとき，等円の直径を求めよ。
正方形の一辺の長さを 2a とする。
楕円の長半径と短半径，中心座標を a, b, (0, 0)
等円の半径と中心座標を r, (a - r, a - r), (0, 0)
とする。楕円の長半径は正方形の一辺の長さの半分，等円の半径と楕円の短半径は等しい。
以下の連立方程式を解く。
しかし，たかが3元連立方程式なのに SymPy では答えを得られない。
include("julia-source.txt");
using SymPy
@syms a::positive, b::positive, r::positive,
x0::positive, y0::positive
@syms a, b, r, x0, y0
a = 1//2
r = b
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*x0/(a^2*y0) + (a - r - x0)/(a - r - y0)
eq3 = (a - r - x0)^2 + (a - r - y0)^2 - r^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 Float64.(v), r.f_converged
end;
function H(u)
(b, x0, y0) = u
return [
4*x0^2 - 1 + y0^2/b^2, # eq1
-4*b^2*x0/y0 + (-b - x0 + 1/2)/(-b - y0 + 1/2), # eq2
-b^2 + (-b - x0 + 1/2)^2 + (-b - y0 + 1/2)^2, # eq3
]
end;
a = 1/2
iniv = BigFloat[0.356, 0.57, 0.295]./2
res = nls(H, ini=iniv)
([0.1785016026524881, 0.27969014973090045, 0.14796196722191884], true)
等円の直径は 正方形の一辺の長さの 0.3570032053049762 倍である。
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
a = 1/2
(b, x0, y0) = [0.1785016026524881, 0.27969014973090045, 0.14796196722191884]
r = b
println("等円の直径は 正方形の一辺の長さの $(2r) 倍である。")
plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:black, lw=0.5)
ellipse(0, 0, a, b, color=:blue)
ellipse(0, 0, b, a, color=:blue)
circle(0, 0, r)
circle4(a - r, a - r, r)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
point(0, b, "r = b", :blue, :center, :bottom, delta=delta)
point(0, 0, "等円:r,(0,0)", :red, :center, delta=-delta)
point(a - r, a - r, "等円:r,(a-r,a-r)", :red, :center, delta=-delta)
point(a, 0, "a ", :blue, :right, :bottom, delta=delta)
end
end;
