算額（その376）

算額（その376）

東京都 住吉神社 嘉永8年(1851)
http://www.wasan.jp/tokyo/sumiyosi.html

厥円（弓型）内に，互いに外接する 3 個の等円が入っている。また，それぞれは外円に内接している。
矢と弦それぞれの長さが与えられたとき，等円の直径を求めよ。

過去に，算額（その11） で掲載したが，それは「矢と弦それぞれの長さが与えられたとき」という条件を無視したものであった。また，この問題においては，区切られた斜線に内接しているのではなく，互いに外接しているようである。
https://blog.goo.ne.jp/r-de-r/e/0dac8f68b8bbc6af25eb9e60427a9c59
今回新たに解を求める。

厥円のもとである外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (0, R - r) および (x, r + a)
とする。
直線と外円の交点座標を (x0, y0), y0 = sqrt(R^2 - x0^2)
とおき，以下の方程式を nlsolve() で解き，数値解を求める。

なお，問では「矢」と「弦」を与えるとあるので，方程式を解くに当たって，「矢」と「弦」から外円の径 R と切り取り位置 a（図参照）を求める式を求めておく。

using SymPy
@syms 矢::positive, 弦::positive, R::positive,  a::positive;

eq1 = R - a - 矢
eq2 = 2sqrt(R^2 - a^2) - 弦
solve([eq1, eq2], (R, a))

   1-element Vector{Tuple{Sym, Sym}}:
    ((弦^2 + 4*矢^2)/(8*矢), (弦 - 2*矢)*(弦 + 2*矢)/(8*矢))

include("julia-source.txt");

using SymPy
@syms 矢::positive, 弦::positive,
     R::positive,  a::positive, b::positive,
     x::positive, r::positibe;

eq1 = x^2 + (r + a)^2 - (R - r)^2
eq2 = x^2 + ((R - r) - (r + a))^2 - 4r^2

res = solve([eq1, eq2], (r, x))

   2-element Vector{Tuple{Sym, Sym}}:
    (R*(-R + a)/(-3*R + a), -sqrt(-1/(-3*R + a))*(-R + a)*sqrt(R + a))
    (R*(-R + a)/(-3*R + a), sqrt(-1/(-3*R + a))*(-R + a)*sqrt(R + a))

どちらも実質的に同じである。左右の円のどちらかを与えることになっている。

等円の半径は R*(-R + a)/(-3*R + a) であるが，これを 矢と弦を使って表そう。

eq = R*(-R + a)/(-3*R + a);

eq の R, a に以下を代入する。
R = (弦^2 + 4*矢^2)/(8*矢)
a = (弦 - 2*矢)*(弦 + 2*矢)/(8*矢)

eq(R => (弦^2 + 4*矢^2)/(8*矢), a => (弦 - 2*矢)*(弦 + 2*矢)/(8*矢)) |> simplify |> println

   矢*(弦^2 + 4*矢^2)/(2*(弦^2 + 8*矢^2))

等円の直径はこれを 2 倍して，矢*(弦^2 + 4*矢^2)/(弦^2 + 8*矢^2) である。

(矢, 弦) = (3, 8)
矢*(弦^2 + 4*矢^2)/(2*(弦^2 + 8*矢^2))

   2.2058823529411766

関数定義
等円直径(矢, 弦) = 矢*(弦^2 + 4*矢^2)/(弦^2 + 8*矢^2);

等円直径(3, 8)

   2.2058823529411766

等円直径(2, 10)

   1.7575757575757576

矢 = 3;  弦 = 8 のとき

   (25//6, 7//6)
   矢 = R - a = 3;  弦 = 2b = 8;  R = 4.16667;  a = 1.16667;  b = 4
   r = 1.10294;  x = 2.05798

using Plots

function draw(矢, 弦, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, a) = (弦^2//(8*矢) + 矢//2, 弦^2//(8*矢) - 矢//2)
   println((R, a))
   b = sqrt(R^2 - a^2)
   θ = round(Int64, atand(a, b))
   (r, x) = (R*(-R + a)/(-3*R + a), -sqrt(-1/(-3*R + a))*(-R + a)*sqrt(R + a))
   @printf("矢 = R - a = %g;  弦 = 2b = %g;  R = %g;  a = %g;  b = %g\n", 矢, 弦, R, a, b)
   @printf("r = %g;  x = %g\n", r, x)
   plot()
   circle(0, 0, R, :black, beginangle=θ, endangle=180 - θ)
   circle(0, R - r, r, :blue)
   circle(x, r + a, r, :blue)
   circle(-x, r + a, r, :blue)
   segment(-b, a, b, a, :black)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, a, " a", :black, delta=-delta/2)
       point(b, a, "(b,a)", :black, :right, delta=-delta/2)
       point(0, R - r, " R - r", :blue)
       point(0, R, " R", :black, :left, :bottom)
       point(x, r + a, "(x,r+a)", :blue, :center, delta=-delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

矢 = 2，弦 = 10 のとき

   (29//4, 21//4)
   矢 = R - a = 2;  弦 = 2b = 10;  R = 7.25;  a = 5.25;  b = 5
   r = 0.878788;  x = 1.74078

算額（その375）

算額（その375）

山形県山形市山寺 立石寺根本中堂
山形算額勝負-湯殿山神社を目指せ-
https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

厥円（弓型）内の 3 個の等円が 2 本の直線で隔てられている。
矢（R - a）と弦（2b）それぞれの長さが与えられたとき，等円の直径を求めよ。

過去に，算額（その11）で掲載したが，それは「矢と弦それぞれの長さが与えられたとき」という条件を無視したものであった。
https://blog.goo.ne.jp/r-de-r/e/0dac8f68b8bbc6af25eb9e60427a9c59
今回新たに解を求める。

厥円のもとである外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (0, R - r) および (x, r + a)
とする。
直線と外円の交点座標を (x0, y0), y0 = sqrt(R^2 - x0^2)
とおき，以下の方程式を nlsolve() で解き，数値解を求める。

なお，問では「矢」と「弦」を与えるとあるので，方程式を解くに当たって，「矢」と「弦」から外円の径 R と切り取り位置 a（図参照）を求める式を求めておく。

using SymPy
@syms 矢::positive, 弦::positive, R::positive,  a::positive;

eq1 = R - a - 矢
eq2 = 2sqrt(R^2 - a^2) - 弦
solve([eq1, eq2], (R, a))

   1-element Vector{Tuple{Sym, Sym}}:
    ((弦^2 + 4*矢^2)/(8*矢), (弦 - 2*矢)*(弦 + 2*矢)/(8*矢))

include("julia-source.txt");

using SymPy
@syms 矢::positive, 弦::positive,
     R::positive,  a::positive, b::positive, y0::positive,
     x::positive, r::positibe, x0::positive;

y0 = sqrt(R^2 - x0^2)
eq1 = x^2 + (r + a)^2 - (R - r)^2
eq2 = distance(0, a, x0, y0, x, r + a) - r^2
eq3 = distance(0, a, x0, y0, 0, R - r) - r^2;

# solve([eq1, eq2, eq3], (r, x, x0))

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, x, x0) = u
   return [
       x^2 - (R - r)^2 + (a + r)^2,  # eq1
       -r^2 + (x - x0*(-a*r + r*sqrt(R^2 - x0^2) + x*x0)/(R^2 + a^2 - 2*a*sqrt(R^2 - x0^2)))^2 + (a + r - (-x0*(a*x + r*x0 - x*sqrt(R^2 - x0^2)) + (a + r)*(R^2 + a^2 - 2*a*sqrt(R^2 - x0^2)))/(R^2 + a^2 - 2*a*sqrt(R^2 - x0^2)))^2,  # eq2
       -r^2 + x0^2*(a - sqrt(R^2 - x0^2))^2*(-R + a + r)^2/(R^2 + a^2 - 2*a*sqrt(R^2 - x0^2))^2 + (R - r - (x0^2*(-R + a + r) + (R - r)*(R^2 + a^2 - 2*a*sqrt(R^2 - x0^2)))/(R^2 + a^2 - 2*a*sqrt(R^2 - x0^2)))^2,  # eq3        
   ]
end;

(矢, 弦) = (3, 8)
(R, a) = (弦^2//(8*矢) + 矢//2, 弦^2//(8*矢) - 矢//2)
iniv = [big"1.0", 2, 2]
res = nls(H, ini=iniv)
println(Float64.(res[1]), ", converge = ",res[2])

   [1.0909090909090908, 2.088931871468374, 1.8031774011398194], converge = true

   矢 = R - a = 3;  弦 = 2b = 8;  R = 4.16667;  a = 1.16667;  b = 4
   r = 1.09091;  x = 2.08893;  x0 = 1.80318;  y0 = 3.75628

(矢, 弦) = (3, 8) のとき，等円の半径は 12/11 ≒ 1.090909 である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (矢, 弦) = (3, 8)
   (R, a) = (弦^2//(8*矢) + 矢//2, 弦^2//(8*矢) - 矢//2)
   println((R, a))
   b = sqrt(R^2 - a^2)
   θ = round(Int64, atand(a, b))
   (r, x, x0) = res[1]
   y0 = sqrt(R^2 - x0^2)
   @printf("矢 = R - a = %g;  弦 = 2b = %g;  R = %g;  a = %g;  b = %g\n", 矢, 弦, R, a, b)
   @printf("r = %g;  x = %g;  x0 = %g;  y0 = %g\n", r, x, x0, y0)
   plot()
   circle(0, 0, R, :black, beginangle=θ, endangle=180 - θ)
   segment(-b, a, b, a, :black)
   circle(0, R - r, r, :blue)
   circle(x, r + a, r, :blue)
   circle(-x, r + a, r, :blue)
   segment(0, a, x0, y0, :black)
   segment(0, a, -x0, y0, :black)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(x0, y0, " (x0,y0)", :black, :left, :bottom)
       point(0, a, " a", :black, delta=-delta/2)
       point(b, a, "(b,a)", :black, :right, delta=-delta/2)
       point(0, R - r, " R - r", :blue)
       point(0, R, " R", :black, :left, :bottom)
       point(x, r + a, "(x,r+a)", :blue, :center, delta=-delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;