算額（その1441）

算額（その1441）

二本松亀谷町坂 観世音 天保10年(1839)
『諸流算額起原一』（佐久間庸軒）

～落書き帳「○△□」～  81. 佐久間文庫
http://streetwasan.web.fc2.com/math15.9.5.html
キーワード：円3個，正三角形，接線
#Julia, #SymPy, #算額, #和算

正三角形内に内接円と甲円，乙円を容れ，3円の共通接線を引く。内接円の直径が 675 寸，甲円の直径が 81 寸のとき，乙円の直径はいかほどか。

最初は，後半に示すような全ての描画パラメータを求める方法で解いた。かなり複雑な7元連立方程式になったので，数値解しか求められなかった。
それでも，乙円は内接円の 1/27 倍の大きさであることは分かった。

その後で，「落書き帳」の著者が「三角関数の加法定理を用いてすんなり？解けます」というヒントを与えていたので，以下のように「簡単に」解けた。図を描くのは別として。

1.　乙円の直径のみを得る

内接円，甲円，乙円の半径を r0, r1, r2 とおく。
甲円の中心 A，内接円の中心 O，内接円と正三角形の斜辺との接点 B を結ぶ直角三角形において，∠AOB = θ1
乙円の中心 C，内接円の中心 O，内接円と正三角形の斜辺との接点 D を結ぶ直角三角形において，∠COD = θ2
とすれば θ1 + θ2 = 120°なので，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms r0, r1, r2, θ1, θ2
θ2 = PI/3 - θ1
eq1 = (r0 + r1)*cos(θ1) - (r0 - r1)
eq2 = (r0 + r2)*cos(θ2) - (r0 - r2)
res = solve([eq1, eq2], (r2, θ1))[2]  # 2 of 2

    (-r0*(sqrt(3)*r0*sqrt(-((r0 - r1)^2 - (r0 + r1)^2)/(r0 + r1)^2) - r0 + sqrt(3)*r1*sqrt(-((r0 - r1)^2 - (r0 + r1)^2)/(r0 + r1)^2) - 3*r1)/(sqrt(3)*r0*sqrt(-((r0 - r1)^2 - (r0 + r1)^2)/(r0 + r1)^2) + 3*r0 + sqrt(3)*r1*sqrt(-((r0 - r1)^2 - (r0 + r1)^2)/(r0 + r1)^2) + r1), acos((r0 - r1)/(r0 + r1)))

res[1](r0 => 675, r1 => 81).evalf() |> println
res[2](r0 => 675, r1 => 81).evalf() |> println

    25.0000000000000
    0.666946344503664

内接円，甲円の直径がそれぞれ 675 寸，81 寸のとき，乙円の直径は 25 寸である。

2.　全ての描画パラメータを得る

内接円の半径と中心座標を r0, (0, r0)
正三角形の一辺の長さを 2a; a = √3r0
甲円の半径と中心座標を r1, (x1, y1)
乙円の半径と中心座標を r2, (x2, y2)
内接円，甲円，乙円の共通接線と正三角形の斜辺との交点座標を (cx, cy), (dx, dy)
cy = √3(a + cx), dy = √3(a - dx)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms a::positive, r0::positive,
      r1::positive, x1::positive, y1::positive,
      r2::positive, x2::positive, y2::positive,
      cx::negative, cy::positive, dx::positive, dy::positive
a = r0*√Sym(3)
cy = √Sym(3)*(a + cx)
dy = √Sym(3)*(a - dx)

eq1 = x1^2 + (y1 -  r0)^2 - (r0 + r1)^2
eq2 = x2^2 + (y2 -  r0)^2 - (r0 + r2)^2
eq3 = dist2(cx, cy, dx, dy, x1, y1, r1)
eq4 = dist2(cx, cy, dx, dy, x2, y2, r2)
eq5 = dist2(cx, cy, dx, dy, 0, a/√Sym(3), r0)
eq6 = dist2(a, 0, 0, √Sym(3)a, x1, y1, r1)
eq7 = dist2(-a, 0, 0, √Sym(3)a, x2, y2, 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 Float64.(v), r.f_converged
end;

function H(u)
    (x1, y1, r2, x2, y2, cx, dx) = u
    return [
        x1^2 + (-r0 + y1)^2 - (r0 + r1)^2,  # eq1
        x2^2 + (-r0 + y2)^2 - (r0 + r2)^2,  # eq2
        (12*cx^2*dx^2 - 12*sqrt(3)*cx^2*dx*r0 - 12*cx^2*dx*x1 + 4*sqrt(3)*cx^2*dx*y1 + 9*cx^2*r0^2 + 6*sqrt(3)*cx^2*r0*x1 - 6*cx^2*r0*y1 - 4*cx^2*r1^2 + 3*cx^2*x1^2 - 2*sqrt(3)*cx^2*x1*y1 + cx^2*y1^2 + 12*sqrt(3)*cx*dx^2*r0 - 12*cx*dx^2*x1 - 4*sqrt(3)*cx*dx^2*y1 - 18*cx*dx*r0^2 + 12*cx*dx*r0*y1 - 4*cx*dx*r1^2 + 6*cx*dx*x1^2 - 2*cx*dx*y1^2 + 9*dx^2*r0^2 - 6*sqrt(3)*dx^2*r0*x1 - 6*dx^2*r0*y1 - 4*dx^2*r1^2 + 3*dx^2*x1^2 + 2*sqrt(3)*dx^2*x1*y1 + dx^2*y1^2)/(4*(cx^2 + cx*dx + dx^2)),  # eq3
        (12*cx^2*dx^2 - 12*sqrt(3)*cx^2*dx*r0 - 12*cx^2*dx*x2 + 4*sqrt(3)*cx^2*dx*y2 + 9*cx^2*r0^2 + 6*sqrt(3)*cx^2*r0*x2 - 6*cx^2*r0*y2 - 4*cx^2*r2^2 + 3*cx^2*x2^2 - 2*sqrt(3)*cx^2*x2*y2 + cx^2*y2^2 + 12*sqrt(3)*cx*dx^2*r0 - 12*cx*dx^2*x2 - 4*sqrt(3)*cx*dx^2*y2 - 18*cx*dx*r0^2 + 12*cx*dx*r0*y2 - 4*cx*dx*r2^2 + 6*cx*dx*x2^2 - 2*cx*dx*y2^2 + 9*dx^2*r0^2 - 6*sqrt(3)*dx^2*r0*x2 - 6*dx^2*r0*y2 - 4*dx^2*r2^2 + 3*dx^2*x2^2 + 2*sqrt(3)*dx^2*x2*y2 + dx^2*y2^2)/(4*(cx^2 + cx*dx + dx^2)),  # eq4
        cx*dx*(3*cx*dx - 2*sqrt(3)*cx*r0 + 2*sqrt(3)*dx*r0 - 3*r0^2)/(cx^2 + cx*dx + dx^2),  # eq5
        9*r0^2/4 - 3*sqrt(3)*r0*x1/2 - 3*r0*y1/2 - r1^2 + 3*x1^2/4 + sqrt(3)*x1*y1/2 + y1^2/4,  # eq6
        9*r0^2/4 + 3*sqrt(3)*r0*x2/2 - 3*r0*y2/2 - r2^2 + 3*x2^2/4 - sqrt(3)*x2*y2/2 + y2^2/4,  # eq7
    ]
end;

r0 = 675/2
r1 = 81/2
iniv = BigFloat[140, 689, 12, -217, 612, -225, 159]
res = nls(H, ini=iniv)

    ([140.29611541307898, 688.5, 12.499999999999389, -216.5063509461113, 612.499999999998, -224.83351829019105, 159.42740387849864], true)

内接円，甲円の直径が 675 寸，81 寸のとき，乙円の直径は 25 寸である。
2r2 = 25，2r0 = 675 なので，乙円の直径は内接円の直径の 1/27 倍である。
ちなみに，甲円の直径は内接円の直径の 3/25 倍になっている。

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r0 = 675/2
    r1 = 081/2
    (x1, y1, r2, x2, y2, cx, dx) = [140.29611541307898, 688.5, 12.499999999999389, -216.5063509461113, 612.499999999998, -224.83351829019105, 159.42740387849864]
    a =  √3r0
    cy = √3(a + cx)
    dy = √3(a - dx)
    plot([a, 0, -a, a], [0, √3a, 0, 0], color=:brown, lw=0.5)
    circle(0, a/√3, a/√3, :orange)
    circle(x1, y1, r1, :blue)
    circle(x2, y2, r2)
    segment(cx, cy, dx, dy, :green)
    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(a, 0, "a", :brown, :left, :bottom, delta=delta)
        point(0, r0, "O", :red, :center, delta=-delta)
        point(x1, y1, "A", :blue, :center, :vcenter, deltax=5delta)
        point(x2, y2, "C", :red, :right, :vcenter, deltax=-2delta)
        (x01, y01) = foot(a, 0, 0, √3a, 0, r0)
        point(x01, y01, " B", :orange, :left, :vcenter)
        (x02, y02) = foot(-a, 0, 0, √3a, 0, r0)
        point(x02, y02, "D ", :orange, :right, :vcenter)
    end
end;

draw(true)