算額（その697）
神壁算法 關流清水與市道香門人 上総國長柄郡上之郷 本石與八利重 寛政七年
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
不等辺三角形の辺の長さが長い順に，大斜，中斜，小斜がそれぞれ 345 寸，322 寸，299 寸である。三角形内に斜線を 2 本引き，区画された領域に甲円と乙円をいれる。甲円の直径が 115 寸のとき，乙円の直径はいかほどか。
三角形の辺の長さを「大斜」，「中斜」，「小斜」，中斜と小斜の交点座標を (x0, y0)，斜線と小斜，中斜との交点座標を (xa, ya), (xb, yb) とする。
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, y2)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms 大斜::positive, 中斜::positive, 小斜::positive,
x0::positive, y0::positive,
xa::positive, xb::positive,
r1::positive, x1::positive,
r2::positive, x2::positive, y2::positive
ya = y0/x0 * xa
yb = y0/(大斜 - x0) * (大斜 - xb)
eq1 = x0^2 + y0^2 - 小斜^2
eq2 = (大斜 - x0)^2 + y0^2 - 中斜^2
eq3 = distance(x0, y0, 大斜, 0, x2, y2) - r2^2
eq4 = distance(0, 0, x0, y0, x2, y2) - r2^2
eq5 = distance(0, 0, xb, yb, x1, r1) - r1^2
eq6 = distance(0, 0, xb, yb, x2, y2) - r2^2
eq7 = distance(xa, ya, 大斜, 0, x1, r1) - r1^2
eq8 = distance(xa, ya, 大斜, 0, x2, y2) - r2^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 v, r.f_converged
end;
function H(u)
(x0, y0, xa, xb, x1, r2, x2, y2) = u
return [
x0^2 + y0^2 - 小斜^2, # eq1
y0^2 - 中斜^2 + (-x0 + 大斜)^2, # eq2
-r2^2 + (x2 - (x0^2*x2 - 2*x0*x2*大斜 + x0*y0*y2 + x2*大斜^2 + y0^2*大斜 - y0*y2*大斜)/(x0^2 - 2*x0*大斜 + y0^2 + 大斜^2))^2 + (-y0*(x0*x2 - x0*大斜 - x2*大斜 + y0*y2 + 大斜^2)/(x0^2 - 2*x0*大斜 + y0^2 + 大斜^2) + y2)^2, # eq3
-r2^2 + (-x0*(x0*x2 + y0*y2)/(x0^2 + y0^2) + x2)^2 + (-y0*(x0*x2 + y0*y2)/(x0^2 + y0^2) + y2)^2, # eq4
-r1^2 + (r1 - y0*(r1*xb^2*y0 - 2*r1*xb*y0*大斜 + r1*y0*大斜^2 + x0*x1*xb^2 - x0*x1*xb*大斜 - x1*xb^2*大斜 + x1*xb*大斜^2)/(x0^2*xb^2 - 2*x0*xb^2*大斜 + xb^2*y0^2 + xb^2*大斜^2 - 2*xb*y0^2*大斜 + y0^2*大斜^2))^2 + (x1 - xb*(r1*x0*xb*y0 - r1*x0*y0*大斜 - r1*xb*y0*大斜 + r1*y0*大斜^2 + x0^2*x1*xb - 2*x0*x1*xb*大斜 + x1*xb*大斜^2)/(x0^2*xb^2 - 2*x0*xb^2*大斜 + xb^2*y0^2 + xb^2*大斜^2 - 2*xb*y0^2*大斜 + y0^2*大斜^2))^2, # eq5
-r2^2 + (x2 - xb*(x0^2*x2*xb - 2*x0*x2*xb*大斜 + x0*xb*y0*y2 - x0*y0*y2*大斜 + x2*xb*大斜^2 - xb*y0*y2*大斜 + y0*y2*大斜^2)/(x0^2*xb^2 - 2*x0*xb^2*大斜 + xb^2*y0^2 + xb^2*大斜^2 - 2*xb*y0^2*大斜 + y0^2*大斜^2))^2 + (-y0*(x0*x2*xb^2 - x0*x2*xb*大斜 - x2*xb^2*大斜 + x2*xb*大斜^2 + xb^2*y0*y2 - 2*xb*y0*y2*大斜 + y0*y2*大斜^2)/(x0^2*xb^2 - 2*x0*xb^2*大斜 + xb^2*y0^2 + xb^2*大斜^2 - 2*xb*y0^2*大斜 + y0^2*大斜^2) + y2)^2, # eq6
-r1^2 + (r1 - xa*y0*(r1*xa*y0 + x0*x1*xa - x0*x1*大斜 - x0*xa*大斜 + x0*大斜^2)/(x0^2*xa^2 - 2*x0^2*xa*大斜 + x0^2*大斜^2 + xa^2*y0^2))^2 + (x1 - (r1*x0*xa^2*y0 - r1*x0*xa*y0*大斜 + x0^2*x1*xa^2 - 2*x0^2*x1*xa*大斜 + x0^2*x1*大斜^2 + xa^2*y0^2*大斜)/(x0^2*xa^2 - 2*x0^2*xa*大斜 + x0^2*大斜^2 + xa^2*y0^2))^2, # eq7
-r2^2 + (x2 - (x0^2*x2*xa^2 - 2*x0^2*x2*xa*大斜 + x0^2*x2*大斜^2 + x0*xa^2*y0*y2 - x0*xa*y0*y2*大斜 + xa^2*y0^2*大斜)/(x0^2*xa^2 - 2*x0^2*xa*大斜 + x0^2*大斜^2 + xa^2*y0^2))^2 + (-xa*y0*(x0*x2*xa - x0*x2*大斜 - x0*xa*大斜 + x0*大斜^2 + xa*y0*y2)/(x0^2*xa^2 - 2*x0^2*xa*大斜 + x0^2*大斜^2 + xa^2*y0^2) + y2)^2, # eq8
]
end;
(大斜, 中斜, 小斜, r1) = (345, 322, 299, 115//2)
iniv = BigFloat[152, 258, 105, 219, 165, 42, 162, 182]
res = nls(H, ini=iniv)
(BigFloat[151.8000000000000000000000000000000000000000000000000000000000000000000000000004, 257.6000000000000000000000000000000000000000000000000000000000000000000000000009, 100.0, 213.75, 160.9999999999999999999999999999999999999999999999999999999999999999999999999956, 42.00000000000000000000000000000000000000000000000000000000000000000000000000111, 156.0, 182.0], true)
乙円の直径は 84 寸である。
その他のパラメータは以下のとおり。
x0 = 151.8; y0 = 257.6; xa = 100; xb = 213.75; x1 = 161; r2 = 42; x2 = 156; y2 = 182
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(大斜, 中斜, 小斜, r1) = (345, 322, 299, 115//2)
(x0, y0, xa, xb, x1, r2, x2, y2) = res[1]
println("乙円の直径 = $(Float64(2r2))")
@printf("x0 = %g; y0 = %g; xa = %g; xb = %g; x1 = %g; r2 = %g; x2 = %g; y2 = %g\n",
x0, y0, xa, xb, x1, r2, x2, y2)
ya = y0/x0 * xa
yb = y0/(大斜 - x0) * (大斜 - xb)
plot([0, 大斜, x0, 0], [0, 0, y0, 0], color=:black, lw=0.5)
circle(x1, r1, r1)
circle(x2, y2, r2, :blue)
segment(0, 0, xb, yb, :green)
segment(xa, ya, 大斜, 0, :green)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
point(大斜, 0, "大斜", :black, :left, :bottom, delta=delta/2)
point(x0, y0, "(x0,y0)", :black, :left, :bottom, delta=delta/2)
point(xa, ya, "(xa,ya)", :black, :right, :bottom, delta=delta/2)
point(xb, yb, "(xb,yb)", :black, :left, :bottom, delta=delta/2)
point(x1, r1, "甲円:r1,(x1,r1)", :red, :center, :top, delta=-delta/2)
point(x2, y2, "乙円:r2,(x2,y2)", :blue, :center, :top, delta=-delta/2)
end
end;
