裏 RjpWiki

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

算額(その98)

2023年11月10日 | Julia

算額(その98)

元の「算額その(98)」の改訂版である

東京都渋谷区渋谷 金王八幡神社 安政6年(1859)4月
https://gunmawasan.web.fc2.com/files/konnou-HP.pdf

金王八幡宮③「宝物館」(渋谷散歩③)
https://wheatbaku.exblog.jp/30049403/

一部が欠けている円内に円弧を置き,甲円 1 個,乙円 3 個を入れる。甲円,乙円の直径がそれぞれ 5 寸,4 寸のとき,円弧の弦の長さはいくらか。

外円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (0, r0 - r1)
乙円の半径と中心座標を r2, (x2, y2)
円弧の一部になる円の半径と中心座標を r3, (0, r0 - 2r1 - r3)
外円と円弧の交点座標を (x, y)

とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive, x2::positive,
      y2::positive, r3::positive;

y = r0 - 2r1 - 2r2
x = sqrt(r0^2 - y^2)
eq1 = x2^2 + y2^2 - (r0 - r2)^2
eq2 = x2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2
eq3 = x2^2 + (y2 -(r0 - 2r1 - r3))^2 - (r2 + r3)^2
eq4 = (r0^2 - y^2) + (r3 - 2r2)^2 - r3^2
res = solve([eq1, eq2, eq3, eq4], (r0, x2, y2, r3))

   1-element Vector{NTuple{4, Sym}}:
    (r1*(r1^2 + r1*r2 - r2^2)/((r1 - r2)*(r1 + r2)), r2*sqrt(4*r1^2 + 4*r1*r2 - 4*r2^2)/r1, -(r1^4 - 3*r1^2*r2^2 - r1*r2^3 + 2*r2^4)/(r1*(r1 - r2)*(r1 + r2)), r1*r2/(r1 - r2))

x の根号の中身 r0^2 - (r0 - 2*r1 - 2*r2)^2 の x0 に
r0 = r1*(r1^2 + r1*r2 - r2^2)/((r1 - r2)*(r1 + r2))
を代入する。

r0 = r1*(r1^2 + r1*r2 - r2^2)/((r1 - r2)*(r1 + r2))
r0^2 - (r0 - 2*r1 - 2*r2)^2 |> simplify |> println

   4*r2^3/(r1 - r2)

甲円と乙円の半径がそれぞれ r1 = 5/2, r2 = 4/2 のとき,弦の長さは 2sqrt(4*r2^3/(r1 - r2)) = 16 である。

draw(5/2, 4/2, true)

   r1 = 2.5;  r2 = 2;  r0 = 8.05556;  x2 = 4.30813;  y2 = 4.25556;  r3 = 10
   x = 8;  y = -0.944444;  弦の長さ = 2x = 16

2sqrt(4*r2^3/(r1 - r2)) |> simplify |> println

   4*r2^(3/2)*sqrt(1/(r1 - r2))

術で言う「弦は 2乙 / sqrt(甲/乙 - 1)」と同じである(甲=2r1, 乙=2r2)。

4r2/sqrt(2r1/2r2 - 1) |> simplify |> println

   4*r2^(3/2)/sqrt(r1 - r2)

function draw(r1, r2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="")#, fontfamily="IPAMincho")
   (r0, x2, y2, r3) = (r1*(r1^2 + r1*r2 - r2^2)/((r1 - r2)*(r1 + r2)), r2*sqrt(4*r1^2 + 4*r1*r2 - 4*r2^2)/r1, -(r1^4 - 3*r1^2*r2^2 - r1*r2^3 + 2*r2^4)/(r1*(r1 - r2)*(r1 + r2)), r1*r2/(r1 - r2))
   @printf("r1 = %g;  r2 = %g;  r0 = %g;  x2 = %g;  y2 = %g;  r3 = %g\n", r1, r2, r0, x2, y2, r3)
   y = r0 - 2r1 - 2r2
   x = sqrt(r0^2 - y^2)
   θ = round(atand(y/x))
   yy = y - (r0 - 2r1 - r3)
   θ2 = round(atand(yy/x))
   @printf("x = %g;  y = %g;  弦の長さ = 2x = %g\n", x, y, 2x)
    plot()
    circle(0, 0, r0, :red, beginangle=θ, endangle=180 - θ)
    circle(0, r0 - r1, r1, :blue)
    circle(0, r0 - 2r1 - r2, r2, :magenta)
    circle(x2, y2, r2, :magenta)
    circle(-x2, y2, r2, :magenta)
   circle(0, r0 - 2r1 - r3, r3, :brown, beginangle=θ2, endangle=180 - θ2)
   segment(-x, y, x, y, :green)
    if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x, y, "(x,y) ", :green, :right, :bottom, delta=delta)
       point(0, r0 - r1, " 甲円:r1\n (0,r0-r1)", :blue, :left, :bottom)
       point(x2, y2, "乙円:r2,(x2,y2)", :magenta, :center, :top, delta=-delta/2)
       point(0, r0 - 2r1 - r2, " 乙円:r2\n (0,r0-2r1-r2)", :black, :left, :vcenter)
       point(0, r0, " r0", :red, :left, :bottom, delta=delta/2)
       point(x/2, -1.3y, "円弧:r3\n(0,r0-2r1-r3)", :brown, mark=false)
    end
end;

draw(6/2, 4/2, true)

   r1 = 3;  r2 = 2;  r0 = 6.6;  x2 = 4.42217;  y2 = 1.26667;  r3 = 6
   x = 5.65685;  y = -3.4;  弦の長さ = 2x = 11.3137

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

算額(その493)

2023年11月10日 | Julia

算額(その493)

東京都渋谷区 金王神社(金王八幡宮) 元治元年(1864)
https://gunmawasan.web.fc2.com/files/konnou-HP.pdf

金王八幡宮③「宝物館」(渋谷散歩③)
https://wheatbaku.exblog.jp/30049403/

正方形の中に小さな正方形が 4 個,等円が 2 個,楕円が 1 個入っている。楕円は小さな正方形の頂点を通る。
等円の直径が 7392 寸のとき,小さな正方形の一辺の長さを求めよ。

小さな正方形の一辺の長さを x とする。
楕円の長径と短径を a, b とする。外側の正方形の一辺の長さは 2a である。
等円の半径を r とする。第一象限において,小さな正方形の頂点は (r, r) である。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r::positive, x::positive;

eq1 = x + r - a
eq2 = a - 2r - b
eq3 = r^2/a^2 + r^2/b^2 - 1
res = solve([eq1, eq2, eq3], (x, a, b))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (r*sqrt(2 + sqrt(5)), r*(1 + sqrt(2 + sqrt(5))), r*(-1 + sqrt(2 + sqrt(5))))

   x = 7607;  a = 11303;  b = 3911;  r = 3696

等円の半径が 7392/2 のとき,小さな正方形の一辺の長さは 7607.000116795436 であり,答,術 とともに一致する。

(7392/2)*sqrt(2 + sqrt(5))

   7607.000116795436

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 7392/2
   (x, a, b) = (r*sqrt(2 + sqrt(5)), r*(1 + sqrt(2 + sqrt(5))), r*(-1 + sqrt(2 + sqrt(5))))
   @printf("x = %g;  a = %g;  b = %g;  r = %g\n", x, a, b, r)
   plot([a, a, -a, -a,  a], [-a, a, a, -a, -a], color=:blue, lw=0.5)
   circle(0, a - r, r, :red)
   circle(0, r - a, r, :red)
   rect(r, r, a, a, :blue)
   rect(-r, r, -a, a, :blue)
   rect(r, -r, a, -a, :blue)
   rect(-r, -r, -a, -a, :blue)
   ellipse(0, 0, a, b, color=:green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, a - r, " a-r", :red, :left, :vcenter)
       point(r, r, " (r,r)", :blue, :left, :bottom, delta=delta)
       point(a, 0, " a=x+r ", :green, :right, :bottom, delta=delta)
       point(0, b, " b=a-2r", :green, :left, :top, delta=-delta)
       point(a - x, 0, " r=a-x", :black, :center, :bottom, delta=delta)
   else
       plot!(showaxis=false)
   end
end;

 

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

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

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