算額（その430）
東京都浅草 池田邸内 伽山
中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html
大円の中に甲円 2 個と乙円 1 個，等円 3 個が入っている。等円は 1 本の元と長さの等しい 3 本の面（隣り合う 3 本の弦）に接している。面の長さが 1 寸のとき，乙円の直径を求めよ。
大円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (r1, y1) ただし y1 = a + r1
乙円の半径と中心座標を r2, (0, r0 - r2)
等円の半径と中心座標を r3, (0, a - r3), (2r3, a - r3)
弦と y 軸の交点座標を a
真ん中の面と y 軸の交点座標を b
面の長さの半分を c とし，以下の連立方程式の数値解を求める。
include("julia-source.txt");
using SymPy
@syms r0::positive, r1::positive, r2::positive,
r3::positive, a::negative;
c = 1//2
b = a - 2r3
d = sqrt(r0^2 - a^2)
y1 = a + r1
eq1 = r1^2 + (r0 - r2 - y1)^2 - (r1 + r2)^2
eq2 = r1^2 + y1^2 - (r0 - r1)^2
eq3 = r0^2 - b^2 - c^2
eq4 = distance(c, b, d, a, 2r3, b+ r3) - r3^2
eq5 = (d - c)^2 + (a - b)^2 - (2c)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r0, r1, r2, r3, a))
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)
(r0, r1, r2, r3, a) = u
return [
r1^2 - (r1 + r2)^2 + (-a + r0 - r1 - r2)^2, # eq1
r1^2 + (a + r1)^2 - (r0 - r1)^2, # eq2
r0^2 - (a - 2*r3)^2 - 1/4, # eq3
-r3^2 + (2*r3 - 2*r3*(-4*a^2 + 4*r0^2 + 4*r3*sqrt(-a^2 + r0^2) + 2*r3 - 4*sqrt(-a^2 + r0^2) + 1)/(-4*a^2 + 4*r0^2 + 16*r3^2 - 4*sqrt(-a^2 + r0^2) + 1))^2 + (a - r3 - (2*r3*(8*r3^2 + 8*r3*sqrt(-a^2 + r0^2) - 4*r3 - 2*sqrt(-a^2 + r0^2) + 1) + (a - 2*r3)*(-4*a^2 + 4*r0^2 + 16*r3^2 - 4*sqrt(-a^2 + r0^2) + 1))/(-4*a^2 + 4*r0^2 + 16*r3^2 - 4*sqrt(-a^2 + r0^2) + 1))^2, # eq4
4*r3^2 + (sqrt(-a^2 + r0^2) - 1/2)^2 - 1, # eq5
]
end;
iniv = [big"65.0", 31, 27, 10, -32] ./ 40
res = nls(H, ini=iniv);
println(res);
(BigFloat[2.176250899482821511100052865997767880198073191589329947230101745924818974139391, 0.9777711616568012201157129796477534133807314665018133584786185122552240330104753, 1.070019583175477769385089300350675253824657209000148549887247280185448490007896, 0.2236067977499789696409173668731276235440618359611525724270897245410529337440141, -1.670820393249936908922752100619382870632185507883457717281269173623140154153112], true)
r0 = 2.17625; r1 = 0.977771; , r2 = 1.07002; , r3 = 0.223607; , a = -1.67082
乙円の直径 = 2.14004; 甲円の直径 = 1.95554; 等円の直径 = 0.447214
面 = 1; 弦 = 2d = 2.78885; 矢 = r0 + a = 0.505431
面の長さが 1 寸のとき，乙円の直径は 2.14004 寸である。
using Plots
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 169.3/2
(r0, r1, r2, r3, a) = res[1]
c = 1//2
b = a - 2r3
d = sqrt(r0^2 - a^2)
y1 = a + r1
@printf("r0 = %g; r1 = %g; , r2 = %g; , r3 = %g; , a = %g\n", r0, r1, r2, r3, a)
@printf("乙円の直径 = %g; 甲円の直径 = %g; 等円の直径 = %g\n", 2r2, 2r1, 2r3)
@printf("面 = %g; 弦 = 2d = %g; 矢 = r0 + a = %g\n",
sqrt((sqrt(r0^2 - a^2) - c)^2 + (a - b)^2), 2d, r0 + a)
plot()
circle(0, 0, r0, :red)
circle(r1, y1, r1, :green)
circle(-r1, y1, r1, :green)
circle(0, r0 - r2, r2, :magenta)
circle(0, a - r3, r3, :orange)
circle(2r3, a - r3, r3, :orange)
circle(-2r3, a - r3, r3, :orange)
segment(-c, b, c, b, :blue)
segment(-d, a, d, a, :blue)
segment(c, b, d, a, :blue)
segment(-c, b, -d, a, :blue)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) / 3 # size[2] * fontsize * 2
point(r1, y1, "甲円:r1,(r1,y1)", :green, :center, :top, delta=-delta)
point(0, r0 - r2, " 乙円:r2,(0,r0-r2)", :magenta, :center, :top, delta=-delta)
point(0, a - r3, "a-r3 ", :black, :right, :vcenter)
point(2r3, a - r3, "(2r3,a-r3)", :black, :center, :bottom, delta=delta)
point(c, b, "(c,b)")
point(d, a, "(d,a)")
vline!([0], color=:black, lw=0.5)
hline!([0], color=:black, lw=0.5)
else
plot!(showaxis=false)
end
end;
術では，乙円径は 1.4221...寸としているが，この値は間違っている。
術の式通りに計算すると下記のように 2.14004...寸となる。
@syms 角, 元, 氏, 房, 面
面 = 1
角 = sqrt(80)
元 = sqrt(角 + 10)
氏 = 元 - sqrt(角 + 5.8 - sqrt(12.8))
房 = sqrt(2氏 * 元) - 氏
乙円径 = (元 * 房 * 面) / (2元 -(2sqrt((元 -2房)元) + 房))
@printf("面 = %g; 角 = %g; 元 = %g; 氏 = %g; 房 = %g\n", 面, 角, 元, 氏, 房)
@printf("乙円径 = %g\n", 乙円径)
面 = 1; 角 = 8.94427; 元 = 4.3525; 氏 = 1.01086; 房 = 1.95554
乙円径 = 2.14004
中村は，術と中村の導いた乙円径を表す式が同じであることを示しているが，(7)の分母の式において（またしても）符号を誤っているのでそれに基づいたのでは正しい答えにならない。"+ 房" は "- 房" である。
@printf("大円径 = %g\n", 元*面)
@printf("h(矢) = %g\n", 氏*面/2)
@printf("(7) の分子 = 甲円径 = %g\n", 房*面)
@printf("(7) の分母 = %g\n", (2元 - 2sqrt(元*(元 - 2房)) - 房)面)
x = (R*d/((2元 - 2sqrt(元*(元 - 2房)) - 房)面))
@printf("乙円径 = %g\n", x.evalf())
大円径 = 4.3525
h(矢) = 0.505431
(7) の分子 = 甲円径 = 1.95554
(7) の分母 = 3.97726
乙円径 = 2.14004
以下は，中村による現代解法である。
@syms R, # 大円径
r, # 等円径
d, # 甲円径
x, # 乙円径
h, # 矢
a # 面
a = 1
r = a/sqrt(Sym(5))
@printf("等円径 = %g\n", r.evalf())
R = sqrt(10 + 4sqrt(Sym(5)))a
@printf("大円径 = %g\n", R.evalf())
h = (sqrt(10 + 4sqrt(Sym(5))) - sqrt(4sqrt(Sym(5)) + 29//5 - 8sqrt(Sym(5))/5))a/2
@printf("h(矢) = %g\n", h.evalf())
eq3 = (d/2)^2 - h*(R - d) + h^2;
d = solve(eq3, d)[1]
@printf("甲円径 = %g\n", d.evalf())
x = (R*(2R - d) + 2R*sqrt(R*(R - 2d)))/(d + 4R)
@printf("乙円径 = %g\n", x.evalf())
等円径 = 0.447214
大円径 = 4.3525
h(矢) = 0.505431
甲円径 = 1.95554
乙円径 = 2.14004
#算額
