裏 RjpWiki

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

算額(その936)

2024年05月09日 | Julia

算額(その936)

群馬県前橋市河原浜町(旧大胡町) 大胡神社 大正4年(1915)
「算額」第三集 全国調査,香川県算額研究会

大円7個が互いに外接し合っている。大円の直径が与えられたとき,小円の直径を求めよ。

算額(その653)では,外円があると仮定して解いた。

算額の図では小円として緑色の 6 個の円が描かれているが,それでは問題にならない。おそらくは図の青い 6 個のことを指しているか,もしくは両方が描いてあって「小円 12 個」という意図だったのかもしれない。

大円の半径と中心座標を r1, (2r1\*cos(π/6), 2r1*sin(π/6))
小円の半径と中心座標を r2, (r1 + r2, 0)
とおき,以下の方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive
eq = (r1√Sym(3) - (r1 + r2))^2 + r1^2 - (r1 + r2)^2
r2 = solve(eq, r2)[1]
r2 |> println

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

小円の半径は大円の半径の 2√3/3 - 1 倍である。
大円の直径が 1 寸のとき,小円の直径は (2√3 - 3)/3 = 0.15470053837925146 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1/2
   r2 = r1*(2√3/3 - 1)
   @printf("大円の直径が %g のとき,小円の直径は %g である。\n", 2r1, 2r2)
   plot()
   circle(0, 0, r1)
   rotate(0, 2r1, r1, angle=60)
   rotate(r1 + r2, 0, r2, :blue, angle=60)
   rotate(2r1√3 - r1 - r2, 0, r2, :green, angle=60)
   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(r1√3, r1, "大円:r1,(r1√3, r1)", :red, :center, delta=-delta/2)
       point(r1 + r2, 0, "小円:r2\n(r1+r2,0)", :blue, :left, delta=-2delta)
   end
end;

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

算額(その935)

2024年05月09日 | Julia

算額(その935)

栃木県真岡市東郷(注) 大前神社 絵馬堂 紀元2536年(明治9年(1876))
「算額」第三集 全国調査,香川県算額研究会

注: 住所は,真岡市荒井町となっていたが,現在は真岡市東郷か。

菱形の頂点と各辺の中点を線分で結び,区画された領域に甲円と乙円を入れる。甲円と乙円の直径を与えたとき,菱形の一辺の長さを求めよ。

菱形の対角線の長さを,長い方を 2a,短い方を 2b
甲円の半径と中心座標を r1, (0, y)
乙円の半径と中心座標を r2, (x, 0)
とおき,以下の連立方程式を解いて,a, b, x, y を求める。
菱形の一辺の長さは sqrt(a^2 + b^2) である。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, x::positive, y::positive,
     r1::positive, r2::positive
eq1 = dist2(a, 0, -a/2, b/2, 0, y, r1)
eq2 = dist2(a/2, -b/2, 0, b, 0, y, r1)
eq3 = dist2(a, 0, -a/2, b/2, x, 0, r2)
eq4 = dist2(a/2, -b/2, 0, b, x, 0, r2);
# res = solve([eq1, eq2, eq3, eq4], (a, b, x, y))

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")

   a^2*b^2 - 6*a^2*b*y - 9*a^2*r1^2 + 9*a^2*y^2 - b^2*r1^2,  # eq1
   a^2*b^2 - 2*a^2*b*y - a^2*r1^2 + a^2*y^2 - 9*b^2*r1^2,  # eq2
   a^2*b^2 - 9*a^2*r2^2 - 2*a*b^2*x - b^2*r2^2 + b^2*x^2,  # eq3
   a^2*b^2 - a^2*r2^2 - 6*a*b^2*x - 9*b^2*r2^2 + 9*b^2*x^2,  # eq4

甲円が斜線に接するという条件式 eq1 を解いて y を求め,eq2 に代入し y を消去する。

ans_y = solve(eq1, y)[2]
ans_y |> println

   (a*b + r1*sqrt(9*a^2 + b^2))/(3*a)

eq2 = eq2(y => ans_y)*(9/4b) |> simplify;

乙円が斜線に接するという条件式 eq3 を解いて x を求め,eq4 に代入し x を消去する。

ans_x = solve(eq3, x)[1]
ans_x |> println

   a - r2*sqrt(9*a^2 + b^2)/b

eq4 = eq4(x => ans_x)/4a |> simplify;

eq2 を解いて b を求め eq4 に代入し b を消去する。

ans_b = solve(eq2, b)[2] |> simplify
ans_b |> println

   3*a^2*r1*sqrt(1/((a - 5*r1)*(a - 4*r1)))/sqrt(a^2 + 9*a*r1 + 20*r1^2)

eq4 = eq4(b => ans_b) |> simplify;

eq4 を解いて a を求める。

ans_a = solve(eq4, a)[5]
ans_a |> println

   4*sqrt(5)*r1*sqrt(r2)*sqrt(-1/(3*r1 - 5*r2))

ans_b, ans_x, ans_y は ans_a を逆順に代入していけばよい。
追加的に手動で式を整理して,以下のような結果を得る。

(r1, r2) = (114, 99)
a = 4r1*sqrt(5r2/(5r2 - 3r1))
b = 3a^2*r1/sqrt(a^4 - 41a^2*r1^2 + 400r1^4)
x = a - r2*sqrt(9*a^2 + b^2)/b
y = (a*b + r1*sqrt(9a^2 + b^2))/(3a)
(a, b, x, y)

   (820.2037049703316, 572.2045280998937, 383.121467453247, 307.7766779931246)

甲円,乙円の半径が 114 寸,99 寸のとき,菱形の一辺の長さは 1000.0760669194523 寸である。先に示した図は,この条件のもとの図である。

sqrt(a^2 + b^2)

   1000.0760669194523

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (114, 99)
   a = 4r1*sqrt(5r2/(5r2 - 3r1))
   b = 3a^2*r1/sqrt(a^4 - 41a^2*r1^2 + 400r1^4)
   x = a - r2*sqrt(9a^2 + b^2)/b
   y = (a*b + r1*sqrt(9a^2 + b^2))/(3a)
   @printf("甲円の直径 = %g,乙円の直径 = %g のとき,菱形の一辺の長さは %g\n", 2r1, 2r2, sqrt(a^2 + b^2))
   plot([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:magenta, lw=0.5)
   plot!([0, a/2, -a, a/2, 0], [-b, b/2, 0, -b/2, b], color=:blue, lw=0.5)
   plot!(-[0, a/2, -a, a/2, 0], [-b, b/2, 0, -b/2, b], color=:blue, lw=0.5)
   circle22(0, y, r1) 
   circle2(x, 0, r2, :green) 
   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(a, 0, " a", :magenta, :left, :bottom, delta=delta/2)
       point(0, b, " b", :magenta, :left, :bottom, delta=delta/2)
       point(a/2, b/2, " (a/2,b/2)", :magenta, :left, :bottom, delta=delta/2)
       point(0, y, " 甲円:r1,(0,y)", :black, :left, :vcenter)
       point(x, 0, " 乙円:r2,(x,0)", :black, :left, :vcenter)
   end
end;

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

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

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