裏 RjpWiki

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

算額(その78)

2022年12月26日 | Julia

算額(その78)

福島県船引町 蚕養国神社 明治24年(1891)2月
http://www.wasan.jp/fukusima/kogaikuni.html

外円の中に 7 個の円が入っている。乙円の径が二寸のとき甲円の径はいくつか。

図のように記号を定め方程式を解く。

外円の半径を 1 とする。甲円の半径 r1 = 1/3,乙円の半径と中心座標を r2,(x2, r1+r2) とする。

using SymPy
@syms r1::positive, r2::positive, x2::positive;
r1 = 1//3
eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x2^2 + (r1 + r2)^2 - (1 - r2)^2;
res = solve([eq1, eq2])

println(res)

   Dict{Any, Any}[Dict(r2 => 2/9, x2 => 2*sqrt(6)/9)]

r1 = 1/3, r2 = 2/9 ゆえ,r1 = 3/2 * r2。算額は,「置く一個五分 乗乙径 甲径を得る」という。和算では「一個五分」とは 1.5 = 3/2 を表すのだなあ。

r2 = 2 寸ならば,r1 = 3 寸。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
 θ = beginangle:0.1:endangle
 x = r.*cosd.(θ)
 y = r.*sind.(θ)
 plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

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, vertical, position, color))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
 plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, x2) = (1/3, 2/9, 2*sqrt(6)/9)
   plot()
   circle(0, 0, 1, :black)
   circle(0, 0, r1, :green)
   circle(0, 1-r1, r1, :green)
   circle(0, r1-1, r1, :green)
   circle(x2, r1+r2, r2, :blue)
   circle(x2, -r1-r2, r2, :blue)
   circle(-x2, r1+r2, r2, :blue)
   circle(-x2, -r1-r2, r2, :blue)
   circle(1-r1, 0, r1, beginangle=270, endangle=450)
   circle(r1-1, 0, r1, beginangle= 90, endangle=270)
   segment(r1-1, r1, 1-r1, r1, :red)
   segment(r1-1, -r1, 1-r1, -r1, :red)
   if more
       println("甲円の径は $(3/2 * 2)")
       point(0, 0, " O")
       point(0, 1-r1, " 1-r1")
       point(r1, 0, " r1")
       point(x2, r1+r2, "(x2,r1+r2)", :magenta, :top)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

   甲円の径は 3.0

 

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

算額(その77)

2022年12月26日 | Julia

算額(その77)

福島県船引町 蚕養国神社 明治24年(1891)2月
http://www.wasan.jp/fukusima/kogaikuni.html

正方形の中に 7 個の円が入っている。乙円の径が一寸二分のとき甲円の径はいくつか。

図のように記号を定め方程式を解く。

甲円の半径 r1 = 1/3,乙円の半径を r2 とする。

using SymPy
@syms r1::positive, r2::positive;
r1 = 1//3
eq1 = (1 - r2)^2 + r1^2 - (2r1 + r2)^2

$\left(1 - r_{2}\right)^{2} - \left(r_{2} + \frac{2}{3}\right)^{2} + \frac{1}{9}$

res = solve(eq1, r2)
println(res)

   Sym[1/5]

r1 = 1/3, r2 = 1/5 ゆえ,r1 = 5/3 * r2。算額の通り,「乙円の径を 5 倍して 3 で割る」。

r2 = 1.2 寸ならば,r1 = 2 寸。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
 θ = beginangle:0.1:endangle
 x = r.*cosd.(θ)
 y = r.*sind.(θ)
 plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

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, vertical, position, color))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (1/3, 1/5)
   plot([-1, 1, 1, -1, -1], [-1, -1, 1, 1, -1], color=:black, lw=0.5)
   circle(0, 0, r1, :green)
   circle(0, 1-r1, r1, :green)
   circle(0, r1-1, r1, :green)
   circle(0, r1, 2r1, :blue)
   circle(0, -r1, 2r1, :blue)
   circle(1 - r2, 0, r2, :magenta)
   circle(r2 - 1, 0, r2, :magenta)
   if more
       println("甲円の直径 = $(5/3 *1.2)")
       point(0, 0, " O")
       point(0, 1-r1, " 1-r1")
       point(r1, 0, " r1")
       point(1-r2, 0, "1-r2 ", :magenta, :right)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

    甲円の直径 = 2.0
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その76)

2022年12月26日 | Julia

算額(その76)

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

一般の三角形において,辺の長い順に大斜,中斜,小斜という。それぞれが 4 寸,2 寸 7 分,1 寸 8 分であるとき,内接する円の径を求めよ。

図のように記号を定め方程式を解けば,半径は求まる。単位は「分」で表しておく。

using SymPy
@syms r::positive, s::positive, S::positive;
(大斜, 中斜, 小斜) = (40, 27, 18)
s = (大斜 + 中斜 + 小斜)//2
eq1 = (大斜 + 中斜 + 小斜)*r//2 - S
eq2 = sqrt(Sym(s) * (s - 大斜) * (s - 中斜) * (s - 小斜)) - S  # ヘロンの公式
solve([eq1, eq2])

   Dict{Any, Any} with 2 entries:
     S => 35*sqrt(527)/4
     r => 7*sqrt(527)/34

7*sqrt(527)/34 * 2  # 直径

   9.452668468557999

もとの単位では 9分4厘5毛である。算額では,径は 1寸7厘5毛となっている。

どちらが正しいか,図を描いて確かめると,得られた解で正しいようである。算額における三角形とは随分異なる。

BF = BE = a, CD = CF = b, AE = AD = c, ∠CAB = α, ∠ABC = β として,プログラム中で必要な座標を独自に計算する。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
 θ = beginangle:0.1:endangle
 x = r.*cosd.(θ)
 y = r.*sind.(θ)
 plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

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, vertical, position, color))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
  plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (大斜, 中斜, 小斜) = (40, 27, 18)
   α = acos((小斜^2 + 大斜^2 - 中斜^2) / (2*小斜*大斜))
   β = acos((中斜^2 + 大斜^2 - 小斜^2) / (2*中斜*大斜))
   a = 中斜/2 + 大斜/2 - 小斜/2
   b = 中斜/2 - 大斜/2 + 小斜/2
   c = -中斜/2 + 大斜/2 + 小斜/2
   s = (大斜 + 中斜 + 小斜)/2
   S = sqrt(s * (s - 大斜) * (s - 中斜) * (s - 小斜))  # ヘロンの公式
   r = 2S/(大斜 + 中斜 + 小斜)
   plot([0, 大斜, 小斜*cos(α), 0], [0, 0, 小斜*sin(α), 0],
       lwd=0.25, linestyle=more ? :dot : :solid,
       xlims=(-3, 大斜+3), ylims=(-2, 小斜*sin(α)+1))
   circle(c, r, r)
   if more
       println("直径 = $(2r)")
       point(0, 0, "  α", :red, :bottom, :left, mark=false)
       point(大斜, 0, "β    ", :red, :bottom, :right, mark=false)
       segment(c, 0, c, r)
       point(c, r, " O", :black)
       point(c*cos(α), c*sin(α), "D   ", :red, :center, :right)
       segment(c*cos(α), c*sin(α), c, r)
       point(c, 0, "   E", :red, :top, :center)
       point(大斜 - a*cos(β), a*sin(β), "    F", :red, :center)
       segment(大斜 - a*cos(β), a*sin(β), c, r)
       point(0, 0, "A ", :black, :right)
       point(大斜, 0, " B", :black)
       point(小斜*cos(α), 小斜*sin(α), "  C", :black, :bottom)
       segment(0, 0, c, 0, :green)
       point(c/2, 0, "c", :green, mark=false)
       segment(大斜, 0, 大斜 - a*cos(β), a*sin(β), :magenta)
       point(大斜 - a*cos(β)/2, a*sin(β)/2, "a ", :magenta, :right, mark=false)
       segment(小斜*cos(α), 小斜*sin(α), c*cos(α), c*sin(α), :blue)
       point((小斜 + c)*cos(α)/2, (小斜+c)*sin(α)/2, "b", :blue, :bottom, :right, mark=false)
   end
end;

   直径 = 9.452668468557997

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

算額(その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でシェアする

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

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