算額（その694）

算額（その694）

神壁算法 關流藤田貞資門人 久世大和守家士 平井彌五太夫正義 寛政八年
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

鈎股弦（直角三角形）内に，斜線，全円，甲円，乙円が入っている。全円，甲円の直径がそれぞれ 19寸7分5厘，3寸1分6厘 のとき，乙円の直径はいかほどか。

直角を挟む2辺（鈎，股）の長さを「鈎」，「股」
斜線と斜辺の交点座標を (x, y)
全円の半径と中心座標を r0, (r0, r0)
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (r2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r0, r1, 鈎, 股, x1, r2, y2, x, y;

eq1 = 鈎 + 股 - sqrt(鈎^2 + 股^2) - 2r0
eq2 = (r0 - x1)^2 + (r0 - r1)^2 - (r0 + r1)^2
eq3 = (r0 - r2)^2 + (r0 - y2)^2 - (r0 + r2)^2
eq4 = dist(0, 0, x, y, x1, r1) - r1^2
eq5 = dist(0, 0, x, y, r2, y2) - r2^2
eq6 = (鈎 + 股 + sqrt(鈎^2 + 股^2))*r0 - 鈎*股
eq7 = 鈎/股 - y/(股 - x);

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)
   (鈎, 股, x1, r2, y2, x, y) = u
   return [
       -2*r0 + 股 + 鈎 - sqrt(股^2 + 鈎^2),  # eq1
       (r0 - r1)^2 - (r0 + r1)^2 + (r0 - x1)^2,  # eq2
       (r0 - r2)^2 - (r0 + r2)^2 + (r0 - y2)^2,  # eq3
       -r1^2 + (r1 - y*(r1*y + x*x1)/(x^2 + y^2))^2 + (-x*(r1*y + x*x1)/(x^2 + y^2) + x1)^2,  # eq4
       -r2^2 + (r2 - x*(r2*x + y*y2)/(x^2 + y^2))^2 + (-y*(r2*x + y*y2)/(x^2 + y^2) + y2)^2,  # eq5
       r0*(股 + 鈎 + sqrt(股^2 + 鈎^2)) - 股*鈎,  # eq6
       -y/(-x + 股) + 鈎/股,  # eq7
   ]
end;

(r0, r1) = (1975, 316) .// 200
iniv = BigFloat[26, 49, 2, 0.6, 5, 5, 24]
res = nls(H, ini=iniv)

   (BigFloat[26.50958690810937781981564610070982992662591012710315095504511369567094160353881, 48.6025397559463070232560931975192400228754697146791710191805408688433795666828, 1.974999999999999999999999999999999999999999999999999999999999999999999999999969, 0.5700002001441297307780475616969606044325776049176168848094355451524229609489859, 5.130001801297167577002428055272645439893198444258551963284919906371806648542264, 5.312669489352098474921664382153576018135850078829535519263800791306165244941906, 23.61186439712043766631850836512700452504822257257571341895022573913851219974303], true)

乙円の直径 = 1.1400004002882596

その他のパラメータは以下のとおり。

鈎 = 26.5096;  股 = 48.6025;  x1 = 1.975;  r2 = 0.57;  y2 = 5.13;  x = 5.31267;  y = 23.6119

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1) = (1975, 316) .// 200
   (鈎, 股, x1, r2, y2, x, y) = res[1]
   println("乙円の直径 = $(Float64(2r2))")
   @printf("鈎 = %g;  股 = %g;  x1 = %g;  r2 = %g;  y2 = %g;  x = %g;  y = %g\n", 鈎, 股, x1, r2, y2, x, y)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(r0, r0, r0, :blue)
   circle(x1, r1, r1)
   circle(r2, y2, r2, :green)
   segment(0, 0, x, y, :orange)
   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(x, y, " (x,y)", :black, :left, :bottom)
       point(r0, r0, " 全円:r0\n (r0,r0)", :blue, :left, :vcenter)
       point(x1, r1, "  甲円:r1,(x1,r1)", :red, :left, :bottom, delta=delta/2)
       point(r2, y2, " 乙円:r2,(r2,y2)", :green, :left, :vcenter)
       point(股, 0, " 股", :black, :left, :bottom, delta=delta/2)
       point(0, 鈎, " 鈎", :black, :left, :bottom, delta=delta/2)
   end
end;