裏 RjpWiki

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

算額(その752)

2024年03月05日 | Julia

算額(その752)

東京都稲城市 穴沢天神社 明治10年
山口正義:やまぶき,第15号

https://yamabukiwasan.sakura.ne.jp/ymbk15.pdf

重なった等円と長方形で区切られた領域に大円と小円が入っている。大円と小円の直径がそれぞれ 5 分 7 厘,2 分  1厘のとき,長方形の長辺の長さはいかほどか。

長方形の長辺と短辺の長さを 2b, a
大円の半径と中心座標を r1, (0, r0)
小円の半径と中心座標を r2, (0, r2), (0, a + r2)
等円の半径と中心座標を r0, (r0 - r1, r0)
とおいて以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r0::positive, r1:: poaitive,
     r2::positive, a::positive, b::positive
eq1 = (r0 - r1)^2 + (r0 - r2)^2 - (r0 + r2)^2
eq2 = (r0 - r1)^2 + (a - r0 + r2)^2 - (r0 - r2)^2
eq3 = b - (2r0 - r1)

res = solve([eq1, eq2, eq3], (a, b, r0))

   4-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (r1 - 2*sqrt(r2)*sqrt(r1 + r2) + r2 - sqrt(-4*r1*r2 + 8*r2^(3/2)*sqrt(r1 + r2) - 8*r2^2 + (r1 - 2*sqrt(r2)*sqrt(r1 + r2) + r2)^2), r1 - 4*sqrt(r2)*sqrt(r1 + r2) + 4*r2, r1 - 2*sqrt(r2)*sqrt(r1 + r2) + 2*r2)
    (r1 - 2*sqrt(r2)*sqrt(r1 + r2) + r2 + sqrt(-4*r1*r2 + 8*r2^(3/2)*sqrt(r1 + r2) - 8*r2^2 + (r1 - 2*sqrt(r2)*sqrt(r1 + r2) + r2)^2), r1 - 4*sqrt(r2)*sqrt(r1 + r2) + 4*r2, r1 - 2*sqrt(r2)*sqrt(r1 + r2) + 2*r2)
    (r1 + 2*sqrt(r2)*sqrt(r1 + r2) + r2 - sqrt(-4*r1*r2 - 8*r2^(3/2)*sqrt(r1 + r2) - 8*r2^2 + (r1 + 2*sqrt(r2)*sqrt(r1 + r2) + r2)^2), r1 + 4*sqrt(r2)*sqrt(r1 + r2) + 4*r2, r1 + 2*sqrt(r2)*sqrt(r1 + r2) + 2*r2)
    (r1 + 2*sqrt(r2)*sqrt(r1 + r2) + r2 + sqrt(-4*r1*r2 - 8*r2^(3/2)*sqrt(r1 + r2) - 8*r2^2 + (r1 + 2*sqrt(r2)*sqrt(r1 + r2) + r2)^2), r1 + 4*sqrt(r2)*sqrt(r1 + r2) + 4*r2, r1 + 2*sqrt(r2)*sqrt(r1 + r2) + 2*r2)

4 組の解が得られるが,最後のものが適解である。

i = 4
@printf("a = %g;  b = %g;  r0 = %g\n",
       res[i][1](r1 => 0.57/2, r2 => 0.21/2),
       res[i][2](r1 => 0.57/2, r2 => 0.21/2),
       res[i][3](r1 => 0.57/2, r2 => 0.21/2))

   a = 1.29841;  b = 1.51444;  r0 = 0.899722

長方形の長辺の長さは 2(r1 + 4*sqrt(r2)*sqrt(r1 + r2) + 4*r2) = 3.028888507587845である。
算額の答では「直長三寸〇九厘二九有奇」とあるが,写し間違いであろう。

r1 = 0.57/2
r2 = 0.21/2
2(r1 + 4*sqrt(r2)*sqrt(r1 + r2) + 4*r2)

   3.028888507587845

その他のパラメータは以下のとおり。

r1 = 0.285;  r2 = 0.105;  a = 1.29841;  b = 1.51444;  r0 = 0.899722

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (57, 21) .// 200
   a = 1.29841;  b = 1.51444;  r0 = 0.899722
   @printf("長方形の長辺の長さ = %g;  r1 = %g;  r2 = %g;  a = %g;  b = %g;  r0 = %g\n", 2b, r1, r2, a, b, r0)
   plot([b, b, -b, -b, b], [0, a, a, 0, 0], color=:black, lw=0.5)
   circle2(r0 -  r1, r0, r0)
   circle(0, r0, r1, :blue)
   circle(0, r2, r2, :green)
   circle(0, a + r2, 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(r0 - r1, r0, "等円:r0,(r0-r1,r0)", :red, :left, delta=-delta)
       point(0, r0, "大円:r1\n(0,r1)", :blue, :center, delta=-delta)
       point(0, r2, "   小円:r2,(0,r2)", :green, :left, :bottom, delta=delta)
       point(0, a + r2, "   小円:r2,(0,a+r2)", :green, :left, :bottom, delta=delta)
       point(0, a, "a", :black, :center, delta=-delta)
       point(b, 0, " b", :black, :left, :bottom, delta=delta)
   end
end;


コメント    この記事についてブログを書く
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする
« 算額(その751) | トップ | 算額(その753) »
最新の画像もっと見る

コメントを投稿

Julia」カテゴリの最新記事