裏 RjpWiki

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

算額(その510)

2023年11月25日 | Julia

算額(その510)

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション

https://dl.ndl.go.jp/pid/3508431/1/11

(16) 京都府京都市右京区山ノ内宮脇町 山王神社 明治23年(1890)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.

正方形内に甲円 1 個,乙円 2 個が入っている。甲円の直径が 9 寸のとき,正方形の一辺の長さはいかほどか。

正方形の一辺の長さを a = 4r2
甲円の半径と中心座標を r1, (0, 4r2 - r1)
乙円の半径と中心座標を r2, (r2, r2)
として方程式を解く。

include("julia-source.txt");

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

   8*r1/9

乙円の半径は甲円の半径の 8/9 である。
甲円の直径が 9 寸ならば乙円の直径は 8 寸,正方形の一辺の長さはその 2 倍,16 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 9/2
   r2 = 8r1/9
   a = 4r2
   @printf("r1 = %g;  r2 = %g;  正方形の一辺の長さ = %g\n", r1, r2, a)
   plot([a/2, a/2, -a/2, -a/2, a/2], [0, a, a, 0, 0], color=:blue, lw=0.5)
   circle(r2, r2, r2, :green)
   circle(-r2, r2, r2, :green)
   circle(0, a - r1, r1)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, a - r1, " 甲円:r1,(0,a-r1)", :red, :left, :vcenter)
       point(r2, r2, " 乙円:r2(r2,r2)", :green, :left, :vcenter)
       point(a/2, 0, "a/2=2r2 ", :blue, :right, :bottom, delta=delta)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その509)

2023年11月25日 | Julia

算額(その509)

岩手県一関市上大槻街路 一関神明社 天保5年(1834)
http://www.wasan.jp/iwate/sinmeisha.html

正方形内に,四分円 2 個,半円 1 個,天円,地円が入っている。
地円の直径が 1 寸のとき,天円の直径はいかほどか。

正方形の一辺の長さを a
天円の半径と中心座標を r1, (r1, a - r1)
地円の半径と中心座標を r2, (x2, r2)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive, x2::positive

eq1 = x2^2 + (a - r2)^2 - (a + r2)^2
eq2 = 2(a - r1)^2 - (a + r1)^2
eq3 = (a - x2)^2 + (a//2 -r2)^2 - (a//2 + r2)^2
res = solve([eq1, eq2, eq3], (r1, r2, x2))

   4-element Vector{Tuple{Sym, Sym, Sym}}:
    (a*(3 - 2*sqrt(2)), -a/2 + a*(2 - sqrt(2)), a*(2 - sqrt(2)))
    (a*(3 - 2*sqrt(2)), -a/2 + a*(sqrt(2) + 2), a*(sqrt(2) + 2))
    (a*(2*sqrt(2) + 3), -a/2 + a*(2 - sqrt(2)), a*(2 - sqrt(2)))
    (a*(2*sqrt(2) + 3), -a/2 + a*(sqrt(2) + 2), a*(sqrt(2) + 2))

4 組の解が得られるが,1 番目のものが適解である。
r1 は r2 の 2 倍である。
地円の直径が 1 寸のとき,天円の直径は 2 寸である。

res[1][1]/res[1][2] |> simplify |> println

   2

ちなみに,地円の直径が 1 寸となるのは,正方形の一辺の長さが 4*sqrt(2) + 6 のときである。

   a = 11.6569;  r1 = 2;  r2 = 1;  x2 = 6.82843

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 4*sqrt(2) + 6
   (r1, r2, x2) = (a*(3 - 2*sqrt(2)), -a/2 + a*(2 - sqrt(2)), a*(2 - sqrt(2)))
   @printf("a = %g;  r1 = %g;  r2 = %g;  x2 = %g\n", a, r1, r2, x2)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:blue, lw=0.5)
   circle(0, a, a, :green, beginangle=270, endangle=360)
   circle(a, 0, a, :green, beginangle=90, endangle=180)
   circle(a, a/2, a/2, :magenta, beginangle=90, endangle=270)
   circle(r1, a - r1, r1)
   circle(x2, r2, r2, :orange)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(0, a, " a", :blue, :left, :bottom, delta=delta/2)
       point(a, a/2, "(a,a/2) ", :magenta, :right, :vcenter)
       point(r1, a - r1, " 天円:r1\n (r1,a-r1)", :red, :left, :vcenter)
       point(x2, r2, "地円:r2,(x2,r2)", :black, :center, :top, delta=-delta/2)
   else
       plot!(showaxis=false)
   end
end;

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

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

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