裏 RjpWiki

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

算額(その496)

2023年11月16日 | Julia

算額(その496)

岐阜県郡上市八幡町 八幡神社 嘉永3年(1850)

早坂平四郎: 算額の一考察,苫小牧工業専門学校紀要
https://www.tomakomai-ct.ac.jp/wp01/wp-content/uploads/2014/06/kiyou5-8.pdf

外円内に半径の等しい 6 個の大円が入っている。小円は外円の同心円で,この 6 個の大円に外接する。2個の大円に外接し外円に内接する円を甲円,2個の大円と小円に外接する円を乙円とする。
甲円と乙円の直径が与えられたとき,外円の直径を求めよ。

外円の半径と中心座標を r0, (0, 0)
大円の半径と中心座標を r1, (0, r0 - r1), (x1, y1)
小円の半径と中心座標を r2, (0, 0); r2 = r0 - 2r1
甲円の半径と中心座標を r3, (x3, y3)
乙円の半径と中心座標を r4, (x4, y4)
として,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive,
     r3::positive, r4::positive;
@syms r0, r1, r2, r3, r4

r2 = r0 - 2r1
x1 = (r0 - r1)cosd(Sym(30))
y1 = (r0 - r1)sind(Sym(30))
x3 = (r0 - r3)cosd(Sym(60))
y3 = (r0 - r3)sind(Sym(60))
x4 = (r2 + r4)cosd(Sym(60))
y4 = (r2 + r4)sind(Sym(60))

eq1 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2
eq2 = (x1 - x4)^2 + (y1 - y4)^2 - (r1 + r4)^2

res = solve([eq1, eq2], (r0, r1))

   2-element Vector{Tuple{Sym, Sym}}:
    ((-sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))*(-56*sqrt(3)*r3^2 + 97*r3^2 - 1448*sqrt(3)*r3*r4 + 2508*r3*r4 - 1616*sqrt(3)*r4^2 + 2799*r4^2 + (-sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))*(-2702*r3 + 1560*sqrt(3)*r3 - 2702*r4 + 1560*sqrt(3)*r4))/(-1351*r3^2 + 780*sqrt(3)*r3^2 - 780*sqrt(3)*r4^2 + 1351*r4^2), -sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))
    ((sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))*(-56*sqrt(3)*r3^2 + 97*r3^2 - 1448*sqrt(3)*r3*r4 + 2508*r3*r4 - 1616*sqrt(3)*r4^2 + 2799*r4^2 + (sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))*(-2702*r3 + 1560*sqrt(3)*r3 - 2702*r4 + 1560*sqrt(3)*r4))/(-1351*r3^2 + 780*sqrt(3)*r3^2 - 780*sqrt(3)*r4^2 + 1351*r4^2), sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))

2組の解が得られるが,最初のものが適解である。

甲円,乙円の直径が 8, 2 のとき,外円の直径は約 86.9559 である。
外円の直径 = 86.9559;  r0 = 43.4779;  r1 = 17.3042;  r2 = 8.86948;  r3 = 4,  r4 = 1

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r3, r4) = (8, 2) .// 2
   (r0, r1) = ((-sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))*(-56*sqrt(3)*r3^2 + 97*r3^2 - 1448*sqrt(3)*r3*r4 + 2508*r3*r4 - 1616*sqrt(3)*r4^2 + 2799*r4^2 + (-sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))*(-2702*r3 + 1560*sqrt(3)*r3 - 2702*r4 + 1560*sqrt(3)*r4))/(-1351*r3^2 + 780*sqrt(3)*r3^2 - 780*sqrt(3)*r4^2 + 1351*r4^2), -sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))
   r2 = r0 - 2r1
   x1 = (r0 - r1)cosd(30)
   y1 = (r0 - r1)sind(30)
   x3 = (r0 - r3)cosd(60)
   y3 = (r0 - r3)sind(60)
   x4 = (r2 + r4)cosd(60)
   y4 = (r2 + r4)sind(60)
   @printf("外円の直径 = %g;  r0 = %g;  r1 = %g;  r2 = %g;  r3 = %g,  r4 = %g\n", 2r0, r0, r1, r2, r3, r4)    
   plot()
   circle(0, 0, r0, :black)
   circle(0, 0, r2, :blue)
   rotate(0, r0 - r1, r1, :green, angle=60)
   rotate(r0 - r3,0, r3, :magenta, angle=60)
   rotate(r2 + r4, 0, r4, :red, angle=60)
   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, r0, " r0", :black, :left, :bottom, delta=delta/2)
       point(0, r2, " r2", :blue, :left, :bottom, delta=delta/2)
       point(0, r0 - r1, " r0-r1", :green, :left, :vcenter)
       point(x1, y1, "    大円:r1,(x1, y1)", :green, :center, :bottom, delta=delta/2)
       point(x3, y3, "甲円:r3 \n(x3, y3) ", :black, :right, :vcenter)
       point(x4, y4, " 乙円:r4,(x4, y4) ", :black, :left, :vcenter)
   else
       plot!(showaxis=false)
   end
end;

 


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

1 コメント

コメント日が  古い順  |   新しい順
ありがとうございます (r-de-r)
2023-11-17 20:59:25
いつも応援ありがとうございます。頑張りますので,今後とも宜しくお願い申し上げます。
返信する

コメントを投稿

Julia」カテゴリの最新記事