算額（その491）

算額（その491）

新潟県長岡市七日市 諏訪神社 嘉永2年(18499)
http://www.wasan.jp/niigata/niigatasuwa.html
みしま観光協会:和算の里みしま,2020/3
https://e-mishima.info/wp-content/uploads/2020/02/74a075f30cb806eef4a9a25c4a71608e.pdf

鈎股弦内に正方形が入っている。正方形の3つの頂点は三辺の上にある。鈎，股の長さが与えられたとき，正方形の面積が最小になるのはどのようなときか。

鈎股弦の辺の長さを鈎，股とする。
正方形の頂点の座標を下，右，上，左の順に (x1, 0), (0, y2), (x3, y3), (x1 + x3, y3 - y2) として以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms x1::positive, y2::positive, x3::positive, y3::positive,
     a, 鈎::positive, 股::positive, 弦::positive;
eq1 = x1^2 + y2^2 - a^2
eq2 = x3^2 + (y3 - y2)^2 - a^2
eq3 = -y2/x1 * (y3 - y2)/x3 + 1
eq4 = (鈎 - y3)/x3 - 鈎//股
eq5 = y3/(股 - x3) - 鈎//股
res = solve([eq1, eq2, eq3, eq5], (a, y2, x3, y3))

   4-element Vector{NTuple{4, Sym}}:
    (-sqrt(2*x1^2*股^2 - 2*x1^2*股*鈎 + x1^2*鈎^2 + 2*x1*股^2*鈎 + 股^2*鈎^2)/(股 - 鈎), 股*(x1 + 鈎)/(股 - 鈎), -股*(x1 + 鈎)/(股 - 鈎), 鈎*(x1 + 股)/(股 - 鈎))
    (sqrt(2*x1^2*股^2 - 2*x1^2*股*鈎 + x1^2*鈎^2 + 2*x1*股^2*鈎 + 股^2*鈎^2)/(股 - 鈎), 股*(x1 + 鈎)/(股 - 鈎), -股*(x1 + 鈎)/(股 - 鈎), 鈎*(x1 + 股)/(股 - 鈎))
    (-sqrt(2*x1^2*股^2 + 2*x1^2*股*鈎 + x1^2*鈎^2 - 2*x1*股^2*鈎 + 股^2*鈎^2)/(股 + 鈎), 股*(-x1 + 鈎)/(股 + 鈎), 股*(-x1 + 鈎)/(股 + 鈎), 鈎*(x1 + 股)/(股 + 鈎))
    (sqrt(2*x1^2*股^2 + 2*x1^2*股*鈎 + x1^2*鈎^2 - 2*x1*股^2*鈎 + 股^2*鈎^2)/(股 + 鈎), 股*(-x1 + 鈎)/(股 + 鈎), 股*(-x1 + 鈎)/(股 + 鈎), 鈎*(x1 + 股)/(股 + 鈎))

解は 4 組得られるが，最後のものが適解である。
a, y2, x3, y3 は x1, 鈎, 股 を含む式で表される。鈎，股は対象とする鈎股弦により決まるので，本質的に未知数なのは x1 である。
a = sqrt(2*x1^2*股^2 + 2*x1^2*股*鈎 + x1^2*鈎^2 - 2*x1*股^2*鈎 + 股^2*鈎^2)/(股 + 鈎) である。
(鈎, 股) = (3, 4) のとき，x1 を横軸，a = f(x1) を縦軸にしてグラフを描くと，x1 が 0.75 前後のときに a が最小になるのがわかる。

(鈎, 股) = (3, 4)
f = sqrt(2*x1^2*股^2 + 2*x1^2*股*鈎 + x1^2*鈎^2 - 2*x1*股^2*鈎 + 股^2*鈎^2)/(股 + 鈎)
plot(f, xlims=(0, 2), xlabel="x1", ylabel="f")

res2 = solve(diff(f, x1))[1]
res2 |> N |>  println
res2.evalf() |> println

   48//65
   0.738461538461539

x1 = 48//65 のときの a を求める。

f(x1 => 48//65) |> println
f(x1 => 48//65).evalf() |> println

   12*sqrt(65)/65
   1.48841681507050

鈎 = 3，股 = 4 のときには，x1 = 0.738462;  y2 = 1.29231;  x3 = 1.29231;  y3 = 2.03077 のとき，正方形の一辺の長さが最小値 12*sqrt(65)/65 = 1.48841681507050 になる。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股) = (3, 4)
   x1 = res2
   (a, y2, x3, y3) = (sqrt(2*x1^2*股^2 + 2*x1^2*股*鈎 + x1^2*鈎^2 - 2*x1*股^2*鈎 + 股^2*鈎^2)/(股 + 鈎), 股*(-x1 + 鈎)/(股 + 鈎), 股*(-x1 + 鈎)/(股 + 鈎), 鈎*(x1 + 股)/(股 + 鈎))
   @printf("a = %g;  x1 = %g;  y2 = %g;  x3 = %g;  y3 = %g\n", a, x1, y2, x3, y3)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:red, lw=0.5)
   plot!([0, x1, x1 + x3, x3, 0], [y2, 0, y3 - y2, y3, y2], color=:blue, lw=0.5)
   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(x1, 0, "x1  ", :blue, :right, :bottom, delta=delta)
       point(x1 + x3, y3 - y2, " (x1+x3,y3-y2)", :blue, :left, :vcenter)
       point(x3, y3, " (x3,y3)", :blue, :left, :vcenter)
       point(0, y2, "  y2", :blue, :left, :vcenter)
       point(0, 鈎, "  鈎", :red, :left, :vcenter)
       point(股, 0, "股", :red, :center, :bottom, delta=2delta)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その490）

算額（その490）

新潟県長岡市上岩井 根立寺 嘉永2年(1849)
http://www.wasan.jp/niigata/konryuji.html
涌田和芳，外川一仁: 三島根立寺の算額，長岡工業高等専門学校研究紀要，第53巻(2017)
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_53/53_17wakuta.pdf

正方形内に 2 本の斜線で分割された領域に大円 1 個，小円 3 個が入っている。
大円の直径が 10 寸のとき，小円の直径はいかほどか。

正方形の左下隅を原点とし，正方形の一辺の長さを a とする
大円の半径と中心座標を r1, (x1, r1)
小円の半径と中心座標を r2, (r2, a - r2), (r2, a - 3r2), (a - r2, a - r2)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r1::positive,
     x1::positive, r2::positive;

eq1 = distance(0, 0, b, a, r2, a - 3r2) - r2^2
eq2 = distance(0, 0, b, a, x1, r1) - r1^2
eq3 = distance(a, 0, b, a, a - r2, a - r2) - r2^2
eq4 = distance(a, 0, b, a, x1, r1) - r1^2;
# res = solve([eq1, eq2, eq3, eq4], (b, r1, x1, r2))

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)
   (b, r1, x1, r2) = u
   return [
       -r2^2 + (-b*(a^2 - 3*a*r2 + b*r2)/(a^2 + b^2) + r2)^2 + (a - a*(a^2 - 3*a*r2 + b*r2)/(a^2 + b^2) - 3*r2)^2,  # eq1
       -r1^2 + (-a*(a*r1 + b*x1)/(a^2 + b^2) + r1)^2 + (-b*(a*r1 + b*x1)/(a^2 + b^2) + x1)^2,  # eq2
       -r2^2 + (a - r2 - (a^3 - a^2*b + a*b^2 + a*b*r2 - b^2*r2)/(2*a^2 - 2*a*b + b^2))^2 + (-a*(a^2 - b*r2)/(2*a^2 - 2*a*b + b^2) + a - r2)^2,  # eq3
       -r1^2 + (x1 - (a*(a^2 - a*r1 - a*x1 + b*r1) + x1*(2*a^2 - 2*a*b + b^2))/(2*a^2 - 2*a*b + b^2))^2 + (-a*(a^2 - a*b + a*r1 - a*x1 + b*x1)/(2*a^2 - 2*a*b + b^2) + r1)^2,  # eq4
   ]
end;
a = 1
iniv = BigFloat[0.61, 0.3, 0.57, 0.14]
res = nls(H, ini=iniv)

   (BigFloat[0.6234898018587335305250048840042398106322575229856815465655969032398117579169351, 0.3079785283699041303721851029979308598027401532891747156054051562155987534336586, 0.5549581320873711914221948710064104810672818660788279185605428752567220611398606, 0.1539892641849520651860925514989654299013683741434972076685407572160588959569412], true)

   a = 1;  b = 0.62349;  r1 = 0.307979;  x1 = 0.554958;  r2 = 0.153989

小円の直径は大円の直径の 1/2 である。
大円の直径が 10 寸のとき，小円の直径は 5 寸である。

res[1][4]/res[1][2]  # r2/r1

   0.4999999999999999999999999999999999999999944720136849757617152696086889978777665

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1
   (b, r1, x1, r2) = res[1]
   @printf("a = %g;  b = %g;  r1 = %g;  x1 = %g;  r2 = %g\n", a, b, r1, x1, r2)
   @printf("大円の直径が 10 寸のとき，小円の直径は %g 寸である\n", 10r2/r1)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
   circle(x1, r1, r1, :red)
   circle(r2, a - 3r2, r2, :green)
   circle(r2, a - r2, r2, :green)
   circle(a - r2, a - r2, r2, :green)
   plot!([0, b, a], [0, a, 0], color=:blue, lw=0.5)
   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(x1, r1, "大円:r1,(x1,r1)", :red, :center, delta=-delta)
       point(r2, a - r2, "小円:r2,(r2,a-r2)", :green, :center, delta=-delta)
       point(r2, a - 3r2, "小円:r2,(r2,a-3r2)", :green, :center, delta=-delta)
       point(a - r2, a - r2, "小円:r2,(a-r2,a-r2)", :green, :center, delta=-delta)
       point(b, a, "(b,a)", :black, :center, :bottom, delta=delta/2)
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
       point(0, a, " a", :black, :left, :bottom, delta=delta/2)
   else
       plot!(showaxis=false)
   end
end;