算額（その1436）

算額（その1436）

福島市立子山 篠葉沢稲荷 明治10年(1877)
街角の数学 ～落書き帳「○△□」～ 68．ウルトラマン・ピーチ
http://streetwasan.web.fc2.com/math15.8.19.html
キーワード：円3個，半円2個，長方形
#Julia, #SymPy, #算額, #和算

外円の中に 2 本の斜線（等斜）を引き，区画された領域に 4 個の等円を容れる。外円の直径が 13 寸のとき，等円の直径はいかほどか。

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

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms R::positive, r::positive, x1::positive, y1::positive, y2::positive, x0::positive
t = sqrt(R^2 - x0^2)
eq1 = x1^2 + y1^2 - (R - r)^2
eq2 = y1* (t + R) - x0*x1  # y1/x1 - x0/(t + R)
eq3 = dist2(0, R, x0, -t, x1, y1, r) |> numerator
eq4 = dist2(0, R, x0, -t, r, y2, r) |> numerator
eq5 = sqrt(x0^2 + (R - t)^2) + sqrt(x0^2 + (R + t)^2) - 2R - 2r;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r, x1, y1, y2, 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=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 Float64.(v), r.f_converged
end;

function H(u)
    (r, x1, y1, y2, x0) = u
    t = sqrt(R^2 - x0^2)
    return [
        x1^2 + y1^2 - (R - r)^2,  # eq1
        -x0*x1 + y1*(R + sqrt(R^2 - x0^2)),  # eq2
        -4*R^3*r^2 + R^3*x0^2 - 4*R^3*x0*x1 + 4*R^3*x1^2 - 4*R^2*r^2*sqrt(R^2 - x0^2) - 2*R^2*x0^2*y1 + R^2*x0^2*sqrt(R^2 - x0^2) + 4*R^2*x0*x1*y1 - 4*R^2*x0*x1*sqrt(R^2 - x0^2) + 4*R^2*x1^2*sqrt(R^2 - x0^2) + 2*R*r^2*x0^2 + 2*R*x0^3*x1 - 3*R*x0^2*x1^2 + R*x0^2*y1^2 - 2*R*x0^2*y1*sqrt(R^2 - x0^2) + 4*R*x0*x1*y1*sqrt(R^2 - x0^2) - 2*x0^3*x1*y1 - x0^2*x1^2*sqrt(R^2 - x0^2) + x0^2*y1^2*sqrt(R^2 - x0^2),  # eq3
        x0*(-4*R^3*r + R^3*x0 + 4*R^2*r*y2 - 4*R^2*r*sqrt(R^2 - x0^2) - 2*R^2*x0*y2 + R^2*x0*sqrt(R^2 - x0^2) - R*r^2*x0 + 2*R*r*x0^2 + 4*R*r*y2*sqrt(R^2 - x0^2) + R*x0*y2^2 - 2*R*x0*y2*sqrt(R^2 - x0^2) - r^2*x0*sqrt(R^2 - x0^2) - 2*r*x0^2*y2 + x0*y2^2*sqrt(R^2 - x0^2)),  # eq4
        -2*R - 2*r + sqrt(x0^2 + (R - sqrt(R^2 - x0^2))^2) + sqrt(x0^2 + (R + sqrt(R^2 - x0^2))^2),  # eq5
    ]
end;

R = 13/2
iniv = BigFloat[1.9, 4.1, 1.7, -3.4, 4.6]
res = nls(H, ini=iniv)

    ([2.0, 4.153846153846154, 1.7307692307692308, -3.5, 4.615384615384615], true)

外円の直径が 13 寸のとき，等円の直径は 4 寸である。

function draw(R, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r, x1, y1, y2, x0) = [2.0, 4.153846153846154, 1.7307692307692308, -3.5, 4.615384615384615]
    t = sqrt(R^2 - x0^2)
    @printf("外円の直径が %g のとき，等円の直径は %g である。\n", 2R, 2r)
    @printf("R = %g;  r = %g;  x1 = %g;  y1 = %g;  y2 = %g;  x0 = %g\n", R, r, x1, y1, y2, x0)
    plot()
    circle(0, 0, R)
    circle2(x1, y1, r, :blue)
    circle2(r, y2, r, :blue)
    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)
        plot!([0, x0, 0, -x0, 0, 0], [R, -t, -R, -t, R, -R], color=:green, lw=0.7)
        point(x1, y1, "等円:r,(x1,y1)", :blue, :center, delta=-delta/2)
        point(r, y2, "等円:r,(r,y2)", :blue, :center, delta=-delta/2)
        point(x0, -t, " (x0,sqrt(R^2-x0^2)", :green, :left, :vcenter)
        point(0, R, "R", :red, :center, :bottom, delta=delta/2)
    end
end;

draw(13/2, true)