算額（その345）

算額（その345）

美濃・飛騨の国の和算の歴史 -算額の問題に挑戦しよう-
http://hamaguri.sakura.ne.jp/minohidawasan.htm

岐阜県郡上市八幡町 郡上八幡神社（小野天満宮） 嘉永3年(1850)

鉤股弦に内接する円がある。鉤 + 股 = 17，(鉤^3 + 股^3 + 弦^3)/(弦 - 円の直径) = 213 のとき，弦の長さはいかほどか。

一般解を求めるために，与えられる 2 つの定数を a, b とおく。
鉤 + 股 = a
(鉤^3 + 股^3 + 円の直径^3)/(弦 - 円の直径) = b
eq3 は鉤股弦に内接する円の直径，eq4 は鉤股弦に関するピタゴラスの定理

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, 鉤::positive, 股::positive, 弦::positive, 円の直径::positive

(a, b) = (17, 213)
eq1 = 鉤 + 股 - a
eq2 = (鉤^3 + 股^3 + 円の直径^3)/(弦 - 円の直径) - b
eq3 = 鉤 + 股 - 弦 - 円の直径
eq4 = 鉤^2 + 股^2 - 弦^2
solve([eq1, eq2, eq3, eq4], (鉤, 股, 弦, 円の直径))

   4-element Vector{NTuple{4, Sym}}:
    (a/2 - sqrt(13*a^2 - 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, a/2 + sqrt(13*a^2 - 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, 2*a - sqrt(3*a^2 - 2*b), -a + sqrt(3*a^2 - 2*b))
    (a/2 + sqrt(13*a^2 - 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, a/2 - sqrt(13*a^2 - 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, 2*a - sqrt(3*a^2 - 2*b), -a + sqrt(3*a^2 - 2*b))
    (a/2 - sqrt(13*a^2 + 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, a/2 + sqrt(13*a^2 + 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, 2*a + sqrt(3*a^2 - 2*b), -a - sqrt(3*a^2 - 2*b))
    (a/2 + sqrt(13*a^2 + 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, a/2 - sqrt(13*a^2 + 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, 2*a + sqrt(3*a^2 - 2*b), -a - sqrt(3*a^2 - 2*b))

鉤 < 股 なので，最初の組のものが適解。
弦は 2*a - sqrt(3*a^2 - 2*b)
a = 17, b = 213 のとき，2*a - sqrt(3*a^2 - 2*b) = 13

(a, b) = (17, 213)
a/2 - sqrt(13*a^2 - 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2,
a/2 + sqrt(13*a^2 - 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2,
2*a - sqrt(3*a^2 - 2*b),
-a + sqrt(3*a^2 - 2*b)

   (5.0, 12.0, 13.0, 4.0)

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鉤, 股, 弦, 円の直径) = (5.0, 12.0, 13.0, 4.0)
   @printf("鉤 = %g, 股 = %g, 弦 = %g, 円の直径 = %g\n", 鉤, 股, 弦, 円の直径)
   plot([0, 股, 0, 0], [0, 0, 鉤, 0], color=:black, lw=0.5)
   circle(円の直径/2, 円の直径/2, 円の直径/2, :blue)
   if more
       point(円の直径/2, 円の直径/2, "", :blue)
       point(股, 0, " 股", :black, :left, :bottom)
       point(0, 鉤, " 鉤", :black, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その344）

算額（その344）

美濃・飛騨の国の和算の歴史 -算額の問題に挑戦しよう-
http://hamaguri.sakura.ne.jp/minohidawasan.htm
米山忠興:等円術(II) : 三斜等円術(一般公式)
https://toyo.repo.nii.ac.jp/?action=repository_uri&item_id=2534&file_id=18&file_no=1

岐阜県郡上市八幡町 郡上八幡神社（小野天満宮） 嘉永3年(1850)

三角形内に界斜を隔てて二等円を入れる。中斜は 257寸、小斜は 68寸、界斜は 40寸である。
大斜はいかほどか。

等円の半径を r, 大斜を x とする他，図のように記号を定め，以下の連立方程式を解く。

solve() で代数解を求めるのは困難なので，nlsolve() で数値解を求める。

include("julia-source.txt");

using SymPy

@syms r, x0, y0, x, x1, x2, x3, dummy

(中斜, 小斜, 界斜) = (257, 68, 40)
eq1 = (小斜 + 界斜 + x3)r +(中斜 + 界斜 +(x - x3))r - x*y0
eq2 = x0^2 + y0^2 - 小斜^2
eq3 = (x - x0)^2 + y0^2 - 中斜^2
eq4 = distance(0, 0, x0, y0, x1, r) - r^2
eq5 = distance(x0, y0, x3, 0, x1, r) - r^2
eq6 = distance(x0, y0, x3, 0, x2, r) - r^2
eq7 = distance(x0, y0, x, 0, x2, r) - r^2;
eq4 = apart(eq4, dummy) |> numerator
eq5 = apart(eq5, dummy) |> numerator
eq6 = apart(eq6, dummy) |> numerator
eq7 = apart(eq7, dummy) |> numerator;

# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7])

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
println(eq5, ",  # eq5")
println(eq6, ",  # eq6")
println(eq7, ",  # eq7")

   r*(x3 + 108) + r*(x - x3 + 297) - x*y0,  # eq1
   x0^2 + y0^2 - 4624,  # eq2
   y0^2 + (x - x0)^2 - 66049,  # eq3
   -r^2*y0^2 - 2*r*x0*x1*y0 + x1^2*y0^2,  # eq4
   -r^2*y0^2 - 2*r*x0*x1*y0 + 2*r*x0*x3*y0 + 2*r*x1*x3*y0 - 2*r*x3^2*y0 + x1^2*y0^2 - 2*x1*x3*y0^2 + x3^2*y0^2,  # eq5
   -r^2*y0^2 - 2*r*x0*x2*y0 + 2*r*x0*x3*y0 + 2*r*x2*x3*y0 - 2*r*x3^2*y0 + x2^2*y0^2 - 2*x2*x3*y0^2 + x3^2*y0^2,  # eq6
   -r^2*y0^2 - 2*r*x^2*y0 + 2*r*x*x0*y0 + 2*r*x*x2*y0 - 2*r*x0*x2*y0 + x^2*y0^2 - 2*x*x2*y0^2 + x2^2*y0^2,  # 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=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r, x0, y0, x, x1, x2, x3) = u
   return [
       r*(x3 + 108) + r*(x - x3 + 297) - x*y0,  # eq1
       x0^2 + y0^2 - 4624,  # eq2
       y0^2 + (x - x0)^2 - 66049,  # eq3
       -r^2*y0^2 - 2*r*x0*x1*y0 + x1^2*y0^2,  # eq4
       -r^2*y0^2 - 2*r*x0*x1*y0 + 2*r*x0*x3*y0 + 2*r*x1*x3*y0 - 2*r*x3^2*y0 + x1^2*y0^2 - 2*x1*x3*y0^2 + x3^2*y0^2,  # eq5
       -r^2*y0^2 - 2*r*x0*x2*y0 + 2*r*x0*x3*y0 + 2*r*x2*x3*y0 - 2*r*x3^2*y0 + x2^2*y0^2 - 2*x2*x3*y0^2 + x3^2*y0^2,  # eq6
       -r^2*y0^2 - 2*r*x^2*y0 + 2*r*x*x0*y0 + 2*r*x*x2*y0 - 2*r*x0*x2*y0 + x^2*y0^2 - 2*x*x2*y0^2 + x2^2*y0^2,  # eq7
   ]
end;

iniv = [big"34.98", 67.27, 117.06, 291.98, 60.54, 119.75, 141.28]
res = nls(H, ini=iniv);
println(res);

 (BigFloat[14.00000000000000000002832495142943662706356368549023146540001887379607263571971, 59.99999999999999999997271770164078264394833470869088845531230198310065464574821, 32.0000000000000000000659299852261495342990267083663162723582624564073423302806, 314.9999999999999999999660433365560846765482392433240243631619025246373680740372, 55.99999999999999999997472857143202733080502439021613717408206542155355708244299, 90.99999999999999999998943873473162951157245262928774364469574387503292726000683, 83.99999999999999999998599488059260011385774611340144279061843511421554191942925], true)

大斜は 315 寸である。
図を見るとわかるが，この三角形は随分と扁平で，算額の図とはイメージが異なる（しかし，これが実際の姿である）。

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x0, y0, x, x1, x2, x3) = res[1]
   @printf("r = %g, x0 = %g, y0 = %g, x = %g, x1 = %g, x2 = %g, x3 = %g\n",
       r, x0, y0, x, x1, x2, x3)
   (中斜, 小斜, 界斜) = (257, 68, 40)
   @printf("大斜 = %g;  術による答え = %g\n", x, sqrt((中斜 + 小斜)^2 - 4界斜^2))
   plot([0, x, x0, 0], [0, 0, y0, 0], color=:black, lw=0.5)
   segment(x0, y0, x3, 0, :orange)
   circle(x1, r, r, :blue)
   circle(x2, r, r, :blue)
   if more
       point(x0, y0, "(x0,y0)", :black, :center, :bottom)
       point(x, 0, " x", :black, :left, :bottom)
       point(x3, 0, "x3", :black, :center, :top)
       point(x1, r, "(x1,r) ", :blue, :right, :vcenter)
       point(x2, r, " (x2,r)", :blue, :left, :vcenter)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;