算額（その1509）

算額（その1509）

六十八 岩手県一関市川崎町薄衣諏訪前 浪分神社 文政5年(1822)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

今有如図 03004
https://w.atwiki.jp/sangaku/pages/202.html

外円の中に弦を 4 本引き，区画された領域に等円 5 個を容れる。等円の直径が与えられたとき，外円の直径を求める術を述べよ。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (0, R - r), (x1, y1), (x2, y2)
弦の端点座標を (x01, sqrt(R^2 - x01^2)), (x02, sqrt(R^2 - x02^2)), (0, -R)
とおき，以下の連立方程式を解く。

R, r1 を記号のままにして SymPy では簡単に解くことができない。数値を代入して解けば，数値解は求まる。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R::positive, r::positive, x1::positive, y1::positive,
      x2::positive, y2::negative,  x01::positive, x02::positive
eq1 = x1^2 + y1^2 - (R - r)^2
eq2 = x2^2 + y2^2 - (R - r)^2
eq3 = dist(0, -R, x01, sqrt(R^2 - x01^2), 0, R - r) - r^2
eq4 = dist(0, -R, x01, sqrt(R^2 - x01^2), x1, y1) - r^2
eq5 = dist(0, -R, x02, sqrt(R^2 - x02^2), x1, y1) - r^2
eq6 = dist(0, -R, x02, sqrt(R^2 - x02^2), x2, y2) - r^2
eq7 = (sqrt(R^2 - x02^2) + R)/x02 * y2/x2 + 1;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (R, x1, y1, x2, y2, x01, x02))  

function H(u)
    (R, x1, y1, x2, y2, x01, x02) = u
    return [
        x1^2 + y1^2 - (R - r)^2,
        x2^2 + y2^2 - (R - r)^2,
        dist(0, -R, x01, sqrt(R^2 - x01^2), 0, R - r) - r^2,
        dist(0, -R, x01, sqrt(R^2 - x01^2), x1, y1) - r^2,
        dist(0, -R, x02, sqrt(R^2 - x02^2), x1, y1) - r^2,
        dist(0, -R, x02, sqrt(R^2 - x02^2), x2, y2) - r^2,
        (sqrt(R^2 - x02^2) + R)/x02 * y2/x2 + 1
    ]
end;
r = 2
iniv = BigFloat[10, 5, 5, 7, -3, 3, 8]
res = nls(H, ini=iniv)

    ([7.509090269827044, 3.8101684797908897, 3.9790315098916964, 4.870544338179145, -2.5744656631880116, 2.2798713379711115, 6.204719439604394], true)

等円の直径が 4 のとき，外円の直径は 2*7.509090269827044 = 15.018180539654088 である。

function draw(r, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (R, x1, y1, x2, y2, x01, x02) = res[1]
    plot()
    circle(0, 0, R, :green)
    circle(0, R - r, r)
    circle2(x1, y1, r)
    circle2(x2, y2, r)
    plot!([-x01, 0, x01], [sqrt(R^2 - x01^2), -R, sqrt(R^2 - x01^2)], color=:blue, lw=0.5)
    plot!([-x02, 0, x02], [sqrt(R^2 - x02^2), -R, sqrt(R^2 - x02^2)], color=:blue, lw=0.5)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(0, R, "R", :green, :center, :bottom, delta=delta)
        point(x1, y1, "r,(x1,y1)", :red, :center, delta=-delta)
        point(x2, y2, "r,(x2,y2)", :red, :center, delta=-delta)
        point(x01, sqrt(R^2 - x01^2), "(x01,y01)", :blue, :left, :bottom, delta=delta)
        point(x02, sqrt(R^2 - x02^2), "(x02,y02)", :blue, :left, :bottom, delta=delta)
        point(0, R - r, " R-r", :blue, :left, :vcenter)
    end
end;

draw(2, true)

術は，「外円の直径を x，等円の直径を d として 以下の d の 5 次式を解くとしている。
-d^5 + 4*d^4*x - 14*d^3*x^2 + 52*d^2*x^3 - 73*d*x^4 + 16*x^5

using SymPy
@syms x, d
A = ((16x - 73d)*x + 52d^2)*x - 14d^3
eq = (A*x + 4d^4)*x - d^5 |> expand
eq |> println

    -d^5 + 4*d^4*x - 14*d^3*x^2 + 52*d^2*x^3 - 73*d*x^4 + 16*x^5

d = 4 のときの解を求める。5 個の解のうち，3 番目のものが適解 15.0181805396541 である。

res = solve(eq(d => 4));
res[3].evalf() |> println

    15.0181805396541