裏 RjpWiki

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

算額(その134)

2023年02月18日 | Julia

算額(その134)

福島県田村郡三春町上舞木 直毘神社 明治28年(1895)4月
http://www.wasan.jp/fukusima/naobi.html

円の中に甲円,乙円,丙円が入っている。甲円の直径を与えて外円,乙円,丙円の直径を求めよ。

他の場合に準じ,外円の半径を 1 とする。

図のように記号を定め,方程式を立て,解く。
右下にある甲円の中心座標,半径を (x1, y1), r1
右上にある乙円の中心座標,半径を (x2, y2), r2
右下にある丙円の中心座標,半径を (x3, y3), r3

using SymPy

@syms r0::positive, r1::positive, x1::positive, y1::negative, r2::positive, r3::positive

r0 = 1
x2 = 2r2 * cos(PI/6)
y2 = 2r2 * sin(PI/6)
x3 = (1 - r3) * cos(PI/6)
y3 = (r3 - 1) * sin(PI/6)
x1 = (1 - r1) * cos(PI/6)
y1 = (r1 - 1) * sin(PI/6)
eq1 = 2r3 - 1 - r2  # 半径間の関連
eq2 = x1^2 + (1 - r3 - y1)^2 - (r1 + r3)^2  # 甲円と丙円が外接
eq3 = x2^2 + (1 - r3 - y2)^2 - (r3 - r2)^2;  # 乙円が丙円に内接

res = solve([eq2, eq3, eq1], (r1, r2, r3))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (1/3, 1/5, 3/5)

外円の径は,甲円の径の 3 倍である。乙円,丙円の径はそれぞれ甲円の径の 0.6 倍,1.8 倍である。

using Plots
using Printf

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; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
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)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = res[1]
   r0 = 1
   x2 = 2r2 * cos(PI/6)
   y2 = 2r2 * sin(PI/6)
   x3 = (r0 - r3) * cos(PI/6)
   y3 = (r3 - r0) * sin(PI/6)
   x1 = (r0 - r1) * cos(PI/6)
   y1 = (r1 - r0) * sin(PI/6)
   @printf("外円の半径 = %.3f,  甲円の半径 = %.3f,  乙円の半径 = %.3f,  丙円の半径 = %.3f\n", r0, r1, r2, r3)
   @printf("x2 = %.3f,  y2 = %.3f,  x3 = %.3f,  y3 = %.3f\n", x2, y2, x3, y3)
   plot()
   circle(0, 0, r0, :black)
   circle(0, r0 - r1, r1, :green)
   circle(x1, y1, r1, :green)
   circle(-x1, y1, r1, :green)
   circle(x2, y2, r2, :blue)
   circle(-x2, y2, r2, :blue)
   circle(0, -2r2, r2, :blue)
   circle(0, 0, r2, :blue)
   circle(0, r0 - r3, r3)
   circle(x3, y3, r3)
   circle(-x3, y3, r3)
   if more == true
       point(x2, y2, "乙(x2,y2,r2)")
       point(x1, y1, "甲(x1,y1,r1)", :green, :center)
       point(x3, y3, "丙(x3,y3,r3)", :red)
       point(0, r0 - r3, "r0-r3 ", :red, :right)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

   外円の半径 = 1.000,  甲円の半径 = 0.333,  乙円の半径 = 0.200,  丙円の半径 = 0.600
   x2 = 0.346,  y2 = 0.200,  x3 = 0.346,  y3 = -0.200

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

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

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