算額(その2009)
(22) 兵庫県養父町左近山 地蔵堂 明治21年(1888)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.
キーワード:円4個,長方形,斜線2本
長方形の中に 1 本の対角線ともう 1 本の斜線を引き,等円 4 個を容れる。長方形の長辺と短辺が与えられたとき,等円の直径を求めよ。
長方形の長辺と短辺を a, b
斜線と長方形の長辺の交点座標を (c, 0), (a - c, b)
等円の半径と中心座標を r, (x1, r), (x2, y2), (r, y3), (x3, b - r)
とおき,以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms a::positive, b::positive, c::positive,
r::positive, x1::positive, x2::positive,
y2::positive, y3::positive, x3::positive
eq1 = dist2(0, b, a, 0, x1, r, r)
eq2 = dist2(0, b, a, 0, x2, y2, r)
eq3 = dist2(0, b, a, 0, r, y3, r)
eq4 = dist2(0, b, a, 0, x3, b - r, r)
eq5 = dist2(c, 0, a - c, b, x2, y2, r)
eq6 = dist2(c, 0, a - c, b, x1, r, r)
eq7 = (x2 - r)^2 + (y2 - y3)^2 - 4r^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)
(c, r, x1, x2, y2, y3, x3) = u
return [
b*(a^2*b - 2*a^2*r - 2*a*b*x1 + 2*a*r*x1 - b*r^2 + b*x1^2), # eq1
a^2*b^2 - 2*a^2*b*y2 - a^2*r^2 + a^2*y2^2 - 2*a*b^2*x2 + 2*a*b*x2*y2 - b^2*r^2 + b^2*x2^2, # eq2
a*(a*b^2 - 2*a*b*y3 - a*r^2 + a*y3^2 - 2*b^2*r + 2*b*r*y3), # eq3
b*(-2*a*r*x3 - b*r^2 + b*x3^2), # eq4
-a^2*r^2 + a^2*y2^2 + 2*a*b*c*y2 - 2*a*b*x2*y2 + 4*a*c*r^2 - 4*a*c*y2^2 + b^2*c^2 - 2*b^2*c*x2 - b^2*r^2 + b^2*x2^2 - 4*b*c^2*y2 + 4*b*c*x2*y2 - 4*c^2*r^2 + 4*c^2*y2^2, # eq5
b*(2*a*c*r - 2*a*r*x1 + b*c^2 - 2*b*c*x1 - b*r^2 + b*x1^2 - 4*c^2*r + 4*c*r*x1), # eq6
-4*r^2 + (-r + x2)^2 + (y2 - y3)^2, # eq7
]
end;
(a, b) = (10, 5)
iniv = BigFloat[2.6, 1.1, 5.2, 3.1, 2.2, 3.2, 4.8]
res = nls(H, ini=iniv)
([2.60785001480277, 1.1239440727077343, 5.238896505102071, 3.1345163522937027, 2.1761341491119186, 3.181420288904903, 4.761103494897929], true)
例えば,長方形の長辺と短辺が 10 寸,5 寸のとき,等円の直径は 2*1.1239440727077343 = 2.2478881454154687 である。
function draw(a, b, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(c, r, x1, x2, y2, y3, x3) = res[1]
plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:green, lw=0.5)
circle(x1, r, r)
circle(x2, y2, r)
circle(r, y3, r)
circle(x3, b - r, r)
segment(0, b, a, 0, :blue)
segment(c, 0, a - c, b, :blue)
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(x1, r, "(x1,r)", :red, :center, delta=-delta/2)
point(x2, y2, "(x2,y2)", :red, :center, delta=-delta/2)
point(r, y3, "(r,y3)", :red, :center, delta=-delta/2)
point(x3, b - r, "(x3,b-r)", :red, :center, delta=-delta/2)
point(c, 0, "c", :blue, :center, delta=-delta/2)
point(a, 0, "a", :blue, :center, delta=-delta/2)
point(0, b, "b ", :blue, :right, :vcenter)
point(a - c, b, "(a-c,b)", :blue, :center, :bottom, delta=delta/2)
point(a, b, "(a,b)", :blue, :right, :bottom, delta=delta/2)
plot!(xlims=(-7delta, a + 5delta), ylims=(-7delta, b + 5delta))
end
end;
draw(a, b, true)