裏 RjpWiki

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

算額(その1606)

2025年02月10日 | Julia

算額(その1606)

長野県須坂市 とある廃寺 元治元年(1864)
深川英俊,トニー・ロスマン:聖なる数学:算額,森北出版株式会社,2010年4月22日.
キーワード:直角三角形,円弧,正方形,最大値
#Julia, #SymPy, #算額, #和算, #数学

直角三角形の中に円弧を描き,頂点が円弧と斜辺に接する正方形を容れる。股が一定のとき,正方形の一辺の長さは鈎の長さに応じて変化する。正方形の一辺の長さが最大になるのはどのようなときか。

鈎,股をそのまま変数名とする。
正方形の一辺の長さを a, 左下の頂点の座標を (x, 0)
とおき,以下の連立方程式を解く。

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

using SymPy
@syms 鈎::positive, 股::positive,
      x::positive, a::positive
eq1 = (股 -x - a)^2 + (鈎 - a)^2 - 鈎^2
eq2 = a/x - 鈎/股
res = solve([eq1, eq2], (a, x))[2]  # 2 of 2

    (股^2*鈎/(股^2 + 2*股*鈎 + 2*鈎^2), 股^3/(股^2 + 2*股*鈎 + 2*鈎^2))

# 正方形の一辺の長さが最大になるときの鈎
ans_鈎 = solve(diff(res[1], 鈎), 鈎)[1]
ans_鈎 |> println

    sqrt(2)*股/2

# 正方形の一辺の長さの最大値
res[1](鈎 => ans_鈎) |> simplify |> factor |> println

    股*(-1 + sqrt(2))/2

鈎 = 股*√2/2 のとき,正方形の一辺の長さは最大値の「股*(√2 - 1)/2」を取る。

# x 正方形の位置
res[2](鈎 => ans_鈎) |> simplify |> factor |> println

    -股*(-2 + sqrt(2))/2

股 = 10 のとき,鈎 = 7.07107;  正方形の一辺の長さ = 2.07107 である。

function draw(股, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    鈎 = √2股/2
    a = 股*(√2 - 1)/2
    x = -股*(-2 + sqrt(2))/2
    θ = atand(股/鈎)
    @printf("股 = %g;  鈎 = %g;  正方形の一辺の長さ = %g\n", 股, 鈎, a)
    plot([0, 股, 股, 0], [0, 0, 鈎, 0], color=:green, lw=0.5)
    plot!([x, x + a, x + a, x, x], [0, 0, a, a, 0], color=:blue, lw=0.5)
    circle(股, 鈎, 鈎, beginangle=270 - θ, endangle=270)
    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, 0, "x", :blue, :center, delta=-delta)
        point(x + a, 0, "x+a", :blue, :center, delta=-delta)
        point(x + a, a, " (x+a,a)", :blue, :left, :vcenter)
        point(x, a, "(x,a)", :blue, :right, :vcenter)
        point(0, 0, "0", :blue, :center, delta=-delta)
        point(股, 0, "股", :blue, :center, delta=-delta)
        point(股, 鈎, "(股,鈎)", :blue, :right, :bottom, delta=delta/2)
        ylims!(-5delta, 股 + delta)
    end
end;

draw(10, true)

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

算額(その1605)

2025年02月10日 | Julia

算額(その1605)

新潟県長岡市 悠久山 享和元年(1801)
深川英俊,トニー・ロスマン:聖なる数学:算額,森北出版株式会社,2010年4月22日.
キーワード:円5個,外円,弦,矢,斜線2本
#Julia, #SymPy, #算額, #和算, #数学

外円の中に,弦,矢,斜線,甲円,乙円を容れる。外円の直径が与えられ,「弦 - 矢」が最大になるとき,乙円の直径はいかほどか。

弦,矢をそのまま変数名とする。
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - 矢 + r1), (0, R - r1)
乙円の半径と中心座標を r2, (x2, R - 矢/2)
斜線と外円の交点座標を (x0, y0)
とおき,以下の連立方程式を解く。

まず初めに,「『弦 - 矢』が最大になるとき」という条件を解決する。
弦 = 2sqrt(R^2 - (R - 矢)^2) なので,
差 = 2sqrt(R^2 - (R - 矢)^2) - 矢 である。
矢について導関数を取り,導関数 = 0 をとき,矢を求める。

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

using SymPy
@syms 矢::positive, 弦::positive, R::positive,
      r1::positive, r2::positive, x2::positive
差 = 2sqrt(R^2 - (R - 矢)^2) - 矢;
ans_矢 = solve(diff(差, 矢), 矢)[1]
ans_矢 |> println

    R*(5 - sqrt(5))/5

矢 = R*(5 - √5)/5 である。弦 = 2sqrt(R^2 - (R - 矢)^2) を条件に加え,以下の連立方程式を解く。

@syms 矢::positive, 弦::positive, R::positive,
      r1::positive, r2::positive, x2::positive,
      x0::positive, y0::positive

矢 = R*(5 - √Sym(5))/5
弦 = 2sqrt(R^2 - (R - 矢)^2) |> simplify
eq1 = dist2(弦/2, R - 矢, 0, R - 矢/2, 0, R - r1, r1)

eq2 = x2^2 + (R - 矢/2)^2 - (R - r2)^2 |> expand |> simplify;
eq3 = r2/x2 - 矢/2/sqrt((弦/2)^2 + (矢/2)^2);

eq4 = (y0 - (R - 矢))/(弦/2 + x0) - 矢/弦
eq5 = x0^2 + y0^2 - R^2

res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, x2, x0, y0))[1];

# r2: 乙円の半径
res[2] |> println

    -R*(-sqrt(2) + sqrt(10))*(-2*sqrt(2590 - 710*sqrt(5)) + sqrt(11 - sqrt(5))*(-5*sqrt(2) + 5*sqrt(10)))/(160*sqrt(11 - sqrt(5)))

R の倍数であるが,式中に R があるので式が簡約化できない。
倍数だけを取り出して簡約化すると,非常に簡潔な式 R*(5 - √5)/20 になる。

@syms d
apart(res[2]/R, d) |> simplify |> factor |> println

    -(-5 + sqrt(5))/20

深川は R/(5 + √5) としているが,分母を有理化すれば同じ式である。

R/(5 + √Sym(5)) |> simplify |> factor |> println

    -R*(-5 + sqrt(5))/20

function draw(R, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    矢 = R*(5 - √5)/5
    弦 = 2sqrt(R^2 - (R - 矢)^2)
    (r1, r2, x2, x0, y0) = (-2*R - 2*sqrt(5)*R/5 - 2*sqrt(145)*R*sqrt(7 - 2*sqrt(5))/(5*(-7 + 2*sqrt(5))), -R*(-sqrt(2) + sqrt(10))*(-2*sqrt(2590 - 710*sqrt(5)) + sqrt(11 - sqrt(5))*(-5*sqrt(2) + 5*sqrt(10)))/(160*sqrt(11 - sqrt(5))), -R*(-2*sqrt(2590 - 710*sqrt(5)) + sqrt(11 - sqrt(5))*(-5*sqrt(2) + 5*sqrt(10)))/80, -11*sqrt(5)*R/145 - R/29 + sqrt(5)*R*sqrt(870 - 350*sqrt(5))/580 + sqrt(5)*R*sqrt(174 - 70*sqrt(5))/116 + 11*R*sqrt(870 - 350*sqrt(5))/580 + 11*R*sqrt(174 - 70*sqrt(5))/116, sqrt(5)*R*sqrt(870 - 350*sqrt(5))/580 + 11*R*sqrt(870 - 350*sqrt(5))/580 + 16*sqrt(5)*R/145 + 12*R/29)
    @printf("R = %g;  r1 = %g;  r2 = %g;  x2 = %g\n", R, r1, r2, x2)
    plot()
    circle(0, 0, R, :green)
    circle2(x2, R - 矢/2, r2)
    circle(0, R - r1, r1, :blue)
    circle(0, R - 矢 + r1, r1, :blue)
    segment(-弦/2, R - 矢, 弦/2, R - 矢, :magenta)
    segment(0, R, 0, R - 矢, :magenta)
    segment(-弦/2, R - 矢, x0, y0, :magenta)
    segment(弦/2, R - 矢, -x0, y0, :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)
        point(0, R, "R", :green, :center, :bottom, delta=delta/2)
        point(0, R - 矢/2, "R-矢/2", :magenta, :right, :vcenter, deltax=-4delta)
        point(0, R - 矢 + r1, "甲円:r1,(0,R-矢+r1)", :blue, :center, delta=-delta)
        point(0, R - r1, "甲円:r1,(0,R-r1)", :blue, :center, :bottom, delta=delta/2)
        point(x2, R - 矢/2, "乙円:r2,(x2,R-矢/2)", :red, :center, delta=-delta/2)
        point(弦/2, R - 矢, "(弦/2,R-矢)", :magenta, :right, delta=-delta/2)
    end
end;

draw(1/2, true)

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

算額(その1604)

2025年02月10日 | Julia

算額(その1604)

大阪住吉神社 弘化3年(1846)
福田理軒:順天堂算譜,明治6年(1873) 
深川英俊,トニー・ロスマン:聖なる数学:算額,森北出版株式会社,2010年4月22日.

福島県郡山市中田町上石 国見山稲荷 明治11年(1878)
五輪教一:黄金比の眠るほこら,日本評論社,2015年7月10日

キーワード:菱形2個,面積
#Julia, #SymPy, #算額, #和算, #数学

図のように,直線 l 上に一辺の長さが a の正方形 ABCD が固定されている。同じ大きさの動く正方形 EFGH が直線(上に載っていて,辺 CD に接してながら回転していて,l 上の接点をE,CD 上の接点を F とする。頂点 H を通り l に垂直な動く直線と l との交点を P,2点 B,G を通る直線との交点を T とする。PT の最大値を a で表せ。

TP = TH + a
TH/BH = GJ/BJ  -->  TH = BH*GJ/BJ
BH = a*sin(θ) + a*cos(θ) + a
GJ = FI + DF - a = a*sin(θ) + a*cos(θ) - a
BJ = GI + a = a*cos(θ) + a
TP = (a*sin(θ) + a*cos(θ) + a)*(a*sin(θ) + a*cos(θ) - a)/(a*cos(θ) + a) + a

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

using SymPy
@syms a::positive, θ::positive
TP = (a*sin(θ) + a*cos(θ) + a)*(a*sin(θ) + a*cos(θ) - a)/(a*cos(θ) + a) + a;

TP は θ の関数で,θ = 1 のあたりで最大値を取ることがわかる。

pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(TP(a => 1), xlims=(0, pi/2), xlabel="θ", ylabel="TP")

# 導関数を求める
g = diff(TP, θ);

# 導関数  = 0 となるときの θ を求める
ans_θ = solve(g, θ)[1]
ans_θ |> println

    2*atan(sqrt(-2 + sqrt(5)))

# 最大値を求める
TP(θ => ans_θ) |> simplify |> factor

SymPy ではこれ以上簡約化されないが,手計算により以下のようになる。

 

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

算額(その1603)

2025年02月10日 | Julia

算額(その1603)

松村吉徳:算法闕疑抄 寛文元年(1661)
深川英俊:例題で知る日本の数学と算額,p.29,森北出版株式会社,1988年2月20日.
キーワード:八角切灯籠,体積
#Julia, #SymPy, #算額, #和算, #数学

立方体の 8 個の角を切り取る。できる面が正三角形と正八角形になるような立体(八角切灯籠)を作った。各辺の長さが 12 のとき,この立体の体積を求めよ。

  

立体の辺を a とすれば,元の立方体の一辺の長さは a + 2a/√2 = a + √2a である。

using SymPy
@syms a
b = a/√Sym(2)
b |> println

    sqrt(2)*a/2

# 三角錐の底面積
S = b^2/2
S |> println

    a^2/4

# 三角錐の体積
V = S*b/3
V |> println

    sqrt(2)*a^3/24

V(a => 12).evalf() |> println

    101.823376490863

# 立方体の体積
V1 = (a + √Sym(2)*a)^3
V1 |> println

    (a + sqrt(2)*a)^3

V1(a => 12).evalf() |> println

    24314.8051789035

# 求める立体の体積
V2 = V1 - 8V
V2 |> println

    -sqrt(2)*a^3/3 + (a + sqrt(2)*a)^3

V2(a => 12).evalf() |> println

    23500.2181669766

立体の体積は 23500.2181669766 である。

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

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

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