裏 RjpWiki

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

算額(その1081)

2024年06月21日 | Julia

算額(その1081)

九十八 岩手県江刺市 雨宝堂 現中善観音堂 文政10年(1827)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

キーワード:円11個,外円,二等辺三角形(直角三角形)

全円の中に圭(二等辺三角形),甲円 3 個,乙円 1 個,丙円 2 個,丁円 4 個を容れる。丙円の直径が 1 寸のとき,乙円の直径はいかほどか。

注:「問」では圭(二等辺三角形)といっているが底辺は全円の直径である。したがって,三角形は二等辺直角三角形である。

全円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1), (r1, -r1)
乙円の半径と中心座標を r2, (0, r2 - R)
丙円の半径と中心座標を r3, (x3, x3)
丁円の半径と中心座標を r4, (x4, y4), (y4, x4)
とおき,以下の連立方程式を解く。

include("julia-source.txt")

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

r1 = res[1] |> sympy.sqrtdenest |> simplify
r2 = res[2] |> sympy.sqrtdenest |> simplify
x3 = res[3] |> sympy.sqrtdenest |> simplify
R = res[4]
(r1, r2, x3, R)

   (2*sqrt(2)*r3, 2*r3*(2 - sqrt(2)), r3*(4 + 3*sqrt(2))/2, 2*r3*(sqrt(2) + 2))

乙円の半径は丙円の半径の 4 - 2√2 = 1.17157287525381 倍である。
丙円の直径が 1 寸のとき,乙円の直径は 1.17157287525381 寸である。

---

これ以降で,図を描くために丁円の半径と中心座標を求める。

パラメータを簡単にするために図を反時計回りに 45° 回転させたものを使う。
まず丁円の半径を求める。算法助術の公式29 を使うまでもなく,以下で求めることができる。中心座標は後で求めるので使わない(使用するためには回転しなければならずかえって面倒になる)。

using SymPy
@syms R::positive, r3::positive, r4::positive, x4::positive, y4::positive
y4 = R - 2r3 + r4
eq5 = x4^2 + y4^2 - (R - r4)^2
eq6 = x4^2 + (r3 - r4)^2 - (r3 + r4)^2
res2 = solve([eq5, eq6], (r4, x4))[1]

   (-r3*(-R + r3)/R, 2*r3*sqrt(R - r3)/sqrt(R))

@syms x4::positive, y4::positive
x3 = r3*(4 + 3√Sym(2))/2
R = 2r3*(√Sym(2) + 2)
r4 = -r3*(-R + r3)/R
eq7 = (x4 - x3)^2 + (y4 - x3)^2 - (r3 + r4)^2 |> expand |> simplify
eq8 = x4^2 + y4^2 - (R - r4)^2 |> expand |> simplify
res3 = solve([eq7, eq8], (x4, y4))[1]

   (-r3*sqrt(47321*sqrt(2) + 66922)/(140 + 99*sqrt(2)) + 5*sqrt(2)*r3/4 + 9*r3/4, r3*sqrt(47321*sqrt(2) + 66922)/(140 + 99*sqrt(2)) + 5*sqrt(2)*r3/4 + 9*r3/4)

丙円の直径が 1 寸のときの描画パラメータは以下のとおりである。

r3= 0.5;  r1 = 1.41421;  r2 = 0.585786;  x3 = 2.06066;  R = 3.41421
r4 = 0.42677669529663687

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 1/2
   (r1, r2, x3, R) = (2*sqrt(2)*r3, 2*r3*(2 - sqrt(2)), r3*(4 + 3*sqrt(2))/2, 2*r3*(sqrt(2) + 2))
   r4 = -r3*(-R + r3)/R
   (x4, y4) = (-r3*sqrt(47321*sqrt(2) + 66922)/(140 + 99*sqrt(2)) + 5*sqrt(2)*r3/4 + 9*r3/4, r3*sqrt(47321*sqrt(2) + 66922)/(140 + 99*sqrt(2)) + 5*sqrt(2)*r3/4 + 9*r3/4)
   @printf("丙円の直径が %g のとき,乙円の直径は %g である。\n", 2r3, 2r2)
   @printf("r3= %g;  r1 = %g;  r2 = %g;  x3 = %g;  R = %g\n", r3, r1, r2, x3, R)
   plot([0, R, -R, 0], [R, 0, 0, R], color=:hotpink, lw=0.5)
   circle(0, 0, R, :blue)
   circle(0, r1, r1)
   circle2(r1, -r1, r1)
   circle(0, r2 - R, r2, :magenta)
   circle2(x3, x3, r3, :darkred)
   circle2(x4, y4, r4, :green)
   circle2(y4, x4, r4, :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(0, R, " R", :blue, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
       point(0, r1, "甲円:r1,(0,r1)", :red, :center, delta=-delta/2)
       point(r1, -r1, "甲円:r1,(r1,-r1)", :red, :center, delta=-delta/2)
       point(0, r2 - R, "乙円:r2\n(0,r2-R)", :magenta, :center, :bottom, delta=delta/2)
       point(x3, x3, " 丙円:r3,(x3,x3)", :darkred, :center, delta=-delta/2, deltax=-4delta)
       point(x4, y4, " 丁円:r4,(x4,y4)", :green, :center, :bottom, delta=delta, deltax=-4delta)
       point(y4, x4, " 丁円:r4,(y4,x4)", :green, :center, delta=-delta/2, deltax=-8delta)
   end
end;


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

コメントを投稿

Julia」カテゴリの最新記事