裏 RjpWiki

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

算額(その1609)

2025年02月12日 | Julia

算額(その1609)

福島県田村市船引町堀越明神前 明石神社(堀越明宮) 明治11年(1878)
http://www.wasan.jp/fukusima/horikosiakasi.html
キーワード:直角三角形,正五角形
#Julia, #SymPy, #算額, #和算, #数学

直角三角形の中に弦(斜辺)と股(底辺)を共有し,残りの頂点が鈎(高さ)上にある正五角形を容れる。
正五角形の一辺の長さが与えられたとき,股,鈎はいかほどか。

正五角形の一辺の長さを a, 正五角形を内接する円の半径を r とする。

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

using SymPy
@syms r::positive, 股::positive, a::positive
s36 = Sym(36)
s18 = s36/2
r = a/2sind(s36);

直角三角形の鈎は a*sqrt(sqrt(5)/8 + 5/8) + 2*a*sqrt(5/8 - sqrt(5)/8),
正五角形の一辺の長さが 1 のとき,鈎 = 2.12662702088010 である。

鈎 = 2a*sind(s36) + a*cosd(s18)
鈎 |> println
鈎(a => 1).evalf() |> println

    a*sqrt(sqrt(5)/8 + 5/8) + 2*a*sqrt(5/8 - sqrt(5)/8)
    2.12662702088010

直角三角形と内部にある直角三角形の相似関係から,股を求める。

eq = 鈎/股 - (r + r*cosd(s36))/(股 - a/2 - a*sind(s18))
ans_股 = solve(eq, 股)[1];

股は √5 が二重根号の中にある複雑な式になり SymPy では自動では簡約化(有理化)できない。
正五角形の一辺の長さが 1 のとき,直角三角形の股(底辺)は 2.92705098312484 である。

ans_股 |> println
ans_股(a => 1).evalf() |> println

    sqrt(5)*a*(sqrt(sqrt(5) + 5) + 2*sqrt(5 - sqrt(5)))/(-3*sqrt(5)*sqrt(5 - sqrt(5)) - sqrt(5)*sqrt(sqrt(5) + 5) + 5*sqrt(5 - sqrt(5)) + 5*sqrt(sqrt(5) + 5))
    2.92705098312484

SymPy の手を借りながら(詳細省略),手動で有理化する。

num = ans_股 |> numerator
den = ans_股 |> denominator
d1 = den.args[1]
d2 = den.args[2]
d3 = den.args[3]
d4 = den.args[4]
den2 = den *(d1 + d2 + d3 + d4) |> expand |> simplify
num2 = num *(d1 + d2 + d3 + d4) |> expand |> simplify
ans_股2 = num2/den2 |> sympy.sqrtdenest |> simplify |> factor;

簡潔な式になった。
股 = a*(5 + 3√5)/4 で,正五角形の一辺の長さが 1 のとき,直角三角形の股(底辺)は 2.92705098312484 である。

ans_股2 |> println
ans_股2(a => 1).evalf() |> println

    a*(5 + 3*sqrt(5))/4
    2.92705098312484

function draw(a, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r = a/2sind(36)
    股 = sqrt(5)*a*(sqrt(sqrt(5) + 5) + 2*sqrt(5 - sqrt(5)))/(-3*sqrt(5)*sqrt(5 - sqrt(5)) - sqrt(5)*sqrt(sqrt(5) + 5) + 5*sqrt(5 - sqrt(5)) + 5*sqrt(sqrt(5) + 5))
    鈎 = 2a*sind(36) + a*cosd(18)
    println("股 = $股, 鈎 = $鈎")
    plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
    plot!([0, a*sind(18), a*sind(18) + a, 2a*sind(18) + a, a*sind(18) + a/2, 0],
        [a*cosd(18), 0, 0, a*cosd(18), r + r*cosd(36), a*cosd(18)], color=:red, lw=0.5)
    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*sind(18), 0, "a*sind(18)", :red, :center, delta=-delta)
        point(a*sind(18) + a/2, 0, "a*sind(18)+a/2", :red, :left, delta=-delta, deltax=-3delta)
        point(股, 0, "股", :red, :center, delta=-delta)
        point(0, a*cosd(18), " a*cosd(18) ", :red, :left, :vcenter)
        point(0, 鈎, " 鈎=a*cosd(18) + 2rsind(36)", :red, :left, :vcenter)
        dimension_line(a*sind(18) + a/2, 0, a*sind(18) + a/2, r + r*cosd(36), " (a*sind(18)+a/2, r+r*cosd(36))", :blue, :center, :vcenter)
        segment(a*sind(18) + a/2, 0, a*sind(18) + a/2, r + r*cosd(36), :gray80)
        dimension_line(a*sind(18) + a/2, 0, 股, 0, " x=股-a*sind(18)-a/2", :blue, delta=2delta, dy=-5delta)
        ylims!(-8delta, 鈎 + 3delta)
    end
end;
draw(1, true)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1608)

2025年02月12日 | Julia

算額(その1608)

福島県田村市船引町堀越明神前 明石神社(堀越明宮) 明治11年(1878)
http://www.wasan.jp/fukusima/horikosiakasi.html

~落書き帳「○△□」~  696. 第10回「街角の問題」―堀越村明宮神社算額(その3)
http://streetwasan.web.fc2.com/math19.9.27.html
キーワード:円7個,直線上
#Julia, #SymPy, #算額, #和算, #数学

直線上に,東円 1 個,西円 2 個,南円 1 個,北円 3 個を載せる。東円の直径が 5 寸のとき北円の直径はいかほどか。

注:horikosiakasi.html の写真では,「北円径一寸四分□東円径幾何」とあり,「答曰東円径五寸」となっている。

東円の半径と中心座標を r1, (0, 2r3 - r1)
西円の半径と中心座標を r2, (r2, r2)
南円の半径と中心座標を r3, (0, r3)
北円の半径と中心座標を r4, (x4, y4), (0, r4)
とおき,以下の連立方程式を解く。

まず,右にある西円と中央の 北円が外接する条件 eq4 で,r2 を求めると,r2 = 4r4 であることがわかる。

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

using SymPy
@syms r1::positive, r2::positive, r3::positive, r4::positive, x4::positive, y4::positive
eq4 = r2^2 + (r2 - r4)^2 - (r2 + r4)^2;
ans_r2 = solve(eq4, r2)[1]
ans_r2 |> println

    4*r4

r2 = 2r4 である。また,y4 = 2r4 でもある。
これを条件として加え,残りの 3 方程式を解く。

horikosiakasi.html の「北円径一寸四分□東円径幾何」に基づき(r4 が与えられたとき),r1, r3, x4 を求める。

r2 = 4r4
y4 = 2r4
eq1 = r2^2 + (2r3 - r1 - r2)^2 - (r1 + r2)^2 |> expand
eq2 = x4^2 + (r3 - y4)^2 - (r3 + r4)^2 |> expand
eq3 = (x4 - r2)^2 + (y4 - r2)^2 - (r4 - r2)^2 |> expand
res = solve([eq1, eq2, eq3], (r1, r3, x4))[2]  # 2 of 2

    (r4*(7*sqrt(5) + 27)/12, 4*r4*(sqrt(5) + 3)/3, r4*(sqrt(5) + 4))

# r1 東円
res[1] |> println

    r4*(7*sqrt(5) + 27)/12

東円の半径 r1 は,北円の半径 r4 の (7√5 + 27)/12 = 3.5543729868748777 倍である。
北円の直径が 1 寸 4 分のとき,東円の直径は 1.4*(7√5 + 27)/12 = 4.976122181624828 寸である。

「術」は「北径を 0.28 で割れ」とあり,1.4/0.28 = 5 である。
実にきれいな結果になるが,0.28 は,r4 に掛ける倍数とすれば 1/0.28 = 3.571428571428571 なので,(7√5 + 27)/12 = 3.5543729868748777 の近似値ではなく,答えが 5 になるように選んだものであろう。

ということで,ちょっと汚いので,「落書き帳」ではこれを逆にしたのであろう。

「連立方程式を解く」というメソッドを取っているので,方程式はそのままで,何を求めるかを変えればよいだけである。

r1 が与えられるとして,r3, r4, x4 を求めればよい。

res2 = solve([eq1, eq2, eq3], (r3, r4, x4))[2]  # 2 of 2

    (24*sqrt(5)*r1/121 + 184*r1/121, -21*sqrt(5)*r1/121 + 81*r1/121, 6*r1*(73/242 - sqrt(5)/242))

# r3
res[1] |> simplify |> println

    r4*(7*sqrt(5) + 27)/12

# r4
res[2] |> simplify |> println

    4*r4*(sqrt(5) + 3)/3

# x4
res[3] |> simplify |> println

    r4*(sqrt(5) + 4)

北円の半径 r4 は,東円の半径 r1 の 3(27 - 7√5)/121 である。
東円の直径が 5 寸のとき,北円の直径は 1.4067178707646453 寸である。

これならば,「小数点以下 1 桁まで取って,北円の直径は 1 寸 4 分である」としても,さほどの違和感はない。(しかし,算額の原文を変更するのはあまり好ましいものとは思えない)

5*3(27 - 7√5)/121

    1.4067178707646453

function draw(r1, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r3 = r1*8(3√5 + 23)/121
    r4 = r1*3(27 - 7√5)/121
    x4 = r1*3(73 - √5)/121
    r2 = 4r4
    y4 = 2r4
    @printf("r1 = %g;  r2 = %g;  r3 = %g;  r4 = %g;  x4 = %g;  y4 = %g\n", r1, r2, r3, r4, x4, y4)
    plot()
    circle(0, r3, r3)
    circle2(r2, r2, r2, :blue)
    circle(0, 2r3 - r1, r1, :green)
    circle2(x4, y4, r4, :magenta)
    circle(0, r4, r4, :magenta)
    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)
        hline!([r4, 2r4, 3r4, 4r4], color=:orange, lw=0.5, style=:dot)
        point(0, 2r3 - r1, "東円:r1,(0,2r3-r1)", :green, :center, delta=-delta/2)
        point(r2, r2, "西円:r2,(r2,r2)", :blue, :center, delta=-delta/2)
        point(0, r3, "南円:r3,(0,r3)", :red, :center, delta=-delta/2)
        point(x4, y4, "北円:r4,(x4,y4)", :magenta, :center, delta=-delta/2)
        point(0, r4, "北円:r4,(0,r3)", :magenta, :center, delta=-delta/2)
    end
end;

draw(5/2, true)

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

高松市鹿角町 中西うどん

2025年02月12日 | さぬきうどん

高松市鹿角町 中西うどん
看板も骨太の太麺。女将さんがこっそり(?)出し始めた細麺もある。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村