算額（その406）

算額（その406）

兵庫県伊丹市 昆陽寺 嘉永7年(1854)，昭和43年(1968)復元奉納
https://toyo.repo.nii.ac.jp/?action=repository_action_common_download&item_id=2534&item_no=1&attribute_id=18&file_no=1

第7問，三角形内に界斜を引き，2 区分された領域に等円を入れる。大斜，中斜，小斜が与えられたとき，界斜はいかほどか。

大斜，中斜，小斜，界斜を変数名として，
等円の半径と中心座標を r, (x1, r), (x2, r)
三角形の頂点座標を (x, y)
界斜と大斜の交点座標を (a, 0) として，以下の連立方程式を立て，nlsolve() により数値解を求める。

include("julia-source.txt");

using SymPy

@syms r::positive, a::positive, x::positive, y::positive, x1::positive,
     x2::positive, 大斜::positive, 中斜::positive, 小斜::positive, 界斜::positive;

eq1 = (界斜 + 中斜 + a)r + (界斜 + 小斜 + 大斜 - a)r - 大斜*y
eq2 = x^2 + y^2 - 中斜^2
eq3 = (大斜 - x)^2 + y^2 - 小斜^2
eq4 = distance(a, 0, x, y, x1, r) - r^2
eq5 = distance(a, 0, x, y, x2, r) - r^2
eq6 = distance(0, 0, x, y, x1, r) - r^2
eq7 = distance(大斜, 0, x, y, x2, r) - r^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7])

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)
   (界斜, r, a, x, y, x1, x2) = u
   return [
       r*(a + 中斜 + 界斜) + r*(-a + 大斜 + 小斜 + 界斜) - y*大斜,  # eq1
       x^2 + y^2 - 中斜^2,  # eq2
       y^2 - 小斜^2 + (-x + 大斜)^2,  # eq3
       -r^2 + (r - y*(a^2 - a*x - a*x1 + r*y + x*x1)/(a^2 - 2*a*x + x^2 + y^2))^2 + (x1 - (a^2*x1 - a*r*y - 2*a*x*x1 + a*y^2 + r*x*y + x^2*x1)/(a^2 - 2*a*x + x^2 + y^2))^2,  # eq4
       -r^2 + (r - y*(a^2 - a*x - a*x2 + r*y + x*x2)/(a^2 - 2*a*x + x^2 + y^2))^2 + (x2 - (a^2*x2 - a*r*y - 2*a*x*x2 + a*y^2 + r*x*y + x^2*x2)/(a^2 - 2*a*x + x^2 + y^2))^2,  # eq5
       -r^2 + (r - y*(r*y + x*x1)/(x^2 + y^2))^2 + (-x*(r*y + x*x1)/(x^2 + y^2) + x1)^2,  # eq6
       -r^2 + (r - y*(r*y + x*x2 - x*大斜 - x2*大斜 + 大斜^2)/(x^2 - 2*x*大斜 + y^2 + 大斜^2))^2 + (x2 - (x2*(x^2 - 2*x*大斜 + y^2 + 大斜^2) + y*(r*x - r*大斜 - x2*y + y*大斜))/(x^2 - 2*x*大斜 + y^2 + 大斜^2))^2,  # eq7   
   ]
end;

(大斜, 中斜, 小斜) = (9//2, 4, 3)
iniv = [big"3.0", 1.1, 2.5, 3, 3, 2, 4]
res = nls(H, ini=iniv);
println(res);
   (BigFloat[2.680951323690902076203539918038571176272552963887644954668323575268635315149453, 0.6975859391932934831098526753294690788534179473670310788732592770458445561358907, 2.43201081695757731639921335154698418305054378580278996601587232707696243601583, 3.027777777777777777777777777777777777777777777777777777777777777777777777777774, 2.613916932191048363498370393355647818939596565227754137625588266018177169955276, 1.875529746633337620097836716754206503388995410957580351169949221547302218391068, 3.306481070324239696301376634792777679661548374845176108710264222377909400696189], true)

大斜, 中斜, 小斜 がそれぞれ 4.5, 4, 3 のとき，
界斜 = 2.68095;  r = 0.697586;  a = 2.43201;  x = 3.02778;  y = 2.61392;  x1 = 1.87553;  x2 = 3.30648
である。
算額の術では 界斜 = sqrt((中斜 + 小斜)^2 - 大斜^2)/2 となっており，正しい答えが求まっていることがわかる。

(大斜, 中斜, 小斜) = (9//2, 4, 3)
sqrt((中斜 + 小斜)^2 - 大斜^2)/2

   2.680951323690902

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (大斜, 中斜, 小斜) = (9//2, 4, 3)
   (界斜, r, a, x, y, x1, x2) = res[1]
   @printf("界斜 = %g;  r = %g;  a = %g;  x = %g;  y = %g;  x1 = %g;  x2 = %g\n", 界斜, r, a, x, y, x1, x2)
   plot([0, 大斜, x, 0], [0, 0, y, 0], color=:black, lw=0.5)
   circle(x1, r, r)
   circle(x2, r, r)
   segment(a, 0, x, y, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(大斜, 0, " 大斜", :black, :left, :bottom, delta=delta)
       point(x, y, " (x,y)", :black, :left, :vcenter, delta=-delta)
       point(x1, r, "(x1,r)", :red, :center, delta=-delta)
       point(x2, r, "(x2,r)", :red, :center, delta=-delta)
       point(a, 0, " a", :black, :left, :bottom, delta=delta)
       point(x/2, 2y/3, "中斜", :black, :left, mark=false)
       point((大斜 + x)/2, 2y/3, "小斜", :black, :left, mark=false)
       point((a + x)/2, 2y/3, "  界斜", :blue, :left, mark=false)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;