裏 RjpWiki

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

算額(その75)

2022年12月26日 | Julia

算額(その75)

大阪府茨木市 井於神社 弘化3年(1846)
http://www.wasan.jp/osaka/iyo.html

鉤が 2 寸 7 分,股が 3 寸 6 分,弦が 4 寸 5 分の直角三角形の内部に正三角形が内接している。正三角形の一辺の長さを求めよ。

寸法を 10 倍し(「分」で表し)図のように記号を定め,方程式を解く。

算額の問題文に明確な記載はないが,∠DBC = 45° が仮定されているのだろうか(これ以外の仮定をする場合の解を後半に示す)。

using SymPy
@syms a::positive, b::positive, x::positive
a = 36*27//(36+27)
eq1 = (a-b)^2 + a^2 - x^2
eq2 = x/sqrt(Sym(2)) - b
solve([eq1, eq2])

   1-element Vector{Dict{Any, Any}}:
    Dict(b => sqrt(2)*(-108*sqrt(2)/7 + 108*sqrt(6)/7)/2, x => -108*sqrt(2)/7 + 108*sqrt(6)/7)

-108*sqrt(Sym(2))/7 + 108*sqrt(Sym(6))/7 |> factor |> println

   -108*(-sqrt(6) + sqrt(2))/7

(√6 - √2)*108/7

   15.972832497755562

正三角形の一辺の長さは (√6 - √2)*108/7 = 15.972832497755562 となる。もとの単位では 1寸5分9厘7毛...であるが,算額では1寸5分5厘8毛になっている。

using Plots

function point(x, y, string="", color=:green, position=:left, vertical=:top; fontsize=10, mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, fontsize, position, color, vertical))
end;

function draw(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="")
    a = 36*27/(36+27)
    b = sqrt(2)*(-108*sqrt(2)/7 + 108*sqrt(6)/7)/2
   x = 108*(sqrt(6) - sqrt(2))/7
   println("x = $x")
    plot([0, 36, 0, 0], [0, 0, 27, 0], xlims=(-2, 37), ylims=(-2, 28))
    plot!([b, a, 0, b], [0, a, b, 0])
   if more == true
       point(0, 27, "A ", :black, :right)
       point(0, 0, "B ", :black, :right)
       point(36, 0, " C", :black)
       point(a, a, "   D (a,a)")
       point(0, b, "b ", :green, :right)
       point(a, 0, " a")
       point(b, 0, " b")
       plot!([0, a], [0, a])
   end
end;

---

上述の解は,左下にできる三角形が二等辺三角形であると仮定したものである。この仮定の代わりに左下にできる三角形が外形の三角形と相似であると仮定した場合の解について考察してみる。

方程式は以下のようになる。正三角形の一辺を x として,5本の方程式を立てて解く。

using SymPy
@syms Dy::positive, Ex::positive, Fx::positive, Fy::positive, x::positive;
eq1 = Ex^2 + Dy^2 - x^2
eq2 = (Fx - Ex)^2 + Fy^2 - x^2
eq3 = Fx^2 + (Fy - Dy)^2 - x^2
eq4 = Dy - x*27//45
eq5 = 27Fx + 36Fy - 27*36;  # 面積の制約
p = solve([eq1, eq2, eq3, eq4, eq5])

   1-element Vector{Dict{Any, Any}}:
    Dict(Ex => -6912/433 + 7200*sqrt(3)/433, Dy => -5184/433 + 5400*sqrt(3)/433, x => -8640/433 + 9000*sqrt(3)/433, Fx => 1008*sqrt(3)/433 + 4644/433, Fy => 8208/433 - 756*sqrt(3)/433)

p[1][x].evalf()

    16.0472454229097

function draw2(more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="")
    p = Dict(Dy => -5184/433 + 5400*sqrt(3)/433, Ex => -6912/433 + 7200*sqrt(3)/433, Fx => 1008*sqrt(3)/433 + 4644/433, x => -8640/433 + 9000*sqrt(3)/433, Fy => 8208/433 - 756*sqrt(3)/433)
    plot([0, 36, 0, 0], [0, 0, 27, 0], xlims=(-2, 37), ylims=(-2, 28))
   plot!([0, p[Ex], p[Fx], 0], [p[Dy], 0, p[Fy], p[Dy]])
   if more
       println("x = $(p[x])\nDy = $(p[Dy])\nEx = $(p[Ex])\n(Fx, Fy) = $((p[Fx], p[Fy]))")
       point(0, p[Dy], "Dy ", :green, :right)
       point(p[Ex], 0, " Ex")
       point(p[Fx], p[Fy], "  (Fx, Fy)")
   end
end
draw2(true)

   x = 16.047245422909686
   Dy = 9.628347253745813
   Ex = 12.837796338327747
   (Fx, Fy) = (14.757291487365883, 15.932031384475586)


コメント    この記事についてブログを書く
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする
« Julia 1.8.4 | トップ | 算額(その76) »
最新の画像もっと見る

コメントを投稿

Julia」カテゴリの最新記事