裏 RjpWiki

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

算額(その950)

2024年05月13日 | Julia

算額(その950)

一四 武州金鑽村 金鑽寺 文化10年(1813)
一五 武州金鑽村 金鑽寺 文化11年(1814)

埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

菱形の中に中央に4円(方併円),左右に2円(縦隔円),上下に2円(横隔円)を入れる。菱形の対角線の長い方が 60 寸,短い方が 32 寸のとき,それぞれの円の直径を求めよ。

菱形の対角線を 2a, 2b; a > b
方併円の半径と中心座標を r1, (r1, r1)
縦隔円の半径と中心座標を r2, (x2, 0)
横隔円の半径と中心座標を r3, (0, y3)
とおき,以下の連立方程式を解く。
eq1, (eq2, eq4), (eq3, eq5) は独立なので,別々に解くことができる。

include("julia-source.txt");

@syms a::positive, b::positive, r1::positive,
     r2::positive, x2::positive,
     r3::positive, y3::positive, c::positive
eq1 = dist2(a, 0, 0, b, r1, r1, r1)
ans_r1 = solve(eq1, r1)[1]
ans_r1 |> println

   a/2 + b/2 - sqrt(a^2 + b^2)/2

# eq2 = dist2(a, 0, 0, b, x2, 0, r2)
eq2 = r2*sqrt(a^2 + b^2) - b*(a - x2)
eq4 = (x2 - r1)^2 + r1^2 - (r1 + r2)^2 |> expand
(ans_r2, ans_x2) = solve([eq2, eq4], (r2, x2))[2]
ans_r2 |> println
ans_x2 |> println

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

# eq3 = dist2(a, 0, 0, b, 0, y3, r3)
eq3 = r3*sqrt(a^2 + b^2) - a*(b - y3)
eq5 = r1^2 + (y3 - r1)^2 - (r1 + r3)^2 |> expand;
(ans_r3, ans_y3) = solve([eq3, eq5], (r3, y3))[2]
ans_r3 |> println
ans_y3 |> println

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

a, b が与えられたとき,r1, x2, r2, y3, r3 を順次求める式は以下のとおりである。

   c = sqrt(a^2 + b^2)
   r1 = a/2 + b/2 - c/2
   x2 = (a^2*r1 - a*b^2 + b^2*r1 - b*r1*c + sqrt(a^4*b^2 - 2a^3*b^2*r1 + 2a^3*b*r1*c + a^2*b^4 + 2a^2*b^2*r1^2 - 2a^2*b*r1^2*c - 2a*b^4*r1 + 2a*b^3*r1*c + 2b^4*r1^2 - 2b^3*r1^2*c))/a^2
   r2 = b*(a - x2)/c
   y3 = (-a^2*b + a^2*r1 - a*r1*c + b^2*r1 + sqrt(a^4*b^2 - 2a^4*b*r1 + 2a^4*r1^2 + 2a^3*b*r1*c - 2a^3*r1^2*c + a^2*b^4 - 2a^2*b^3*r1 + 2a^2*b^2*r1^2 + 2a*b^3*r1*c - 2a*b^2*r1^2*c))/b^2
   r3 = a*(b - y3)/c

a = 60/2, b = 32/2 のとき,中央の4個の円の直径(方併円径),左右にある2個の円の直径(縦隔円径),上下にある2個の円の直径(横隔円径)はそれぞれ,12 寸,12.5237 寸,5.91265 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (60, 32) ./2
   c = sqrt(a^2 + b^2)
   r1 = a/2 + b/2 - c/2
   x2 = (a^2*r1 - a*b^2 + b^2*r1 - b*r1*c + sqrt(a^4*b^2 - 2a^3*b^2*r1 + 2a^3*b*r1*c + a^2*b^4 + 2a^2*b^2*r1^2 - 2a^2*b*r1^2*c - 2a*b^4*r1 + 2a*b^3*r1*c + 2b^4*r1^2 - 2b^3*r1^2*c))/a^2
   r2 = b*(a - x2)/c
   y3 = (-a^2*b + a^2*r1 - a*r1*c + b^2*r1 + sqrt(a^4*b^2 - 2a^4*b*r1 + 2a^4*r1^2 + 2a^3*b*r1*c - 2a^3*r1^2*c + a^2*b^4 - 2a^2*b^3*r1 + 2a^2*b^2*r1^2 + 2a*b^3*r1*c - 2a*b^2*r1^2*c))/b^2
   r3 = a*(b - y3)/c

   @printf("方併円径 = %g;  縦隔円径 = %g;  横隔円径 = %g\n", 2r1, 2r2, 2r3)
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  y3 = %g\n", r1, r2, x2, r3, y3)
   plot([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:black, lw=0.5)
   circle4(r1, r1, r1)
   circle2(x2, 0, r2, :green)
   circle22(0, y3, r3, :blue)
   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", :black, :left, :bottom, delta=delta/2)
       point(0, b, " b", :black, :left, :bottom, delta=delta/2)
       point(r1, r1, "方併円:r1\n(r1,r1)", :red, :center, delta=-delta)
       point(x2, 0, "縦隔円:r2\n(x2,0)", :green, :center, delta=-delta)
       point(0, y3, " 横隔円:r3,(0,y3)", :blue, :left, :vcenter)
   end
end;

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

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

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