算額（その440）
東都下谷摩利支天堂 天明8年戊申3月(1788)
東都浅草新堀端住 三浦伴次郎成喜
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
求めるものが小頭か大頭かの違いだけで，算額（その185）と同じ問題である。プログラミングスタイルを最近のものに変えた。
等脚台形の隔斜（対角線）で区切られた領域に甲，乙，丙の円を入れる。甲円の直径が 100 寸，乙円の直径が 28 寸，丙円の直径が 45 寸のとき，大頭（台形の下底）はいかほどか。
以下のように記号を定め方程式を解く。
台形の右下隅 A，右上隅 b， の座標をそれぞれ (a, 0)，(b, y) と置く。
甲円，乙円，丙円の半径をそれぞれ，r1，r2，r3 とする。
甲円，乙円，丙円の中心から隔斜（対角線）までの距離が円の半径に等しいという連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms a::positive, b::positive, x1::positive, y1::positive, y::positive,
r1::positive, r2::positive, r3::positive;
eq1 = distance(b, y, a, 0, x1, y1) - r3^2
eq2 = distance(-b, y, a, 0, x1, y1) - r3^2
eq3 = distance(-a, 0, b, y, x1, y1) - r3^2
eq4 = distance(-a, 0, b, y, 0, y - r2) - r2^2
eq5 = distance(-b, y, a, 0, 0, r1) - r1^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)
(a, b, x1, y1, y) = u
return [
-r3^2 + (x1 - (a^2*x1 - 2*a*b*x1 + a*y^2 - a*y*y1 + b^2*x1 + b*y*y1)/(a^2 - 2*a*b + b^2 + y^2))^2 + (-y*(a^2 - a*b - a*x1 + b*x1 + y*y1)/(a^2 - 2*a*b + b^2 + y^2) + y1)^2, # eq1
-r3^2 + (x1 - (a^2*x1 + 2*a*b*x1 + a*y^2 - a*y*y1 + b^2*x1 - b*y*y1)/(a^2 + 2*a*b + b^2 + y^2))^2 + (-y*(a^2 + a*b - a*x1 - b*x1 + y*y1)/(a^2 + 2*a*b + b^2 + y^2) + y1)^2, # eq2
-r3^2 + (x1 - (x1*(a^2 + 2*a*b + b^2 + y^2) - y*(a*y - a*y1 - b*y1 + x1*y))/(a^2 + 2*a*b + b^2 + y^2))^2 + (-y*(a^2 + a*b + a*x1 + b*x1 + y*y1)/(a^2 + 2*a*b + b^2 + y^2) + y1)^2, # eq3
-r2^2 + y^2*(-a*r2 - b*r2 + b*y)^2/(a^2 + 2*a*b + b^2 + y^2)^2 + (-r2 - y*(a^2 + a*b - r2*y + y^2)/(a^2 + 2*a*b + b^2 + y^2) + y)^2, # eq4
-r1^2 + y^2*(-a*r1 + a*y - b*r1)^2/(a^2 + 2*a*b + b^2 + y^2)^2 + (r1 - y*(a^2 + a*b + r1*y)/(a^2 + 2*a*b + b^2 + y^2))^2, # eq5
]
end;
(r1, r2, r3) =(100, 28, 45) ./ 2
iniv = [big"151.0", 44, 39, 112, 144]
res = nls(H, ini=iniv)
names = ("a", "b", "x1", "y1", "y")
if res[2]
for (name, x) in zip(names, res[1])
@printf("%s = %g; ", name, round(Float64(x), digits=6))
end
println()
else
println("収束していない")
end
a = 150; b = 42; x1 = 37.5; y1 = 112.5; y = 144;
大頭の長さは 2a = 300（寸）である。
using Plots
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, r3) =(100, 28, 45) ./ 2
(a, b, x1, y1, y) = res[1]
@printf("r1 = %g; r2 = %g; r3 = %g\n", r1, r2, r3)
@printf("a = %g; b = %g; x1 = %g; y1 = %g; y = %g\n", a, b, x1, y1, y)
plot([a, b, -b, -a, a], [0, y, y, 0, 0], color=:orange, lw=0.5)
circle(0, r1, r1, :green)
circle(x1, y1, r3, :blue)
circle(-x1, y1, r3, :blue)
circle(0, y - r2, r2, :red)
segment(-b, y, a, 0, :orange)
segment(-a, 0, b, y, :orange)
if more == true
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(0, r1, " 甲円:r1,(0,r1)", :black, :left, :vcenter)
point(0, y - r2, " 乙円:r2,(0,y-r2)", :black, :left, :vcenter)
point(x1, y1, " 丙円:r3,(x1,y1,r3)", :black, :left, :vcenter)
point(b, y, " B(b,y)", :black, :left, :bottom, delta=delta/2)
point(a, 0, " A(a,0)", :black, :left, :bottom, delta=delta/2)
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
