裏 RjpWiki

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

算額(その2018)

2024年08月14日 | Julia

算額(その2018)

(17) 兵庫県姫路市飾磨区英賀宮町 英賀神社 明治12年(1879)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.
キーワード:円1個,正方形5個

外円の中に合同な正方形を 5 個容れる。外円の直径が 10 寸のとき,正方形の一辺の長さはいかほどか。

正方形の一辺の長さを a
外円の半径と中心座標を R, (0, 0)
とおき,以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, R::positive
eq1 = √Sym(2)a + a/2 - R
res = solve(eq1, a)[1] |> simplify
res |> println
res(R => 10/2).evalf() |> println

   2*R*(-1 + 2*sqrt(2))/7
   2.61203874963741

正方形の一辺の長さ a は,外円の半径 R の 2(2√2 - 1)/7 倍である。
外円の直径が 10 寸のとき,正方形の一辺の長さは 5*2(2√2 - 1)/7 = 2.6120387496374144 寸である。

function draw_diamond(R, a, θ)
   (y0, x0) = (R - a/√2).*sincos(pi*θ/180)
   θ2 = 0:90:360
   x = x0 .+ a/√2*cosd.(θ2)
   y = y0 .+ a/√2*sind.(θ2)
   plot!(x, y, color=:blue, lw=0.5)
end
       
function diamond(R, a)
   for θ = 0:90:270
       draw_diamond(R, a, θ)
   end
end
       
function draw(R, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = R*2(2√2 - 1)/7
   @printf("外円の直径が %g のとき,正方形の一辺の長さは %g である。\n", 2R, a)
   plot()
   circle(0, 0, R)
   rect(-a/2, -a/2, a/2, a/2, :blue)
   diamond(R, a)
   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", :red, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(a/2, a/2, "(a/2,a/2)", :blue, :center, :bottom, delta=delta/2)
       point(R - a/√2, a/√2, "(R-a/√2,a/√2)", :blue, :center, :bottom, delta=delta/2)
   end
end;

draw(10/2, true)

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

算額(その2017)

2024年08月14日 | Julia

算額(その2017)

(17) 兵庫県姫路市飾磨区英賀宮町 英賀神社 明治12年(1879)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.
キーワード:円4個,正三角形

正三角形の中に,等円 4 個を容れる。正三角形の一辺の長さが 10 寸のとき,等円の直径はいかほどか。

正三角形の一辺の長さを 2a
等円の半径と中心座標を r, (0, √3a/3), (0, √3a/3 + 2r), (a - √3r, r)
とおき,以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, r::positive
eq1 = √Sym(3)a - (√Sym(3)a/3 + 2r) - 2r
res = solve(eq1, r)[1]
res |> println

   sqrt(3)*a/6

等円の直径は,正三角形の一辺の長さの √3/6 倍である。
正三角形の一辺の長さが 10 寸のとき,等円の直径は 2.8867513459481287 寸である。

function draw(a, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = √3a/6
   @printf("三角形の一辺の長さが %g のとき,等円の直径は %g である。\n", 2a, 2r)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:blue, lw=0.5)
   circle(0, √3a/3, r)
   circle(0, √3a/3 + 2r, r)
   circle2(a - √3r, r, r)
   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, √3a/3, "(0,√3a/3)", :red, :center, delta=-delta/2)
       point(0, √3a/3+2r, "(0,√3a/3+2r)", :red, :center, delta=-delta/2)
       point(a - √3r, r, "(a-√3r,r)", :red, :center, delta=-delta/2)
       point(0, √3a, "√3a", :blue, :center, :bottom, delta=delta/2)
       point(a, 0, " a", :blue, :center, :bottom, delta=delta/2)
   end
end;

draw(10/2, true)

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

算額(その2016)

2024年08月14日 | Julia

算額(その2016)

(19) 兵庫県姫路市広峯山 広峰神社 明治18年(1885)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.
キーワード:円4個,楕円

楕円の中に,甲円 4 個,乙円 2 個を容れる。楕円の長径と短径が与えられたとき,甲円の直径はいかほどか。

注:甲円,乙円,楕円は共通接点で接する。
楕円に内接する 2 個の等円については算法助術の公式84を適用する。

楕円の長半径,短半径,中心座標を a, b, (0, 0)
甲円の半径と中心座標を r1, (r1, 0)
乙円の半径と中心座標を r2, (x2, r2)
甲円,乙円,楕円の共通接点を (x0, y0)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r1::positive,
     r2::positive, x2::positive,
     x0::positive, y0::positive
eq1 = (2r1)^2/4 - (a^2 - b^2)*(b^2 - r1^2)/b^2  # 算法助術の公式84
eq2 = (x2 - r1)^2 + r2^2 - (r1 - r2)^2
eq3 = (x0 - r1)^2 + y0^2 - r1^2
eq4 = (x0 - x2)^2 + (y0 - r2)^2 - r2^2
eq5 = x0^2/a^2 + y0^2/b^2 - 1;

連立方程式を解くと,かなり長い式で表される解が得られる。

res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, x2, x0, y0))[4]  # 4 0f 4

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

乙円の半径は r1 = b*(a^2 - b^2)*sqrt(1/(a - b))/(a*sqrt(a + b)) となるが,eq1 を r1 についてとくと,もう少し簡約化された解 r1 = b*sqrt(a^2 - b^2)/a が得られる。

ans_r1 = solve(eq1, r1)[1]
ans_r1 |> println

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

甲円の式も長いものであるが,たとえば長径,短径がそれぞれ 11170 寸,4204 寸のとき,式に代入することにより甲円の直径は 1934.000000079764 寸である。

res[2] |> println

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

2res[2](a => 11170/2, b => 4204/2).evalf() |> println

   1934.00000007983

その他のパラメータは以下のとおりである

  a = 5585;  b = 2102;  r1 = 1947.44280674992, r2 = 967.000000039889, x2 = 2109.24;  x0 = 2268.82;  y0 = 1920.74

function draw(a, b, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = b*sqrt(a^2 - b^2)/a
   println("r1 = $r1")
   (r1, r2, x2, x0, y0) = (b*(a^2 - b^2)*sqrt(1/(a - b))/(a*sqrt(a + b)), (a^4*sqrt((a^2 - 2*b^2)/(a^2 - b^2)) - a^4*sqrt(a + b)*sqrt(1/(a - b)) + a^3*b*sqrt(a + b)*sqrt(1/(a - b)) - 2*a^2*b^2*sqrt((a^2 - 2*b^2)/(a^2 - b^2)) + 2*a^2*b^2*sqrt(a + b)*sqrt(1/(a - b)) - 2*a*b^3*sqrt(a + b)*sqrt(1/(a - b)) + b^4*sqrt((a^2 - 2*b^2)/(a^2 - b^2)))/b^3, (-a^2*sqrt((a^2 - 2*b^2)/(a^2 - b^2)) + a^2*sqrt(a + b)*sqrt(1/(a - b)) - a*b*sqrt(a + b)*sqrt(1/(a - b)) + b^2*sqrt((a^2 - 2*b^2)/(a^2 - b^2)))/b, a*b*sqrt(1/(a - b))/sqrt(a + b), b*sqrt((a^2 - 2*b^2)/(a^2 - b^2)))
   println("r1 = $r1")
   @printf("長径,短径がそれぞれ %g,%g のとき,甲円の直径は %.15g 寸である。\n", 2a, 2b, 2r2)
   @printf("a = %g;  b = %g;  r1 = %.15g, r2 = %.15g, x2 = %g;  x0 = %g;  y0 = %g\n", a, b, r1, r2, x2, x0, y0)
   plot()
   ellipse(0, 0, a, b, color=:green)
   circle2(r1, 0, r1)
   circle4(x2, r2, r2, :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(x2, r2, "甲円:r2\n(x2,r2)", :blue, :center, delta=-delta)
       point(r1, 0, " 乙円:r1\n (r1,0)", :red, :center, delta=-4delta)
       point(x0, y0, "(x0,y0)", :green, :left, :bottom, delta=delta)
       point(a, 0, " a", :green, :left, :bottom, delta=delta)
       point(0, b, " b", :green, :left, :bottom, delta=delta)
   end
end;

draw(11170/2, 4204/2, true)

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

算額(その2015

2024年08月14日 | Julia

算額(その2015)

(20) 兵庫県青垣町遠坂字後岶 熊野神社 明治18年(1885)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.
キーワード:円5個,長方形

長方形の中に等円 4 個,都円 1 個を容れる。長方形の長辺は等円の直径の 3 倍である(短辺は図より等円の 2 倍である)。等円の直径が 10 寸のとき,都円の直径はいかほどか。

都円の半径と中心座標を r1, (0, 0)
等円の半径と中心座標を r2, (a - r2, r2)
長方形の長辺と短辺を 2a, 2b とおく。a = 3r2, b = 2r2
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r1::positive, r2::positive
a = 3r2
b = 2r2
eq1 = (a - r2)^2 + r2^2 - (r1 + r2)^2
res = solve(eq1, r1)[1]
res |> println

   r2*(-1 + sqrt(5))

都円の半径 r1 は,等円の半径 r2 の (√5 - 1) 倍である。
等円の直径が 10 寸のとき,都円の直径は 10(√5 - 1) = 12.360679774997898 寸である。

算額の「答」では 12.227362 寸としているが,「術」を見る限り,正確な値ではない。実際に図を描いても,「答」では等円と都円が外接しない。

function draw(r2, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 3r2
   b = 2r2
   r1 = (√5 - 1)r2
   @printf("r2 = %g;  a = %g;  b = %g;  r1 = %g\n", r2, a, b, r1)
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:green, lw=0.5)
   circle(0, 0, r1)
   circle4(a - r2, r2, r2, :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(0, 0, "都円:r1,(0,0)", :red, :center, delta=-delta)
       point(a - r2, r2, "等円:r,(a-r2,r2)", :blue, :center, delta=-delta)
       point(a, b, "(a,b)", :green, :right, :bottom, delta=delta)
   end
end;

draw(10/2, true)

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

算額(その2014)

2024年08月14日 | Julia

算額(その2014)

(20) 兵庫県青垣町遠坂字後岶 熊野神社 明治18年(1885)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.
キーワード:円2個,台形,水平線

台形の中に大円と水平線を描き,できる円弧の中に月円を容れる。大頭(下底),小頭(上底)がそれぞれ 3 寸,2 寸のとき,月円の直径はいかほどか。

大頭,小頭を b, a
大円の半径と中心座標を R, (0, R)
月円の半径と中心座標を r, (0, a + r)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, R::positive, r::positive
eq1 = dist2(R, a, -R, b, 0, R, R)
eq2 = (2R - a) - 2r
res = solve((eq1, eq2), (r, R))[1]

   (-a*(a - b)/(2*(a + b)), a*b/(a + b))

月円の半径は a*(b - a)/(2a + 2b) である。

大頭(下底),小頭(上底)がそれぞれ 3 寸,2 寸のとき,月円の直径は 0.4 寸(4 分)である。

ちなみに,大円の直径は 2.4 寸である。

a = 2
b = 3
a*(b - a)/(2a + 2b) |> println
a*b/(a + b) |> println

   0.2
   1.2

function draw(a, b, more=false)
   pyplot(size=(300, 300), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, R) = (-a*(a - b)/(2*(a + b)), a*b/(a + b))
   @printf("a = %g;  b = %g;  r = %g;  R = %g\n", a, b, r, R)
   plot([R, R, -R, -R, R], [0, a, b, 0, 0], color=:blue, lw=0.5)
   circle(0, R, R)
   circle(0, a + r, r, :green)
   segment(-R, a, R, a, :magenta)
   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,(0,R)", :red, :center, delta=-delta/2)
       point(0, a + r, "月円:r,(0,a+r)", :green, :center, delta=-r)
       point(R, a, " (R,a)", :blue, :left, :vcenter)
       point(-R, b, " (-R,b)", :blue, :left, :bottom)
   end
end;

draw(2, 3, true)

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

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

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