裏 RjpWiki

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

算額(その339)

2023年07月19日 | Julia

算額(その339)

和算で遊ぼう!! 「三春まちなか寺子屋」2017レポート
https://miharukoma.com/wp-content/uploads/2018/01/%E4%B8%89%E6%98%A5%E3%81%BE%E3%81%A1%E3%81%AA%E3%81%8B%E5%AF%BA%E5%AD%90%E5%B1%8B2017%E3%83%AC%E3%83%9D%E3%83%BC%E3%83%88.pdf

3月 田村大元神社
外円の中に 2 個の等円があり,その交わった部分に共通弦を設ける。
外円の直径が 9 寸,弦の長さが 3 寸のとき,等円の直径を求めよ。


外円と等円の半径を r,r1 とし,以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms r::positive, r1::positive, 弦::positive;

eq = (r - r1)^2 + (弦//2)^2 - r1^2
r1 = solve(eq, r1)[1]
2r1 |> println

   r + 弦^2/(4*r)

等円の直径は,「外円の半径に弦の二乗を外円の直径の2倍で割ったものを足す」

外円の直径が 9寸,弦の長さが 3寸のとき,等円径は 5寸 である。

r = 9//2
弦 = 3
等円径 = r + 弦^2/(4*r)
@printf("等円径 = %g\n", 等円径)

   等円径 = 5

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, r1) = (9//2, 5//2)
   plot()
   circle(0, 0, r, :blue)
   circle(0, r - r1, r1, :red)
   circle(0, r1 - r, r1, :red)
   segment(-3//2, 0, 3//2, 0, :green, lw=2)
   if more
       point(3/2, 0, "  3/2", :magenta, :left, :bottom)
       point(r, 0, "r0=9/2 ", :blue, :right, :bottom)
       point(0, r - r1, " r-r1", :red)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その338)

2023年07月19日 | Julia

算額(その338)

和算で遊ぼう!! 「三春まちなか寺子屋」2017レポート
https://miharukoma.com/wp-content/uploads/2018/01/%E4%B8%89%E6%98%A5%E3%81%BE%E3%81%A1%E3%81%AA%E3%81%8B%E5%AF%BA%E5%AD%90%E5%B1%8B2017%E3%83%AC%E3%83%9D%E3%83%BC%E3%83%88.pdf
2月 諏訪神社

一辺の長さが 9寸の正方形の中に図のように直角三角形を入れる。中央に小さな正方形ができるが,この面積が直角三角形1個の面積と等しくなるようにするとき(すなわち外側の正方形を 5 等分する),直角三角形の鉤,股(それぞれ直角三角形の直角を挟む短い方の辺と長い方の辺)の長さを求めよ。

外側の正方形の一辺の長さを a とし,正方形の中心を原点とする。第1象限内にある頂点 A の座標を (x, y) とする。鉤は AB, 股は AC である。以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms x::positive, y::positive, AB::positive, AC::positive, a::positive

eq1 = (a//2 - x)^2 + (a//2 - y)^2 - AB^2
eq2 = (a//2 + x)^2 + (a//2 - y)^2 - AC^2
eq3 = AB*AC//2 - a^2//5
eq4 = AB^2 + AC^2 - a^2

solve([eq1, eq2, eq3, eq4], (x, y, AB, AC))

   2-element Vector{NTuple{4, Sym}}:
    (3*a/10, a/10, sqrt(5)*a/5, 2*sqrt(5)*a/5)
    (3*a/10, 9*a/10, sqrt(5)*a/5, 2*sqrt(5)*a/5)

最初の方が適解。
x = 3*a/10,y = a/10,鉤 = a/√5,股 = 2a/√5 である。

a = 9寸 のとき,鉤 = 4.024922359499621寸,股 = 8.049844718999243寸である。

a = 9
println("鉤 = $(a/√5)")
println("股 = $(2a/√5)")

   鉤 = 4.024922359499621
   股 = 8.049844718999243

直角三角形の面積は a/√5 * 2a/√5 / 2 = 2a^2/5 = 16.2, 正方形の面積の 1/5 は 81/5 = 16.2

a/√5 * 2a/√5 / 2, 2a^2/5 / 2, a^2/5

   (16.2, 16.2, 16.2)

図の全体を描くには,⊿ABC を90°回転して4個描けばよい。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 9
   (x, y) = (3*a/10, a/10)
   plot()
   # rect(-a/2, -a/2, a/2, a/2, :black)
   A = [0.0 -1.0; 1.0 0]
   xy = [[x, a/2, -a/2, x], [y, a/2, a/2, y]]
   for i = 1:4
       plot!(xy[1], xy[2], color=:blue, lw=1)
       xy = A * xy
   end
   if more
       point(x, y, " A:(x,y)\n(3*a/10,a/10)")
       point(a/2, a/2, "B:(a/2,a/2) ", :green, :right, :bottom)
       point(-a/2, a/2, " C:(-a/2,a/2)", :green, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その337)

2023年07月19日 | Julia

算額(その337)

和算で遊ぼう!! 「三春まちなか寺子屋」2017レポート
https://miharukoma.com/wp-content/uploads/2018/01/%E4%B8%89%E6%98%A5%E3%81%BE%E3%81%A1%E3%81%AA%E3%81%8B%E5%AF%BA%E5%AD%90%E5%B1%8B2017%E3%83%AC%E3%83%9D%E3%83%BC%E3%83%88.pdf
1月 龍穏院

長方形に三角形を内接させると,図のように甲,乙,丙の直角三角形ができる。長方形の面積が 16歩,甲と乙の面積がそれぞれ 4歩,2歩のとき,三角形の面積を求めよ。

中学レベルだと長方形の面積(直)から甲,乙,丙の3個の直角三角形の面積を引いた残りが求める三角形の面積とするであろう。

図のように,横 = a,縦 = b の長方形があり,c, d が横,縦を区分するとき,以下のようになる。

include("julia-source.txt");

using SymPy

@syms a, b, c, d, 甲, 乙, 丙, 直, 三角

直 = a*b
甲 = a*d/2
乙 = b*c/2
丙 = (a - c)*(b - d)/2;
三角 = 直 - 甲 - 乙 - 丙 |> expand;

これを SymPy で簡約化すると,a*b/2 - c*d/2 になる。見慣れない式になるが,図のように等積変形を行うと,このような式になる。

三角 = 直 - 甲 - 乙 - 丙
⊿ODF = □OAEB - ⊿OAD - ⊿OBF - ⊿DEF
     = □OAEB - ⊿OAD - ⊿OBX - ⊿DEX
     = ⊿BEX
     = a*(b - x)/2
     = a*(b - c*d/a)/2
     = (a*b - a*c*d/a)/2
     = a*b/2 - c*d/2
     = a*b/2 - (2乙/b)*(2甲/a)/2
     = 直/2 - 2甲*乙/直

三角 |> println

   a*b/2 - c*d/2

更に,甲,乙の面積の計算式から c, d を代入すると,以下のようになる。

直/2 - (2甲/a)*(2乙/b)/2 |> expand |> println

   a*b/2 - c*d/2

整理すると算額の術に書かれている式になる。
直/2 - 2甲*乙/直
直,甲,乙だけを用いるトリッキーな式であるが,実際には適用しにくい(甲,乙,丙のどれを使ってもよいのでないのは明らかだが)。最初に述べた中学レベルの解法のほうが間違いもない。

直/2 - 2甲*乙/直 |> expand |> println

   a*b/2 - c*d/2

s(直, 甲, 乙) = 直/2 - 2甲*乙/直;  # 関数定義

s(16, 4, 2)

   7.0

using Plots

function polygon(xys, color=:blue; lw=0.5, alpha=0.2, fill=false)
   (x, y) = ([], [])
   for xy in xys
       append!(x, xy[1])
       append!(y, xy[2])
   end
   append!(x, x[1])
   append!(y, y[1])
   if fill
       plot!(x, y, linecolor=color, lw=lw, seriestype=:shape, fillcolor=color, alpha=alpha) 
   else
       plot!(x, y, color=color, lw=lw) 
   end
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, c, d) = (7, 5, 2, 3)
   #(a, b, c, d) = (8, 2, 2, 1)
   x = d/a * c
   plot()
   rect(0, 0, a, b, :black, fill=false)
   polygon([(0, 0), (a, 0), (a, d)], :cadetblue4, alpha=0.8, fill=true)  # 甲
   polygon([(0, 0), (c, b), (0, b)], :chocolate1, alpha=0.8, fill=true)  # 乙
   polygon([(0, 0), (c, x), (0, b)], :chocolate1, alpha=0.8, fill=true)  # 乙
   polygon([(a, d), (a, b), (c, b)], :coral3, alpha=0.8, fill=true)  # 丙
   polygon([(a, d), (a, b), (c, x)], :coral3, alpha=0.8, fill=true)  # 丙
   polygon([(0, 0), (c, b), (a, d)], :blue4, alpha=0.3, fill=false)  # 三角
   segment(c, 0, c, b, :red, linestyle=:dash)
   if more
       annotate!([
               (0, 0, ("O ", 10, :right, :top)),
               (a, 0, (" A:a", 10, :left, :top)),
               (0, b, ("B:b ", 10, :right, :bottom)),
               (a, b, (" E:(a,b)", 10, :bottom)),
               (a, d, (" D:d", 10, :left)),
               (c, 0, (" C:c", 10, :left, :top)),
               (c, b, ("F", 10, :bottom)),
               (c, x, (" X:(c,x), x=c*d/a", 10, :white, :left, :top)),
               (4a/5, d/2, ("甲")),
               (0.5, b/2, ("乙")),
               (a - 1, (b+d)/2, ("丙"))
               ])
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       plot!(xlims=(-0.8, a + 0.5), ylims=(-0.5, b + 0.5))
   else
       plot!(showaxis=false)
   end
end;

 

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

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

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