算額（その1409）

算額（その1409）

百五十 群馬県多野郡新町 稲荷神社 文政3年(1820)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：球3個，四角台
#Julia, #SymPy, #算額, #和算

四角台の中に，上球，中球，下球が互いに接し，側面，上面，下面にも接している。
上球，下球の直径が与えられたとき，中球の直径はいかほどか。

以下の図は四角台の側面を除いて，横倒しにした状態である。

算額（その1405）の拡張なので，おなじように解いてもよい。

しかし，単に球の直径を求めるだけならば，算額（その1393）などに述べた問題の3次元版なので，もっと簡単になる。

2直線に挟まれた互いに外接する円の直径は等比数列をなす。
同じく，底辺が正方形の四角錐の 4 平面に挟まれた互いに外接する球の直径は等比数列をなす。
後者は，1 つの側面と平行な丙面に射影すれば，2直線に挟まれた円が見えるであろうことから同じ問題であることがわかる。

よって，上球，中球，下球 の直径を d1, d2, d3，公比を p とすれば，d2 = d1*p，d3 = d2*p = d1*p^2 である。
d1 * d3 = d1^2 * p^2 = d2^2 なので，d2 = sqrt(d1 * d3) である。

中球の直径は，上球の直径と下球の直径を掛け合わせ，平方根をとれば求まる。

 

 

算額（その1408）

算額（その1408）

百四十九 群馬県多野郡新町 稲荷神社 文政3年(1820)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，半円，直角三角形
#Julia, #SymPy, #算額, #和算

半円の中に三角形（直角三角形）を作り，甲円，乙円，丙円，丁円を容れる。大斜（半円の直径），小斜（三角形の一番短い辺）が与えられたとき，丁円の直径はいかほどか。

半円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, r3)
丁円の半径と中心座標を r4, (x4, r4)
円周上にある直角三角形の頂点の座標を (x0, sqrt(R^2 - x0))
とおき，以下の連立方程式を解く。

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

@syms R, x0, r1, x1, r2, x2, r3, x3, r4, x4,
      大斜, 小斜, 中斜
R = 大斜/2
eq1 = sqrt((R - x0)^2 + (R^2 - x0^2)) - 小斜
eq2 = r2/(x2 + R) - r1/(x1 + R)
eq3 = r3/(R - x3) - r1/(R - x1)
eq4 = x3 - x1 - 2sqrt(r1*r3)
eq5 = x1 - x4 - 2sqrt(r1*r4)
eq6 = x4 - x2 - 2sqrt(r2*r4)
eq7 = x1 - x2 - 2sqrt(r1*r2)
eq8 = dist2(-R, 0, x0, sqrt(R^2 - x0^2), x1, r1, r1)
eq9 = 小斜 + sqrt(大斜^2 - 小斜^2) - 大斜 - 2r1;
# eqxx = (大斜 + sqrt(大斜^2 - 小斜^2) + 小斜)*r1 - 小斜*中斜
# eqxx = dist2(R, 0, x0, sqrt(R^2 - x0^2), x1, r1, r1);
# eqxx = dist2(-R, 0, x0, sqrt(R^2 - x0^2), x2, r2, r2)
# eqxx = dist2(R, 0, x0, sqrt(R^2 - x0^2), x3, r3, r3);
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq10], (x0, r1, x1, r2, x2, r3, x3, r4, x4))

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)
    (x0, r1, x1, r2, x2, r3, x3, r4, x4) = u
    return [
        -小斜 + sqrt(-x0^2 + 大斜^2/4 + (-x0 + 大斜/2)^2),  # eq1
        -r1/(x1 + 大斜/2) + r2/(x2 + 大斜/2),  # eq2
        -r1/(-x1 + 大斜/2) + r3/(-x3 + 大斜/2),  # eq3
        -x1 + x3 - 2*sqrt(r1*r3),  # eq4
        x1 - x4 - 2*sqrt(r1*r4),  # eq5
        -x2 + x4 - 2*sqrt(r2*r4),  # eq6
        x1 - x2 - 2*sqrt(r1*r2),  # eq7
        (r1^2*x0 - r1*x1*sqrt(-4*x0^2 + 大斜^2) - x0*x1^2 + 大斜*(-4*r1^2 - 4*r1*sqrt(-4*x0^2 + 大斜^2) - 8*x0*x1 - 2*x0*大斜 + 4*x1^2 + 4*x1*大斜 + 大斜^2)/8)/大斜,  # eq8
        -2*r1 - 大斜 + 小斜 + sqrt(大斜^2 - 小斜^2),  # eq9
    ]
end;

大斜 = 10
R = 大斜/2
小斜 = 4
iniv = BigFloat[3.45, 1.58, 2.58, 1.05, 0.01, 0.46, 4.29, 0.32, 1.16]
res = nls(H, ini=iniv)

    ([3.4, 1.5825756949558398, 2.582575694955839, 1.0456116707250647, 0.009826491126086901, 0.46246227038447935, 4.293577213300662, 0.31816569849213683, 1.163390997458993], true)

たとえば，大斜，小斜，中斜が 10, 4, 9.16515 のとき，丁円の直径は 0.6363313969842737 である。

全てのパラメータは以下のとおりである。

   x0 = 3.4;  r1 = 1.58258;  x1 = 2.58258;  r2 = 1.04561;  x2 = 0.00982649;  r3 = 0.462462;  x3 = 4.29358;  r4 = 0.318166;  x4 =1.16339

function draw(大斜, 小斜, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = 大斜/2
    中斜 = sqrt(大斜^2 - 小斜^2)
    (x0, r1, x1, r2, x2, r3, x3, r4, x4) = res[1]
    @printf("大斜，小斜，中斜が %g, %g, %g のとき，丁円の直径は %g である。\n", 大斜, 小斜, 中斜, 2r4)
    @printf("x0 = %g;  r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  r4 = %g;  x4 =%g\n", x0, r1, x1, r2, x2, r3, x3, r4, x4)
    y0 = sqrt(R^2 - x0^2)
    plot([-R, x0, R], [0, y0, 0], color=:blue, lw=0.5)
    circle(0, 0, R, beginangle=0, endangle=180)
    circle(x1, r1, r1)
    circle(x2, r2, r2, :green)
    circle(x3, r3, r3, :magenta)
    circle(x4, r4, r4, :orange)
    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(x0, sqrt(R^2 - x0^2), "(x0,sqrt(R^2-x0^2))", :blue, :center, :bottom, delta=delta)
        point(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=-delta)
        point(x2, r2, "乙円:r2\n(x2,r2)", :green, :center, delta=-delta)
        point(x3, r3, "丙円:r3,(x3,r3)", :magenta, :right, delta=-7delta, deltax=13delta)
        point(x4, r4, "丁円:r4,(x4,r4)", :orange, :right, delta=-5delta, deltax=13delta)
        point(R, 0, "大斜/2", :red, :left, :bottom, delta=delta)
        point(0, 2.6r2, "中斜", :blue, :center, mark=false)
        point((R + x0)/2, r1, "小斜", :blue, :center, mark=false, deltax=7delta)
        plot!(xlims=(-R - 3delta, R + 15delta), ylims!(-10delta, R + 3delta))
    end
end;

draw(10, 4, true)

算額（その1407）

算額（その1407）

百四十三 群馬県榛名町榛名山 榛名神社 明治33年(1900)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：楕円，直角三角形
#Julia, #SymPy, #算額, #和算

直角三角形の中に斜線を 2 本描き，区画された領域に楕円を容れる。鈎が 3 寸，股が 4 寸のとき，楕円の長径，短径はいかほどか。

注：問には明記されておらず，図も正確ではないが，斜線と鈎，股の交点は鈎，股の中点である。

鈎，股，楕円の長径，短径をそのまま変数「鈎」，「股」，「長径」，「短径」とする。

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

三角形に内接する楕円の長径・短径に関する「算法助術」の公式97を援用する。
三角形の底辺, 斜辺1, 斜辺2, 高さを与えて，長径，短径の関係式を返す関数 equation を定義する。

@syms 底辺::positive, 斜辺1::positive, 斜辺2::positive, 高さ::positive, 長径::positive, 短径::positive
equation(底辺, 斜辺1, 斜辺2, 高さ) = -(斜辺1^2 - 斜辺2^2)^2*高さ*短径^2 + (斜辺1^2 - 斜辺2^2)^2*短径^3 + 底辺^4*高さ*(2§高さ - 短径)^2 - 底辺^4*短径*(2高さ - 短径)^2 - 底辺^2*高さ*長径^2*(2高さ - 短径)^2;

同じ楕円が 2 つの直角三角形に内接しているので，以下の連立方程式を解く。

@syms 鈎, 股
eq1 = equation(股, sqrt(股^2 + (鈎/2)^2), 鈎/2, 鈎/2)
eq2 = equation(股/2, sqrt((股/2)^2 + 鈎^2), 鈎, 鈎)
res = solve([eq1, eq2], (長径, 短径))[4];

# 長径
p = res[1] |> simplify |> sympy.sqrtdenest |> simplify
p |> println
p(鈎 => 3, 股 => 4).evalf() |> println

    股*(3 - sqrt(5))/2
    1.52786404500042

# 短径
q = res[2] |> simplify
q |> println
q(鈎 => 3, 股 => 4).evalf() |> println

    鈎*(3 - sqrt(5))/2
    1.14589803375032

長径は股の (3 - √5)/2 倍，短径は鈎の (3 - √5)/2 倍である。
鈎，股が 3 寸，4 寸のとき，長径は 4(3 - √5)/2 = 1.5278640450004204，短径は 3(3 - √5)/2 = 1.1458980337503153 である。

function draw(鈎, 股, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (a, b) = (股, 鈎) .* ((3 - √5)/ 2)  ./ 2
    @printf("鈎，股が %g，%g のとき，楕円の長径，短径は %g，%g である。\n", 鈎, 股, 2a, 2b)
    plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
    ellipse(a, b, a, b, color=:red)
    segment(0, 鈎/2, 股, 0, :green)
    segment(0, 鈎, 股/2, 0, :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(0, 鈎, "鈎 ", :blue, :right, :vcenter)
        point(0, 鈎/2, "鈎/2 ", :blue, :right, :vcenter)
        point(股, 0, "股", :blue, :center, delta=-delta)
        point(股/2, 0, "股/2", :blue, :center, delta=-delta)
        point(a, b, "楕円:a,b,(a,b)", :red, :center, delta=-delta)
        plot!(xlims=(-10delta, 股 + 5delta), ylims=(-7delta, 鈎 + 5delta))
    end
end;

draw(3, 4, true)

算額（その1406）

算額（その1406）

百四十三 群馬県榛名町榛名山 榛名神社 明治33年(1900)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円7個，外円，楕円，曲率円
#Julia, #SymPy, #算額, #和算

外円の中に，楕円と 2 個の交わる甲円を容れ，それらによって区切られた領域に乙円，丙円を 2 個ずつ容れる。乙円の直径が 1 寸で，楕円の短径がもっとも大きくなる（乙円が曲率円になる）ときの丙円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, r3)
乙円の半径と中心座標を r2, (R - r2, 0)
丙円の半径と中心座標を r3, (R - r3)
楕円の長半径，短半径と中心座標を R, R - 2r3, (0, 0)
とおき，以下の連立方程式を解く。

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

@syms R, r1, r2, r3
r1 = R - r3 
eq1 = (R - r2)^2 + r3^2 - (r1 + r2)^2  # 乙円は甲円と外接
eq2 = (R - 2r3)^2/R - r2  # 乙円は楕円の曲率円
res = solve([eq1, eq2], (r3, R))[3]  # 3 of 3

    (r2*(sqrt(2) + 2)/2, r2*(2*sqrt(2) + 3))

丙円の半径 r3 は，乙円の半径 r2 の (√2 + 2)/2 倍である。
乙円の直径が 1 寸で，楕円の短径がもっとも大きくなる（乙円が曲率円になる）とき，丙円の直径は (√2 + 2)/2 = 1.7071067811865475 寸である。

術では「天 = √8 + 3; (天 - √天)*乙円径/2」と複雑な計算式を提示しているが，「(√2/2 + 1)*乙円径」と簡約化できる。

全てのパラメータは以下のとおりである。

r2 = 0.5;  r3 = 0.853553;  R = 2.91421;  r1 = 2.06066

function draw(r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = r2*(2√2 + 3)
    r3 = r2*(√2 + 2)/2
    r1 = R - r3
    @printf("乙円の直径が %g のとき，丙円の直径は %g である。\n", 2r2, 2r3)
    @printf("r2 = %g;  r3 = %g;  R = %g;  r1 = %g\n", r2, r3, R, r1)
    plot()
    circle(0, 0, R)
    circle22(0, R - r3, r3, :blue)
    circle22(0, r3, r1, :green)
    circle(R - r2, 0, r2, :magenta)
    ellipse(0, 0, R, R - 2r3, color=:orange)
    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, r3, "甲円:r1,(0,r3)", :green, :center ,delta=-delta/2)
        point(R - r2, 0, "乙円:r2\n(R-r2,0)", :magenta, :center, :bottom, delta=delta/2)
        point(0, R - r3, "丙円:r3,(0,R-r3)", :blue, :center, delta=-delta/2)
        point(0, R, "R", :red, :center, :bottom, delta=delta/2)
        point(R, 0, " R", :red, :left, :vcenter)
        point(0, R - 2r3, "R-2r3", :orange, :center, :bottom, delta=delta/2)
    end
end;

draw(1/2, true)