裏 RjpWiki

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

算額(その163)

2023年03月14日 | Julia

算額(その163)

岐阜県大垣市西外側町 大垣八幡神社 天保年間
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

第14問: 2 個の紫円を黃円に入れ,その間に 4 個の赤円を入れる。紫円の直径を知って黃円の直径を求めよ。

一般性を失わず黃円の半径を 1 とする。紫円,赤円の径をそれぞれ r1, r2 として連立方程式を解く。

using SymPy

@syms r1::positive, r2::positive;

eq1 = (1 - r2)^2 + (1 - r1)^2 - (r1 + r2)^2
eq2 = r2^2 + (1 - r1)^2 - (r1 - r2)^2

res = solve([eq1, eq2], (r1, r2))
res |> println

   Tuple{Sym, Sym}[(1/8 + sqrt(17)/8, 5/4 - sqrt(17)/4)]

r1 = res[1][1]
r2 = res[1][2];

1/r1 |> simplify |> println  # 紫円 r1 を知って,黃円の径を得るときの倍数

   -1/2 + sqrt(17)/2

r2/r1 |> simplify |> println  # 紫円 r1 を知って,赤円 r2 の径を得るときの倍数

   -11/4 + 3*sqrt(17)/4

1/r2 |> simplify |> println  # 赤円 r2 を知って,黃円の径を得るときの倍数

   sqrt(17)/2 + 5/2

r1/r2 |> simplify |> println  # 赤円 r2 を知って,紫円 r1 の径を得るときの倍数

   11/8 + 3*sqrt(17)/8

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
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 draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = res[1]
   println("r1 = $r1 = $(r1.evalf())\nr2 = $r2 = $(r2.evalf())")
   plot()
   circle(0, 0, 1, :yellow, fill=true)
   circle(0, 1 - r1, r1, :orchid1, fill=true)
   circle(0, r1 - 1, r1, :orchid1, fill=true)
   circle(r2, 0, r2, :tomato1, fill=true)
   circle(-r2, 0, r2, :tomato1, fill=true)
   circle(1 - r2, 0, r2, :tomato1, fill=true)
   circle(r2 - 1, 0, r2, :tomato1, fill=true)
   circle(0, 1 - r1, r1, :orchid3)
   circle(0, r1 - 1, r1, :orchid3)
   if more == true
       point(0, 1 - r1, " 1-r1", :black)
       point(r2, 0, "r2", :black, :top)
       point(1 - r2, 0, "1-r2", :black, :top)
       point(1, 0, " 1", :black, :left, :top)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

   r1 = 1/8 + sqrt(17)/8 = 0.640388203202208
   r2 = 5/4 - sqrt(17)/4 = 0.219223593595585

 

 

 

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

算額(その162)

2023年03月14日 | Julia

算額(その162)

岐阜県大垣市西外側町 大垣八幡神社 天保年間
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

第13問: 長方形内に正三角形と青円 2 個,赤円を入れる。青円の直径を知って,赤円の径を求めよ。

青円の半径とその中心の x 座標を r1,x1 として方程式を解く。r1 = 1 として一般性を失わない。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
   p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
   l = sympy.Line(p1, p2)
   l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す!!
end;

@syms r1::positive, r2::positive, x1::positive;
r1 = 1
eq1 = x1^2 + (r1-r2)^2 - (r1+r2)^2  # 青円と赤円が外接する
eq2 = distance((2r1 - 2r2)/sqrt(Sym(3)), 0, 0, 2r1 - 2r2, x1, r1) - r1^2  # 青円が正三角形の斜辺に接する

res = solve([eq1, eq2], (x1, r2))
res |> println

   Tuple{Sym, Sym}[(sqrt(3)*(-123*sqrt(3) - 60*(3 - 3*sqrt(3)/2)^2 + 8*(3 - 3*sqrt(3)/2)^3 + 231)/12, 3 - 3*sqrt(3)/2)]

res[1][1] |> simplify |> println  # x1
res[1][2] |> simplify |> println  # r2

   3 - sqrt(3)
   3 - 3*sqrt(3)/2

赤円の径は青円の径の 3 - 3*√3/2 倍である。

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
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 draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1
   (x1, r2) = ((3 - √3)r1, (3 - 3*√3/2)r1)
   println("青円径 = $r1, 赤円径 = $r2;  x1 = $x1")
   plot([x1 + r1, x1 + r1, -x1 - r1, -x1 - r1, x1 + r1], [0, 2r1, 2r1, 0, 0], color=:black, lw=0.5)
   circle(0, 2r1 - r2, r2, :orangered, fill=true)
   circle(x1, r1, r1, :steelblue1, fill=true)
   circle(-x1, r1, r1, :steelblue1, fill=true)
   plot!([(2r1-2r2)/√3, 0, -(2r1-2r2)/√3, (2r1-2r2)/√3], [0, 2r1 - 2r2, 0, 0],
         linecolor=:olivedrab1, linewidth=0.5, seriestype=:shape, fillcolor=:olivedrab1)
   if more == true
       point(x1, r1, "C:(x1,r1)", :black, :center)
       point(0, 2r1 - r2, " 2r1-r2", :black)
       point(0, 2r1 - 2r2, " A:2r1-2r2", :black)
       point((2r1-2r2)/√3, 0, " B:(2r1-2r2)/√3", :black, :center, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

   青円径 = 1, 赤円径 = 0.401923788646684;  x1 = 1.2679491924311228

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

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

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