裏 RjpWiki

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

算額(その792)

2024年03月19日 | Julia

算額(その792)

文化元年甲子九月 藤田貞資門人 北越村上藩 杉村勘助道定
藤田貞資(1807):続神壁算法

http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf

圓壔(えんとう;円柱)から弧環(トーラス?)で欠き取った鼓状の物体内に大球 2 個と,数個の小球が入っている。小球は大球と外側面および両隣の小球に接する(小球間の隙間はない)。円柱の底面(円)の直径が 35 寸 9 分のとき,円柱の高さはいかほどか。

円柱の高さ,底面の半径と中心座標を 4r2, r1, (0, 0, 2r2)
大球の半径と中心座標を r2, (0, 0, r2)
小球の半径と中心座標を r4, (x4, 0, 0)
とおき,以下の連立方程式を解くが,小球の個数を試行錯誤で決める必要がある(?)。
n = 10 のときに,条件を満たす解が得られた。

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

using SymPy
@syms r1::positive, r2::positive, r3::positive,
     x3::positive, r4::positive, x4::positive,
     h::positive, n::integer
h = 2r2
n = 10
eq1 = x3^2 + r2^2 - (r2 + r3)^2
eq2 = (x3 - r1)^2 + h^2 - r3^2
eq3 = x4^2 + r2^2 - (r2 + r4)^2
eq4 = x4 + r4 + r3 - x3
eq5 = x4*sind(Sym(180)//n) - r4
res = solve((eq1, eq2, eq3, eq4, eq5), (r2, r3, x3, r4, x4));

2 組の解が得られるが,2 番目のものが適解である。

res[2][1] |> println
res[2][4] |> println
res[2][5] |> println

   r1*(29*sqrt(5 + 5*sqrt(5)) + 65*sqrt(1 + sqrt(5)) + 125 + 56*sqrt(5))/(2*(40*sqrt(1 + sqrt(5)) + 18*sqrt(5)*sqrt(1 + sqrt(5)) + 105 + 47*sqrt(5)))
   -3*sqrt(5)*r1*sqrt(5 + 5*sqrt(5))/5 - sqrt(5)*r1*sqrt(1 + sqrt(5)) - r1/4 + 3*sqrt(5)*r1/20 + 9*r1*sqrt(1 + sqrt(5))/4 + 27*r1*sqrt(5 + 5*sqrt(5))/20
   -sqrt(5)*r1/10 + r1*sqrt(1 + sqrt(5))/(2*sqrt(5) + 5) + r1/2

円柱の底面(円)の直径が 35 寸 9 分のとき,円柱の高さは 49.0003 寸,小球の個数は 10 個で,小球の半径は 2.58655 寸,小球の中心を通る円の半径は 8.37025 寸である。

その他のパラメータは以下のとおりである。

r2 = 12.2501;  r3 = 46.4137;  x3 = 57.3705;  r4 = 2.58655;  x4 = 8.37025

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 35.9/2
   u = √5
   t = sqrt(1 + u)
   s = u*t
   (r2, r3, x3, r4, x4) = (
       359(29s + 65t + 125 + 56u)/(40(18s + 40t + 47u + 105)),
       359(30254s + 67650t + 58301u + 130365)/(200(987s + 2207t + 2584u + 5778)),
       359(987s + 2207t + 4253 + 1902u)/(10(521s + 1165t + 1364u + 3050)),
       (2513s -5385t + 1077u - 1795)/400,
       -359u/200 + 359t/(20(2u + 5)) + 359/40)
   h = 2r2
   n = 10
   @printf("円柱の高さ = %g;  小球の個数 = %g;  小球の半径 = %g\n", 4r2, n, r4)
   @printf("r2 = %g;  r3 = %g;  x3 = %g;  r4 = %g;  x4 = %g\n", r2, r3, x3, r4, x4)
   θ = atand(h, x3 - r1)
   plot()
   circle22(0, r2, r2, :blue)
   circle(x3, 0, r3, beginangle=180 - θ, endangle=180 + θ)
   circle(x4, 0, r4, :green)
   segment(-r1, h, r1, h, :black)
   segment(r1, h, x3, 0)
   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, r2, " 大球:r2\n (0, r2)", :blue, :left, :vcenter)
       point(x4, 0, "   小球:r4,(x4, 0)", :green, :left, :vcenter)
       point(r1, h, "底面:r1,(0,2r2)", :black, :center, :bottom, delta=delta)
       point(x3, 0, "弧環:r3,(x3,0)", :red, :right, delta=-delta/2)
   end

   function rotate2(ox, oy, r, color=:red; angle=120, beginangle=0, endangle=360, by=0.5, n=0)
       for deg in 0:angle:360-1
           (ox2, oy2) = [cosd(deg) -sind(deg); sind(deg) cosd(deg)] * [ox; oy]
           @printf("(x-%g)^2+(y-%g)^2+z^2=%g^2\n", ox2, oy2, r)
           circle(ox2, oy2, r, color; beginangle, endangle, by, n)
           beginangle += angle
           endangle += angle
       end
   end;
    pyplot(size=(300, 300), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   circle(0, 0, r1, :black, lw=0.1)
   circle(0, 0, r2, :blue, lw=0.1)
   circle(0, 0, x4, :orange)
   rotate2(x4, 0, r4, :green, angle=2*18)
   point(x4, 0, "   小球:r4\n   (x4, 0)", :green, :left, :vcenter)
   point(0, r1, "圓壔", :black, :center, :vcenter, mark=false)
   point(0, r2, "大球", :blue, :center, :vcenter, mark=false)
end;

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

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

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