算額（その1607）

算額（その1607）

宮城県石巻市尾崎宮下 久須師神社 明治20年(1887)
徳竹亜紀子，谷垣美保：2021年度の算額調査，仙台高等専門学校名取キャンパス研究紀要，第 58 号（2022），p.7-28.
https://www.sendai-nct.ac.jp/natori-library/wp-content/uploads/2022/03/kiyo2022-2.pdf
キーワード：円2個，円弧2個，斜線
#Julia, #SymPy, #算額, #和算, #数学

交差する円弧の中に等円を 2 個と共通接線を入れる。短径，長径，等円の直径が与えられたとき，共通接線と円弧の交点間の距離（斜の長さ）を求める術を述べよ。

等円の半径を r
円弧を構成する円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x, y)
共通接線と円弧の交点座標を (x0, y0)
短径を 2K1, 長径を 2K2
とおく。

AB= K1 = R - y
BC = K2 = sqrt(R^2 - y^2) である。
以下の連立方程式を解き，(x0, y0, y, x, R) を求める。
これらは，K1, K2 を含む式で表される。

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

using SymPy
@syms R::positive, x::positive,
      y::positive, r::positive,
      x0::positive, y0::positive,
      K1::positive, K2::positive
eq1 = x^2 + y^2 - (R - r)^2
eq2 = r/x - (y0 - y)/sqrt(x0^2 + (y0 - y)^2)
eq3 = x0^2 + y0^2 - R^2
eq4 = (R - y) - K1
eq5 = sqrt(R^2 - y^2) - K2
res = solve([eq1, eq2, eq3, eq4, eq5], (x0, y0, y, x, R))[4]

    (K2^2*sqrt(-1/(K1*r - K2^2))*sqrt((K1^2*r - K1*K2^2 + K2^2*r)/(K1*r - K2^2))/sqrt(K1), (-K1^3*r + K1^2*K2^2 - K1*K2^2*r - K2^4)/(2*K1*(K1*r - K2^2)), (-K1^2 + K2^2)/(2*K1), sqrt(-(K1 - r)*(K1*r - K2^2))/sqrt(K1), (K1^2 + K2^2)/(2*K1))

斜の長さは 2sqrt(x0^2 + (y0 - y)^2) である。

# x0^2
x2 = res[1]^2 |> simplify
x2 |> println

    K2^4*(-K1^2*r + K1*K2^2 - K2^2*r)/(K1*(K1*r - K2^2)^2)

# (y0 - y)^2
y2 = (res[2] - res[3])^2 |> simplify
y2 |> println

    K2^4*r^2/(K1^2*r^2 - 2*K1*K2^2*r + K2^4)

# 斜 = 2sqrt(x0^2 + (y0 - y)^2) 
斜 = x2 + y2 |> simplify |> x -> sqrt(x) |> simplify |> x -> 2x

斜(K1 => 8/2, K2 => 16/2, r => 5/2) |> println

    10.6666666666667

徳竹も方法（経路）は違うが，結果として同じ式を導いている。

# 短径，長径，等円径はそれぞれの 1/2 で指定する
@syms 短径, 長径, 等円径
徳竹 = 長径^2*sqrt((短径 - 等円径)/(短径*(長径^2 - 短径*等円径)))

2*徳竹(短径 => 8/2, 長径 => 16/2, 等円径 => 5/2)|> println

    10.6666666666667

function draw(K1, K2, r, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (x0, y0, y, x, R) = (K2^2*sqrt(-1/(K1*r - K2^2))*sqrt((K1^2*r - K1*K2^2 + K2^2*r)/(K1*r - K2^2))/sqrt(K1), (-K1^3*r + K1^2*K2^2 - K1*K2^2*r - K2^4)/(2*K1*(K1*r - K2^2)), (-K1^2 + K2^2)/(2*K1), sqrt(-(K1 - r)*(K1*r - K2^2))/sqrt(K1), (K1^2 + K2^2)/(2*K1))
    斜 = sqrt(x0^2 + (y0 - y)^2)
    @printf("短径 = %g，長径 = %g，等円径 = %g のとき，斜の長さ = %g\n", 2K1, 2K2, 2r, 2斜)
    plot()
    θ = atand(y/sqrt(R^2 - y^2))
    circle(0, 0, R, beginangle=θ, endangle=180-θ, n=1000)
    circle(0, 2y, R, beginangle=180+θ, endangle=360-θ, n=1000)
    circle(x, y, r, :green)
    circle(-x, y, r, :green)
    segment(-x0, y - (y0 - y), x0, y0, :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)
        point(x, y, "等円:r,(x,y)", :green, :center, delta=-delta)
        point(0, R, "R", :red, :center, :bottom, delta=delta)
        point(0, y, " y", :blue, :left, :vcenter)
        point(x0, y0, "(x0,y0)", :blue, :left, :bottom, delta=delta)
        dimension_line(0, y, sqrt(R^2 - y^2), y, "K2", :magenta, :center, :bottom, delta=2delta)
        dimension_line(0, R, 0, y, "K1 ", :orange, :right, :vcenter)
    end
end;

draw(8/2, 16/2, 5/2, true)png)

 

 

算額（その616）改訂版

算額（その616）改訂版

より鮮明な写真により，大円の中の小円は互いに外接しているのではないようなので，算額（その616）の改訂版を書く。

福島県白河市明神 境明神 万延元年(1860)
http://www.wasan.jp/fukusima/sakai1L2.html
拡大図
http://www.wasan.jp/fukusima/sakai1-2.png

団扇の中に小円 1 個と長方形が内接し，長方形の中に大円 2 個，小円 7 個が入っている。
小円の径を 1 としたとき，団扇の径はいかほどか。

団扇（外円）の半径と中心座標を R, (R0, 0); R0 < 0
大円の半径と中心座標を r2, (r2 - r, 0), (r - r2, 0)
小円の半径と中心座標を r，(x, r2 - r), (2r2 - 2r, 0), (0, 0)
長方形の右上の頂点座標は (2r2 - r, r2) である。

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

using SymPy
@syms r::positive, x::positive, r2::positive
eq1 = (x + r2 - r)^2 + (r2 - r)^2 - (r2 + r)^2
eq2 = (2r2 - 2r - x)^2 + (r2 - r)^2 - 4r^2
res = solve([eq1, eq2], (r2, x))[3]

    (-(-700*r^3*(8*sqrt(6)/25 + 22/25) - 58*r^3*(8*sqrt(6)/25 + 22/25)^2 + 40*r^3 + 175*r^3*(8*sqrt(6)/25 + 22/25)^3)/(192*r^2), r*(8*sqrt(6)/25 + 22/25))

# r2
ans_r2 = res[1] |> simplify
ans_r2 |> println

    3*r*(4*sqrt(6) + 11)/25

# x
ans_x = res[2] |> factor
ans_x |> println

    2*r*(4*sqrt(6) + 11)/25

外円の半径は (r2 + 2r - R0) で長方形の頂点を通ることから，中心座標（y 座標） R0 < 0 を求める。

@syms R0
eq3 = (2r2 - r)^2 + (r2 - R0)^2 - (r2 + 2r - R0)^2;
ans_R0 = solve(eq3, R0)[1];

# R0
ans_R0 |> println

    3*r/4 + 2*r2 - r2^2/r

# R0  r だけを含む式
ans_R0 = ans_R0(r2 => ans_r2) |> simplify
ans_R0 |> println

    3*r*(221 - 256*sqrt(6))/2500

# R
R = ans_r2 + 2r - ans_R0 |> simplify
R |> println

    r*(1968*sqrt(6) + 7637)/2500

R(r => 1/2).evalf() |> println
2*R(r => 1/2).evalf() |> println

    2.49151916275946
    4.98303832551892

小円の直径が 1 寸のとき，外円の直径は 4.98303832551892 寸である。

答は「答曰団扇径五寸」，術は「術曰置小径五之得団扇径合問」と術とは言えないものである。

この算額は問題の順番が答えになるように作られているので，4.98303832551892 の概数を答えにしたものであろうか。

function draw(r, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r2 = 3r*(4√6 + 11)/25
    x2 = 2r2/3
    println("r2 = $r2")
    R0 = 3r/4 + 2r2 - r2^2/r
    R = r2 + 2r - R0
    println("外円の直径は $(2R) である。")
    a = 2r2 - r
    b = r2
    plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:green, lw=0.5)
    circle(0, R0, R, :magenta)
    circle2(r2 - r, 0, r2, :blue)
    circle4(x2, r2 - r, r)
    circle2(2r2 - 2r,  0, r)
    circle(0, 0, r)
    circle(0, r2 + r, r)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        vline!([0], color=:gray80, lw=0.5)
        hline!([0], color=:gray80, lw=0.5)
        point(0, R0, "R0")
        point(a, b, "(2r2-r,r2)", :green, :left, :bottom, delta=delta/2)
        point(0, r2 + 2r, "(0,r2+2r)", :magenta, :center, :bottom, delta=delta/2)
        point(0, r2 + r, "小円:r\n(0,r2+r)", :red, :center, delta=-delta/2)
        point(x2, r2 - r, "小円:r\n(x2,r2-r)", :red, :center, delta=-delta/2)
        point(2r2 - 2r, 0, "小円:r\n(2r2-2r,0)", :red, :center, delta=-delta/2)
        point(r2 - r, 0, "大円:r2\n(r2-r,0)", :blue, :center, delta=2.5delta)
        point(r - r2, 0, "大円:r2\n(r-r2,0)", :blue, :center, delta=2.5delta)
    end
end;

draw(1/2, true)