算額（その769）
神壁算法 天明5年乙巳 秋九月 第二術 關流四傅藤田権平定資六弟子 奧州一關 梶山主水次俊
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
大円，小円，甲円，乙円がある。上下の二線（上線，下線）と斜線を引く。
記述が重複するが正確を期すため，これらは以下のように接する。
大円は二線，斜線，小円，甲円，乙円に接する。
小円は二線，大円に接する。
甲円は上線，斜線，大円に接する。
乙円は下線，斜線，大円に接する。
大円，小円，甲円の直径がそれぞれ 245寸，196寸，25寸のとき，乙円の直径はいかほどか。
大円の半径と中心座標を r1, (x1, r1)
小円の半径と中心座標を r2, (0, r2)
甲円の半径と中心座標を r3, (x3, y3)
乙円の半径と中心座標を r4, (x4, r4)
上線は二点 (x02, y02), (x03, y03) を通る
下線は二点 (x02, y02), (x01, 0) を通る
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms r1::positive, x1::positive, r2::positive,
r3::positive, x3::positive, y3::positive,
r4::positive, x4::positive,
x01::positive, x02::positive, y02::positive,
x03::positive, y03::positive
eq1 = dist(x02, y02, x03, y03, x1, r1) - r1^2
eq2 = dist(x02, y02, x03, y03, 0, r2) - r2^2
eq3 = dist(x02, y02, x03, y03, x3, y3) - r3^2
eq4 = dist(x01, 0, x02, y02, x1, r1) - r1^2
eq5 = dist(x01, 0, x02, y02, x3, y3) - r3^2
eq6 = dist(x01, 0, x02, y02, x4, r4) - r4^2
eq7 = sqrt((x02 - x03)^2 + (y02 - y03)^2) - x1
eq8 = x1^2 + (r1 - r2)^2 - (r1 + r2)^2
eq9 = (x1 - x3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq10 = (x1 - x4)^2 + (r1 - r4)^2 - (r1 + r4)^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)
(x1, x3, y3, r4, x4, x01, x02, y02, x03, y03) = u
return [
-r1^2 + (r1 - y02 - (-y02 + y03)*((r1 - y02)*(-y02 + y03) + (-x02 + x03)*(-x02 + x1))/((-x02 + x03)^2 + (-y02 + y03)^2))^2 + (-x02 + x1 - (-x02 + x03)*((r1 - y02)*(-y02 + y03) + (-x02 + x03)*(-x02 + x1))/((-x02 + x03)^2 + (-y02 + y03)^2))^2, # eq1
-r2^2 + (-x02 - (-x02 + x03)*(-x02*(-x02 + x03) + (r2 - y02)*(-y02 + y03))/((-x02 + x03)^2 + (-y02 + y03)^2))^2 + (r2 - y02 - (-y02 + y03)*(-x02*(-x02 + x03) + (r2 - y02)*(-y02 + y03))/((-x02 + x03)^2 + (-y02 + y03)^2))^2, # eq2
-r3^2 + (-x02 + x3 - (-x02 + x03)*((-x02 + x03)*(-x02 + x3) + (-y02 + y03)*(-y02 + y3))/((-x02 + x03)^2 + (-y02 + y03)^2))^2 + (-y02 + y3 - (-y02 + y03)*((-x02 + x03)*(-x02 + x3) + (-y02 + y03)*(-y02 + y3))/((-x02 + x03)^2 + (-y02 + y03)^2))^2, # eq3
-r1^2 + (r1 - y02*(r1*y02 + (-x01 + x02)*(-x01 + x1))/(y02^2 + (-x01 + x02)^2))^2 + (-x01 + x1 - (-x01 + x02)*(r1*y02 + (-x01 + x02)*(-x01 + x1))/(y02^2 + (-x01 + x02)^2))^2, # eq4
-r3^2 + (-y02*(y02*y3 + (-x01 + x02)*(-x01 + x3))/(y02^2 + (-x01 + x02)^2) + y3)^2 + (-x01 + x3 - (-x01 + x02)*(y02*y3 + (-x01 + x02)*(-x01 + x3))/(y02^2 + (-x01 + x02)^2))^2, # eq5
-r4^2 + (r4 - y02*(r4*y02 + (-x01 + x02)*(-x01 + x4))/(y02^2 + (-x01 + x02)^2))^2 + (-x01 + x4 - (-x01 + x02)*(r4*y02 + (-x01 + x02)*(-x01 + x4))/(y02^2 + (-x01 + x02)^2))^2, # eq6
-x1 + sqrt((x02 - x03)^2 + (y02 - y03)^2), # eq7
x1^2 + (r1 - r2)^2 - (r1 + r2)^2, # eq8
(-r1 + y3)^2 - (r1 + r3)^2 + (x1 - x3)^2, # eq9
(r1 - r4)^2 - (r1 + r4)^2 + (x1 - x4)^2, # eq10
]
end;
(r1, r2, r3) = (245, 196, 25) .// 2
iniv = BigFloat[219.13, 118.51, 212.5, 24.5, 109.57, 82.18, 107.08, 222.73, -106.65, 174.33]
res = nls(H, ini=iniv)
([219.1346617949794, 118.51160280748886, 212.5, 24.5, 109.5673308974897, 82.17549817311728, 107.0771642861831, 222.72727272727272, -106.6467651187968, 174.33221099887766], true)
乙円の半径は 24.5寸 である（直径 49寸）。
その他のパラメータは以下のとおり。
r1 = 122.5; r2 = 98; r3 = 12.5; x1 = 219.135; x3 = 118.512; y3 = 212.5; r4 = 24.5; x4 = 109.567; x01 = 82.1755; x02 = 107.077; y02 = 222.727; x03 = -106.647; y03 = 174.332
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, r3) = (245, 196, 25) .// 2
(x1, x3, y3, r4, x4, x01, x02, y02, x03, y03) = res[1]
@printf("大円，小円，甲円の直径がそれぞれ %g寸, %g寸, %g寸 のとき，乙円の直径は %g寸\n", 2r1, 2r2, 2r3, 2r4)
@printf("r1 = %g; r2 = %g; r3 = %g; x1 = %g; x3 = %g; y3 = %g; r4 = %g; x4 = %g; x01 = %g; x02 = %g; y02 = %g; x03 = %g; y03 = %g\n", r1, r2, r3, x1, x3, y3, r4, x4, x01, x02, y02, x03, y03)
plot()
circle(x1, r1, r1)
circle(0, r2, r2, :blue)
circle(x3, y3, r3, :green)
circle(x4, r4, r4, :magenta)
segment(x01, 0, x02, y02)
y04 = (y02 - y03)/(x02 - x03)*(1.25x1 - x03) + y03
segment(x03, y03, 1.25x1, y04)
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(x01, 0, "x01 ", :black, :right, :bottom, delta=delta)
point(x02, y02, "(x02,y02)", :black, :right, :bottom, delta=delta)
point(x03, y03, "(x03,y03)", :black, :left, :bottom, delta=4delta)
point(x1, r1, " 大円:r1,(x1,r1)", :red, :center, delta=-delta)
point(0, r2, " 小円:r2,(0,r2)", :blue, :center, delta=-delta)
point(x3, y3, " 甲円:r3,(x3,y3)", :green, :left, :vcenter)
point(x4, r4, " 乙円:r4,(x4,r4)", :magenta, :left, :vcenter)
segment(x03, 0, 1.4x1, 0)
end
end;
