裏 RjpWiki

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

算額(その432)

2023年09月15日 | Julia

算額(その432)

神奈川県横須賀市 浦賀叶大明神 文化11年(1814)

中村信弥(2001):幻の算額
http://www.wasan.jp/maborosi/maborosi.html

大円に等円 2 個が内接し,その上下に甲円と乙円が大円に内接し,等円に外接している。大円の直径が 1 尺,甲円の直径が 5 寸のとき乙円の直径はいかほどか。

大円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (0, r0 - r1)
等円の半径と中心座標を r2, (r2, y2)
乙円の半径と中心座標を r3, (0, r3 - r0)
とし,以下の連立方程式の解を求める。

include("julia-source.txt");

using SymPy
@syms r0::positive, r1::positive, r2::positive,
     y2::negative, r3::positive;
eq1 = r2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2
eq2 = r2^2 + (y2 - r3 + r0)^2 - (r2 + r3)^2
eq3 = r2^2 + y2^2 - (r0 - r2)^2
res = solve([eq1, eq2, eq3], (r2, r3, y2))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (4*r0*r1*(r0 - r1)/(r0 + r1)^2, r0*(r0 - r1)/(r0 + 3*r1), (r0^2 - 3*r0*r1)/(r0 + r1))

術および中村の現代解法は乙円の半径についてかなり複雑な式を提示している。
SymPy での結果は,乙円の半径は,r0*(r0 - r1)/(r0 + 3*r1) すなわち,「大円の半径と大円の半径と甲円の半径の差をかけたものを大円の半径と甲円の半径の3倍の和で割る」となり,藤井貞夫が「長野県和算研究会報No.5」で指摘したものになった。

大円の直径が 1 尺,甲円の直径が 5 寸のとき乙円の直径は 2 寸である。

2res[1][2](r0 => 10/2, r1 => 5/2).evalf() |> float |> println

   2.0

   r2 = 2.22222;  r3 = 1;  y2 = -1.66667;  乙円の直径 = 2

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 10//2
   r1 = 5//2
   (r2, r3, y2) = (4*r0*r1*(r0 - r1)/(r0 + r1)^2, r0*(r0 - r1)/(r0 + 3*r1), (r0^2 - 3*r0*r1)/(r0 + r1))
   @printf("r2 = %g;  r3 = %g;  y2 = %g;  乙円の直径 = %g\n", r2, r3, y2, 2r3)
   plot()
   circle(0, 0, r0, :blue)
   circle(0, r0 - r1, r1, :red)
   circle(r2, y2, r2, :green)
   circle(-r2, y2, r2, :green)
   circle(0, r3 - r0, r3, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r0 - r1, " 甲円:r1,(0,r0-r1)", :red, :left, :vcenter)
       point(r2, y2, "(等円:r2,(r2,y2)", :green, :center, delta=-delta)
       point(0, r3 - r0, " 乙円:r3\n (0,r3-r0)", :black, :left, :vcenter)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 


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

コメントを投稿

Julia」カテゴリの最新記事