裏 RjpWiki

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

算額(その940)

2024年05月10日 | ブログラミング

算額(その940)

福岡県朝倉市上秋月(旧甘木市秋月町) 秋月八幡宮絵馬堂 明治7年(1874)
「算額」第三集 全国調査,香川県算額研究会

外円に内接する正五角形に一本の対角線を引き,中円と小円を入れる。中円と小円の直径が与えられたとき,対角線の長さを求めよ。

外円の大きさを確定するために必要なのは中円の直径だけであり,小円は無関係である。

算額(その937)の一部分で,求めるものが違うだけである。

座標の定義に使う三角関数の値
(x18, y18) = (cosd(18), sind(18))
(x54, y54) = (cosd(54), sind(54))

正五角形が内接する外円の半径と中心座標を R, (0, 0)
中円の半径と中心座標を r1,(2r1*x18, 2r1*y18)
とおき,以下の方程式より R を得る。

外円の半径は中円の半径の (1 + √5) 倍である。

include("julia-source.txt");

@syms x18::positive, y18::positive, x54::positive, y54::positive
(x18, y18) = (cosd(Sym(18)), sind(Sym(18)))
(x54, y54) = (cosd(Sym(54)), sind(Sym(54)));

@syms R::positive, r1::positive
eq1 = dist2(0, R, R*x54, -R*y54, 2r1*x18, 2r1*y18, r1)
R = solve(eq1, R)[1] |> sympy.sqrtdenest |> simplify
R |> println

   r1*(1 + sqrt(5))

外円の半径が R のとき,求める対角線の長さは RE = AB = 2R*cosd(18)である。

最終的に,

中円の半径が 6.5 のとき,R = 6.5*(1 + sqrt(5)),RE = 2 * 6.5*(1 + sqrt(5)) * cosd(18) = 40.00988598327829 である。

2 * 6.5*(1 + sqrt(5)) * cosd(18)

   40.00988598327829

function draw(more=false)
   pyplot(size=(600, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x18, y18) = (cosd(18), sind(18))
   (x54, y54) = (cosd(54), sind(54))
   r1 = 6.5
   (x1, y1) = 2r1.*(x18, y18)
   R = r1*(1 + √5)
   l = 2 * r1*(1 + sqrt(5)) * cosd(18)
   println("中円の直径が $(2r1) のとき,対角線の長さは $l")
   x = R .* [x18, 0, -x18, -x54, x54, x18]
   y = R .* [y18, 1, y18, -y54, -y54, y18]
   plot(x, y, color=:blue, lw=0.5)
   circle(0, 0, R, :gray70)
   segment(0, R, R*x54, -R*y54, :green)
   circle(2r1*x18, 2r1*y18, r1)
   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(R*x18, R*y18, " A", :black, :left, :vcenter)
       point(-R*x18, R*y18, "B ", :black, :right, :vcenter)
       point(R*x54, -R*y54, " E", :black, :left, :vcenter)
       point(0, R, " R", :blue, :left, :bottom, delta=delta/2)
       point(R*x18, R*y18, "(R*x18,R*y18)  ", :blue, :right, :bottom, delta=delta/2)
       point(R*x54, -R*y54, " (R*x54,-R*y54)   ", :blue, :right, :bottom, delta=delta/2)
       point(2r1*x18, 2r1*y18, "甲円:r1\n(2r1*x18,2r1*y18)", :red, :center, delta=-delta/2)
       plot!(xlims=(-1.1R*x18, 1.1R*x18))
   end
end;

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

算額(その939)

2024年05月10日 | Julia

算額(その939)

福岡県朝倉市上秋月(旧甘木市秋月町) 秋月八幡宮 明治4年(1871)
「算額」第三集 全国調査,香川県算額研究会

半円の中に大円,大円の中に 8 個の小円を入れる。半円の直径が 2 寸のとき,小円の直径はいかほどか。

半円は重要ではない。

大円と小円の半径を r1, r2 とおくと,(r1 - r2)sin(π/8) = r2 という関係が成り立つ。
方程式を解いて,小円の半径は大円の半径の sqrt(2 - √2)/(sqrt(2 - √2) + 2) 倍である。
大円の直径が 1 寸のとき,小円の直径は 0.2767686539141552 

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive
eq1 = (r1 - r2)sin(PI/8) - r2
r2 = solve(eq1, r2)[1]
r2 |> println

   r1*sqrt(2 - sqrt(2))/(sqrt(2 - sqrt(2)) + 2)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1/2
   r2 = r1*sqrt(2 - √2)/(sqrt(2 - √2) + 2)
   @printf("半円の直径が %g 寸のとき小円の直径は %g 寸である。\n", 4r1, 2r2)
   plot()
   circle(0, 0, 2r1, :blue, beginangle=0, endangle=180)
   circle(0, r1, r1)
   for i = 1:8
       theta = pi*(i/4 - 1/8)
       circle((r1 - r2)cos(theta), r1 + (r1 - r2)sin(theta), r2, :green)
   end
   plot!((r1 - r2) .* [0, cos(pi/8), cos(pi/8), 0], r1 .+ (r1 - r2) .* [0, 0, sin(pi/8), 0], color=:black, lw=0.5)
   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(2r1, 0, "2r1 ", :blue, :right, :bottom, delta=delta/2)
       point(0, r1, "大円:r1,(0,r1)", :red, :center, delta=-delta/2)
       
   end
end;

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

算額(その938)

2024年05月10日 | Julia

算額(その938)

福岡県朝倉市上秋月(旧甘木市秋月町) 秋月八幡宮 明治4年(1871)
「算額」第三集 全国調査,香川県算額研究会

外円の中に大円,中円,小円が入っている。外円の直径が 5 寸のとき,中円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, R - r1)
中円の半径と中心座標を r2, (0, R - r2)
小円の半径と中心座標を r3, (x3, r3)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive,
     r3::positive, x3::positive
eq1 = x3^2 + r3^2 - (R - r3)^2
eq2 = x3^2 + (R - r1 - r3)^2 - (r1 + r3)^2
eq3 = 2r2 + r3 - R
eq4 = 2r1 -r3 - R
solve([eq1, eq2, eq3, eq4], (r1, r2, r3, x3))

   1-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (3*R/5, 2*R/5, R/5, sqrt(15)*R/5)

中円の半径 r2 は,外円の半径 R の 2/5 倍である。
外円の直径が 5 寸のとき,中円の直径は 2 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 5/2
   (r1, r2, r3, x3) = R .* (3/5, 2/5, 1/5, √15/5)
   @printf("外円の直径が %g 寸のとき,中円の直径は %g 寸である。\n", 2R, 2r2)
   plot()
   circle(0, 0, R, :blue)
   circle22(0, R - r1, r1)
   circle22(0, R - r2, r2, :green)
   circle(0, 0, r3, :brown)
   circle4(x3, r3, r3, :brown)
   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, R - r1, "大円:r1,(0,R-r1)", :red, :center, delta=-delta/2)
       point(0, R - r2, "中円:r2,(0,R-r2)", :green, :center, delta=-delta/2)
       point(x3, r3, "小円:r3\n(x3,r3)", :brown, :center, delta=-delta/2)
       point(0, R, " R", :blue, :left, :bottom, delta=delta/2)
   end
end;

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

算額(その937)

2024年05月10日 | Julia

算額(その937)

千葉県富津市西大和田(旧大佐和町絹) 吾妻神社 文久2年(1862)
「算額」第三集 全国調査,香川県算額研究会

正五角形の中に対角線 2 本を引き,区分された領域に 甲円 2 個,乙円 2 個,丙円 1 個を入れる。甲円,乙円の直径がそれぞれ 1.6180 寸,1.3819 寸のとき,丙円の直径はいかほどか。

座標の定義に使う三角関数の値
(x18, y18) = (cosd(18), sind(18))
(x54, y54) = (cosd(54), sind(54))

正五角形が内接する外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1,(0, R*y18 - r1)
乙円の半径と中心座標を r2, (x2, 2R*y18 - R)
丙円の半径と中心座標を r3, (0, r3 - 2R*y54)
とする。

算額の答えを得るには,実は乙円は関係ない(外円の大きさ R も関係ない)。

2 つの三角形 ABC と DEC は相似で,相似比は R*x18 : R*x54 = cosd(18) : cosd(54) である。

したがって,それに内接する甲円と丙円も相似で,その半径の比も r1:r3 = cosd(18) : cosd(54) なので,
r3 = r1 * cosd(54)/cosd(18)
  = r1 * sqrt(5/8 - sqrt(5)/8) / sqrt(sqrt(5)/8 + 5/8)
  = 1.6180/2 * (-1/2 + sqrt(5)/2)
  = 0.499989496898665

すなわち,甲円の直径が 1.6180 寸のとき,丙円の直径は 0.99997899379733 寸(ほぼほぼ 1 寸)である。

1.6180/2 * (-1/2 + sqrt(5)/2)

   0.499989496898665

以下では,図を描くためにパラメータを得る手順を示す。

include("julia-source.txt");

using SymPy
@syms x18::positive, y18::positive, x54::positive, y54::positive
(x18, y18) = (cosd(Sym(18)), sind(Sym(18)))
(x54, y54) = (cosd(Sym(54)), sind(Sym(54)));

まず,外円の大きさを求める。外円の半径 R は,甲円の半径 r1 の (1 + √5) 倍である。 

@syms R::positive, r1::positive
eq1 = dist2(0, R, R*x18, R*y18, 0, R*y18 + r1, r1)
R = solve(eq1, R)[1] |> sympy.sqrtdenest |> simplify
R |> println

   r1*(1 + sqrt(5))

次に,丙円の大きさを求める。丙円の半径 r3 は(前にも示したが),甲円の半径 r1 の (√5 - 1)/2 倍である。

@syms r3::positive
eq4 = dist2(R*x18, R*y18, -R*x54, -R*y54, 0, r3 - R*y54, r3);
eq4 = eq4 |> simplify |> numerator
r3 = solve(eq4, r3)[2] |> simplify |> sympy.sqrtdenest |> simplify
r3 |> println

   r1*(-1 + sqrt(5))/2

最後に,乙円の中心の x 座標を求める。

@syms r2::positive, x2::positive, y2::negative
eq2 = dist2(R*x18, R*y18, R*x54, -R*y54, x2, 2R*y18 - R, r2) |> simplify |> numerator;
ans_x2 = solve(eq2, x2)[1] |> simplify;
ans_x2 |> println

   r1*(15*sqrt(2)*sqrt(sqrt(5) + 5) + 11*sqrt(10)*sqrt(5 - sqrt(5)) + 25*sqrt(2)*sqrt(5 - sqrt(5)) + 7*sqrt(10)*sqrt(sqrt(5) + 5))/(2*(11*sqrt(5) + 25)) - 2*sqrt(-850*r1^2 - 380*sqrt(5)*r1^2 + 190*r1^2*sqrt(5 - sqrt(5))*sqrt(sqrt(5) + 5) + 85*sqrt(5)*r1^2*sqrt(5 - sqrt(5))*sqrt(sqrt(5) + 5) + 152*sqrt(5)*r2^2 + 340*r2^2)/(11*sqrt(5) + 25)

function draw(more=false)
   pyplot(size=(600, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x18, y18) = (cosd(18), sind(18))
   (x54, y54) = (cosd(54), sind(54))
   r1 = 1.6180/2
   R = r1*(1 + √5)
   r3 = r1*(√5 - 1)/2
   r2 = 1.3819/2
   s = sqrt(5 - √5)
   t = sqrt(√5 + 5)
   x2 = r1*((15√2 + 7√10)t + (11√10 + 25√2)s)/(22√5 + 50) - 2sqrt(-(850 + 380√5)r1^2 + (190 + 85√5)r1^2*s*t + (152√5 + 340)r2^2)/(11√5 + 25)
   y2 = 2R*y18 - R
   @printf("2r1 = %g;  2R = %g;  2r3 = %g\n", 2r1, 2R, 2r3)
   x = R .* [x18, 0, -x18, -x54, x54, x18]
   y = R .* [y18, 1, y18, -y54, -y54, y18]
   plot(x, y, color=:blue, lw=0.5)
   segment(-R*x18, R*y18, R*x18, R*y18, :green)
   segment(-R*x54, -R*y54, R*x18, R*y18, :green)
   segment(R*x54, -R*y54, -R*x18, R*y18, :green)
   circle(0, R*y18 + r1, r1)
   circle(0, R*y18 - r1, r1)
   circle(0, r3 - R*y54, r3, :brown)
   circle2(x2, y2, r2, :magenta)
   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(R*x18, R*y18, " A", :black, :left, :vcenter)
       point(-R*x18, R*y18, "B ", :black, :right, :vcenter)
       point(0, 2R*y18 - R, "C  ", :black, :right, :vcenter)
       point(-R*x54, -R*y54, "D ", :black, :right, :vcenter)
       point(R*x54, -R*y54, " E", :black, :left, :vcenter)
       point(0, R, " R", :blue, :left, :bottom, delta=delta/2)
       point(R*x18, R*y18, "(R*x18,R*y18)  ", :blue, :right, :bottom, delta=delta/2)
       point(R*x54, -R*y54, " (R*x54,-R*y54)   ", :blue, :right, :bottom, delta=delta/2)
       point(0, R*y18 + r1, "甲円:r1\n(0,R*y18+r1)", :red, :center, delta=-delta/2)
       point(0, R*y18 - r1, "甲円:r1\n(0,R*y18-r1)", :red, :center, delta=-delta/2)
       point(x2, 2R*y18 - R, "乙円:r2\n(x2,2R*y18-R)", :magenta, :center, delta=-delta/2)
       point(0, r3 - R*y54, "丙円:r3\n(0,r3-2R*y54)", :brown, :center, :bottom, delta=delta/2)
       plot!(xlims=(-1.1R*x18, 1.1R*x18))
   end
end;

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

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

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