裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

算額(その1467)

2024年12月13日 | Julia

算額(その1467)

福島県田村郡堀越村(現田村市船引町堀越明神前) 明石明神(明石神社) 奉納年不明
街角の数学 Street Wasan ~落書き帳「○△□」~ 339.村瀬義益の算額(その4)

http://streetwasan.web.fc2.com/math17.11.18.html

福島県木幡村(現二本松市木幡治家) 弁財神社 嘉永3年?(1850)
街角の数学 Street Wasan ~落書き帳「○△□」~ 340.村瀬義益の算額(その5)
http://streetwasan.web.fc2.com/math17.11.20.html

キーワード:円,直角三角形
#Julia, #SymPy, #算額, #和算

1. 等円の直径を求める場合

直角三角形の中に弦に接する 5 個の円を容れる。両端の円はそれぞれ鈎,股にも接している。鈎,股が 30 寸,40 寸のとき,等円の直径はいかほどか。

鈎,股,弦はそのまま変数名とする。
円の個数を n とする。
等円の半径と中心座標を r, (r, y), (x, r) とおき,以下の連立方程式を解く。

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

using SymPy
@syms r, x, y, 鈎, 股, 弦, n
# 弦はまえもって計算しておく
eq1 = (x - r)^2 + (y - r)^2 - (2(n - 1)*r)^2
eq2 = (y - r)/(x - r) -  鈎/股
eq3 = 鈎*r +  股*y + 弦*r - 鈎*股;
res = solve([eq1, eq2, eq3], (r, x, y))[1];

r, x, y ともに,分母は共通で (-4*n^2*股^2*鈎^2 + 8*n*股^2*鈎^2 + 弦^2*股^2 + 弦^2*鈎^2 + 2*弦*股^3 + 2*弦*股^2*鈎 + 2*弦*股*鈎^2 + 2*弦*鈎^3 + 股^4 + 2*股^3*鈎 - 2*股^2*鈎^2 + 2*股*鈎^3 + 鈎^4) である。

# r
ans_r = res[1](sqrt(股^2 + 鈎^2) => 弦)
ans_r |> println

    股*鈎*(-2*n*弦*股*鈎 + 弦*股^2 + 2*弦*股*鈎 + 弦*鈎^2 + 股^3 + 股^2*鈎 + 股*鈎^2 + 鈎^3)/(-4*n^2*股^2*鈎^2 + 8*n*股^2*鈎^2 + 弦^2*股^2 + 弦^2*鈎^2 + 2*弦*股^3 + 2*弦*股^2*鈎 + 2*弦*股*鈎^2 + 2*弦*鈎^3 + 股^4 + 2*股^3*鈎 - 2*股^2*鈎^2 + 2*股*鈎^3 + 鈎^4)

明石神社の場合,等円の個数が 5 股,鈎,股が 30 寸,40 寸のとき,弦は 50 寸で,等円の半径は 50/13 = 3.84615384615385 寸である。

ans_r(n => 5, 鈎 => 30, 股 => 40, 弦 => 50) |> println
ans_r(n => 5, 鈎 => 30, 股 => 40, 弦 => 50).evalf() |> println

    50/13
    3.84615384615385

弁財神社の場合,等円の個数が 4 股,鈎,股が 5 寸,12 寸のとき,弦は 13 寸で,等円の半径は 26/25 = 1.04 寸である。

ans_r(n => 4, 鈎 => 5, 股 => 12, 弦 => 13) |> println
ans_r(n => 4, 鈎 => 5, 股 => 12, 弦 => 13).evalf() |> println

    26/25
    1.04000000000000

2. 等円の個数を求める場合

弁財神社の算額は,等円の直径を与えて個数を聞く問題になっている。



上で得られた等円の半径 r を表す式を n について解けばよい。

eq = r - 股*鈎*(-2*n*弦*股*鈎 + 弦*股^2 + 2*弦*股*鈎 + 弦*鈎^2 + 股^3 + 股^2*鈎 + 股*鈎^2 + 鈎^3)/(-4*n^2*股^2*鈎^2 + 8*n*股^2*鈎^2 + 弦^2*股^2 + 弦^2*鈎^2 + 2*弦*股^3 + 2*弦*股^2*鈎 + 2*弦*股*鈎^2 + 2*弦*鈎^3 + 股^4 + 2*股^3*鈎 - 2*股^2*鈎^2 + 2*股*鈎^3 + 鈎^4);
res2 = solve(eq, n);

解は 2 個得られるが,鈎股弦,等円の半径を代入したとき整数になるものが適解である。

res_n_n = res2[1]
res_n_p = res2[2];

res_n_n(鈎 => 30, 股 => 40, 弦 => 50, r => 50/13) |> println
res_n_p(鈎 => 30, 股 => 40, 弦 => 50, r => 50/13) |> println

    3.50000000000000
    5.00000000000000

半径として概数を与えたときは両方とも整数にならないことがある。

res_n_n(鈎 => 30, 股 => 40, 弦 => 50, r => 3.846) |> println
res_n_p(鈎 => 30, 股 => 40, 弦 => 50, r => 3.846) |> println

    3.50000000000000
    5.00026001040042

res_n_n(鈎 => 5, 股 => 12, 弦 => 13, r => 26/25) |> println
res_n_p(鈎 => 5, 股 => 12, 弦 => 13, r => 26/25) |> println

    4.00000000000000
    4.25000000000000

解は 1 + 弦/4r ± sqrt(D1*r^2 + D2*r + D3)/(4r*鈎*股) の形をしているが,根号の中は全体としてかなり長い式である。

D1 = 4*r^2*弦^2*股^2 + 4*r^2*弦^2*鈎^2 + 8*r^2*弦*股^3 + 8*r^2*弦*股^2*鈎 + 8*r^2*弦*股*鈎^2 + 8*r^2*弦*鈎^3 + 4*r^2*股^4 + 8*r^2*股^3*鈎 + 8*r^2*股^2*鈎^2 + 8*r^2*股*鈎^3 + 4*r^2*鈎^4
D2 = - 4*r*弦*股^3*鈎 - 4*r*弦*股*鈎^3 - 4*r*股^4*鈎 - 4*r*股^3*鈎^2 - 4*r*股^2*鈎^3 - 4*r*股*鈎^4
D3 = 弦^2*股^2*鈎^2;

SymPy では簡約化できないが,res_n_n, res_n_p を以下のように簡約化できる。

res_n_n2 = 1 + 弦/4r - sqrt(4r*(股^2 + 鈎^2)*(弦 + 股 + 鈎)*(r*(弦 + 股 + 鈎) - 股*鈎) + (弦*股*鈎)^2)/(4r*股*鈎)
res_n_p2 = 1 + 弦/4r + sqrt(4r*(股^2 + 鈎^2)*(弦 + 股 + 鈎)*(r*(弦 + 股 + 鈎) - 股*鈎) + (弦*股*鈎)^2)/(4r*股*鈎);

明石神社の算額では適解は n = 5 である。

res_n_n2(r => 50/13, 鈎 => 30, 股 => 40, 弦 => 50) |> println
res_n_p2(r => 50/13, 鈎 => 30, 股 => 40, 弦 => 50) |> println

    3.50000000000000
    5.00000000000000

弁財神社の算額では適解は n = 4 である。

res_n_n2(r =>26/25, 鈎 => 5, 股 => 12, 弦 => 13) |> println
res_n_p2(r =>26/25, 鈎 => 5, 股 => 12, 弦 => 13) |> println

    4.00000000000001
    4.24999999999999

3. その他のパラメータ

以下は,図を描くために必要な等円の中心座標である。

# x
res[2](sqrt(股^2 + 鈎^2) => 弦) |> println

    股*鈎*(-4*n^2*股^2*鈎 + 2*n*弦^2*股 + 2*n*弦*股^2 + 8*n*股^2*鈎 - 2*弦^2*股 - 弦*股^2 + 弦*鈎^2 + 股^3 - 3*股^2*鈎 + 股*鈎^2 + 鈎^3)/(-4*n^2*股^2*鈎^2 + 8*n*股^2*鈎^2 + 弦^2*股^2 + 弦^2*鈎^2 + 2*弦*股^3 + 2*弦*股^2*鈎 + 2*弦*股*鈎^2 + 2*弦*鈎^3 + 股^4 + 2*股^3*鈎 - 2*股^2*鈎^2 + 2*股*鈎^3 + 鈎^4)

# y
res[3](sqrt(股^2 + 鈎^2) => 弦) |> println

    股*鈎*(-4*n^2*股*鈎^2 + 8*n*股*鈎^2 + 弦*股^2 + 弦*鈎^2 + 2*弦*鈎*(n - 1)*(弦 + 鈎) + 股^3 + 股^2*鈎 - 3*股*鈎^2 + 鈎^3)/(-4*n^2*股^2*鈎^2 + 8*n*股^2*鈎^2 + 弦^2*股^2 + 弦^2*鈎^2 + 2*弦*股^3 + 2*弦*股^2*鈎 + 2*弦*股*鈎^2 + 2*弦*鈎^3 + 股^4 + 2*股^3*鈎 - 2*股^2*鈎^2 + 2*股*鈎^3 + 鈎^4)

function draw(n, 鈎, 股, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    弦 = sqrt(鈎^2 + 股^2)
    denominator = (-4*n^2*股^2*鈎^2 + 8*n*股^2*鈎^2 + 弦^2*股^2 + 弦^2*鈎^2 + 2*弦*股^3 + 2*弦*股^2*鈎 + 2*弦*股*鈎^2 + 2*弦*鈎^3 + 股^4 + 2*股^3*鈎 - 2*股^2*鈎^2 + 2*股*鈎^3 + 鈎^4)
    r = 股*鈎*(-2*n*弦*股*鈎 + 弦*股^2 + 2*弦*股*鈎 + 弦*鈎^2 + 股^3 + 股^2*鈎 + 股*鈎^2 + 鈎^3)/denominator
    x = 股*鈎*(-4*n^2*股^2*鈎 + 2*n*弦^2*股 + 2*n*弦*股^2 + 8*n*股^2*鈎 - 2*弦^2*股 - 弦*股^2 + 弦*鈎^2 + 股^3 - 3*股^2*鈎 + 股*鈎^2 + 鈎^3)/denominator
    y = 股*鈎*(-4*n^2*股*鈎^2 + 8*n*股*鈎^2 + 弦*股^2 + 弦*鈎^2 + 2*弦*鈎*(n - 1)*(弦 + 鈎) + 股^3 + 股^2*鈎 - 3*股*鈎^2 + 鈎^3)/denominator
    @printf("鈎 = %g;  股 = %g;  弦 = %g;  等円の直径 = %.15g\n", 鈎, 股, 弦, 2r)
    p1 = plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw = 0.5)
    for i in 0:n - 1
        circle(r + i* (x - r)/(n - 1), y - i*(y - r)/(n - 1), r)
    end
    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, r, "等円:r1,(x,r)", :red, :right, delta=-delta/2, deltax=-r)
        point(r, y, "等円:r1,(r,y)", :red, :left, :bottom, delta=4delta, deltax=r)
        point(股, 0, " 股", :blue, :left, :bottom, delta=delta/2)
        point(0, 鈎, " 鈎", :blue, :left, :bottom, delta=delta/2)
    end
    return p1
end;

draw(5, 30, 40, true)


コメント    この記事についてブログを書く
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする
« 本場 讃岐うどん たも屋 三条店 | トップ |   
最新の画像もっと見る

コメントを投稿

Julia」カテゴリの最新記事