裏 RjpWiki

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

算額(その1017)

2024年05月31日 | Julia

算額(その1017)

四 岩手県花巻市南笹間 東光寺 慶応2年(1866)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html

外円の中に交差する大円 2 個,中円 2 個,小円 6 個を入れる。外円の直径が与えられたとき,小円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, R - r1)
中円の半径と中心座標を r2, (0, R - r2)
小円の半径と中心座標を r3, (x3, y3), (R - r3, 0)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

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

   1-element Vector{NTuple{5, Sym{PyCall.PyObject}}}:
    (2*R/3, R/3, R/5, 4*sqrt(3)*R/15, 2*R/5)

外円の半径を R とすれば,大円,中円,小円の半径はその 2/3, 1/3, 1/5 である。
例えば,外円の直径が 15 寸のとき,大円,中円,小円の直径は 10 寸,5 寸,3 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 15/2
   (r1, r2, r3, x3, y3) = (2*R/3, R/3, R/5, 4*sqrt(3)*R/15, 2*R/5)
   @printf("R = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  y3 = %g\n", R, r1, r2, r3, x3, y3)
   plot()
   circle(0, 0, R)
   circle(0,0, r2, :gray90)
   circle22(0, R - r1, r1, :magenta)
   circle22(0, R - r2, r2, :blue)
   circle4(x3, y3, r3, :green)
   circle2(R - r3, 0, r3, :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 - r1, "大円:r1\n(0,R-r1)", :red, :center, delta=-delta/2)
       point(0, R - r2, "中円:r2\n(0,R-r2)", :blue, :center, delta=-delta/2)
       point(x3, y3, "小円:r3\n(x3,y3)", :green, :center, delta=-delta/2)
       point(R-r3, 0, "小円:r3\n(R-r3,0)", :green, :center, delta=-delta/2)
       point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(0, R, " R", :red, :left, :bottom, delta=delta/2)
   end
end;

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

算額(その1016)

2024年05月31日 | Julia

算額(その1016)

四 岩手県花巻市南笹間 東光寺 慶応2年(1866)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html

外円の中に水平な弦を 2 本引き,区画された領域に,上円,中円,下円を入れる。上円,下円の直径が与えられたときに,外円の直径を求める術を述べよ。

外円の半径と中心座標を R, (0, 0)
上円の半径と中心座標を r1, (0, R - r1)
中円の半径と中心座標を r2, (2r2, R - 2r1 - r2), (r2, 2r3 - R + r2), (3r2, 2r3 - R + r2)
下円の半径と中心座標を r3, (0, r2 - R)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive, r3::positive
eq1 = 9r2^2 + (2r3 - R + r2)^2 - (R - r2)^2
eq2 = r2^2 + ((R - 2r1 - r2) - (2r3 - R + r2))^2 - 4r2^2
eq3 = 2r1 + 2r3 + √Sym(3)*r2 + 2r3 - 2R
res = solve([eq1, eq2], (R, r2))[2];  # 2 0f 2

# R
res[1] |> simplify |> println

   r1 + sqrt(r3)*sqrt(12*r1 + r3)/6 + 7*r3/6 + sqrt(36*r1*r3 + 6*r3^(3/2)*sqrt(12*r1 + r3) + 6*r3^2)/9

# r2
res[2] |> simplify |> println

   sqrt(36*r1*r3 + 6*r3^(3/2)*sqrt(12*r1 + r3) + 6*r3^2)/9

r1, r3 が与えられたとき,R(r2)は以下のように計算できる。
   t = sqrt(12r1*r3 + r3^2)
   r2 = sqrt(6r3*(6r1 + t + r3))/9
   R = r1 + r2 + (t + 7r3)/6

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r3) = (3, 5)./2
   t = sqrt(12r1*r3 + r3^2)
   r2 = sqrt(6r3*(6r1 + t + r3))/9
   R = r1 + r2 + (t + 7r3)/6
   plot()
   circle(0, 0, R)
   circle(0, R - r1, r1, :magenta)
   circle2(2r2, R - 2r1 - r2, r2, :blue)
   circle2(r2, 2r3 - R + r2, r2, :blue)
   circle2(3r2, 2r3 - R + r2, r2, :blue)
   circle(0, r3 - R, r3, :green)
   y1 = R - 2r1
   x1 = sqrt(R^2 - y1^2)
   segment(-x1, y1, x1, y1, :orange)
   y2 = 2r3 - R
   x2 = -sqrt(R^2 - y2^2)
   segment(-x2, y2, x2, y2, :orange)
   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 - r1, "上円:r1\n(0,R-r1)", :magenta, :center, :bottom, delta=delta/2)
       point(2r2, R - 2r1 - r2, "中円:r2\n(2r2,R-2r1-r2)", :blue, :center, :bottom, delta=delta/2)
       point(r2, 2r3 - R + r2, "中円:r2\n(r2,2r3-R+r2)", :blue, :center, :bottom, delta=delta/2)
       point(3r2, 2r3 - R + r2, "中円:r2\n(3r2,2r3-R+r2)", :blue, :center, :bottom, delta=delta/2)
       point(0, r3 - R, "下円:r3\n(0,r3-R)", :green, :center, :bottom, delta=delta/2)
       point(0, R, " R", :red, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(0, R - 2r1, " R-2r1", :orange, :center, delta=-delta/2)
       point(0, 2r3 - R, " 2r3-R", :orange, :center, :bottom, delta=delta/2)
   end
end;

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

算額(その1015)

2024年05月31日 | Julia

算額(その1015)

三 岩手県花巻市太田 音羽山清水観世音堂 明治25年(1892)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

半円と斜線で区切られた領域に大円,中円,小円を入れる。小円の直径が与えられたとき,中円の直径を求める術を述べよ。

半円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, r1)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (0, -r3)
小円の中心と中円の中心を結ぶ線分と x 軸との交点座標を (x0, 0)
とする。

⊿OCE∽⊿DBE∽OEAより,r3 が与えられたときに,まず r1 を求め,最終的に r2 を求める。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive,
     r2::positive, x2::positive,
     r3::positive, x3::positive, y3::positive,
     x0::positive
R = 2r1
eq1 = dist2(0, -r3, x3, y3, 0, r1, r1 - r3)
eq3 = x2^2 + r2^2 - (R - r2)^2 |> expand
eq4 = x3^2 + y3^2 - (R - r3)^2 |> expand
eq5 = x3^2 + (y3 - r1)^2 - (r1 + r3)^2 |> expand;

1. 小円と中円の関係

eq6 = r3*(x2 - x0)/x0 - r2  # = r3/x0 - r2/(x2 - x0)
eq6 |> println

   -r2 + r3*(-x0 + x2)/x0

2. 小円と大円の関係

eq7 = r3/x0 - x0/r1;
eq7 |> println

   r3/x0 - x0/r1

eq7 を解いて x0 を求める。

ans_x0 = solve(eq7, x0)[1]
ans_x0 |> println

   sqrt(r1)*sqrt(r3)

eq3 を解いて x2 を求める。

ans_x2 = solve(eq3, x2)[1];
ans_x2 |> println

   2*sqrt(r1)*sqrt(r1 - r2)

3. 中円の半径を求める

eq = eq6(x0 => ans_x0, x2 => ans_x2)
eq |> println

   -r2 + sqrt(r3)*(-sqrt(r1)*sqrt(r3) + 2*sqrt(r1)*sqrt(r1 - r2))/sqrt(r1)

ans_r2 = solve(eq, r2)[1]
ans_r2 |> println

   2*sqrt(r3)*sqrt(r1 + 2*r3) - 3*r3

eq4, eq5 から x3, y3 を求め,eq1 に代入し,eq1 を解いて r1 を求め, 上の式に代入する。

(ans_x3, ans_y3) = solve([eq4, eq5], (x3, y3))[1];
ans_x3 |> println
ans_y3 |> println

   2*sqrt(2)*sqrt(r3)*sqrt(r1 - r3)
   2*r1 - 3*r3

eq1 = eq1(x3 => ans_x3, y3 => ans_y3) |> simplify
eq1 |> println

   -4*r1^4 + 16*r1^3*r3 + 8*r1^2*r3^2 - 16*r1*r3^3 - 4*r3^4

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

   r3*(2 + sqrt(5))

r2 = ans_r2(r1 => ans_r1) |> simplify
r2 |> println

   r3*(-3 + 2*sqrt(sqrt(5) + 4))

中円の半径 r2 は,小円の半径 r3 の (2sqrt(√5 + 4) - 3) 倍である。
小円の直径が 1 寸のとき 中円の直径は 1.99442408191367 寸である。

2r2(r3 => 0.5).evalf() |> println

   1.99442408191367

術では,sqrt(sqrt(80) + 16) - 3 倍となっている。

(sqrt(sqrt(Sym(80)) + 16) - 3) |> println

   -3 + sqrt(4*sqrt(5) + 16)

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

r3 = 0.5;  r1 = 2.11803;  R = 4.23607;  r2 = 0.997212;  x2 = 3.08152;  x3 = 2.54404;  y3 = 2.73607;  x0 = 1.02909
y00 = 1.30902;  x01 = 3.20585;  y01 = 2.76889

#=
r1 = r3*(2 + √5)
R = 2r1
x3 = 2sqrt(2r3*(r1 - r3))
y3 = 2r1 - 3r3
r2 = r3*(2sqrt(√5+ 4) - 3)
x2 = 2*sqrt(r1)*sqrt(r1 - r2)
x0 = sqrt(r1)*sqrt(r3)
=#

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 0.5
   r1 = r3*(2 + √5)
   R = 2r1
   x3 = 2sqrt(2r3*(r1 - r3))
   y3 = 2r1 - 3r3
   r2 = r3*(2sqrt(√5+ 4) - 3)
   x2 = 2sqrt(r1*(r1 - r2))
   x0 = sqrt(r1*r3)
   tanθ = (r3 + y3)/x3
   y00 = x0*tanθ
   x01 = (tanθ*y00 + sqrt(4*r1^2*tanθ^2 + 4*r1^2 - y00^2))/(tanθ^2 + 1)
   y01 = (tanθ*sqrt(4*r1^2*tanθ^2 + 4*r1^2 - y00^2) - y00)/(tanθ^2 + 1)
   @printf("小円の直径が %g のとき,中円の直径は %g である。\n", 2r3, 2r2)
   @printf("r3 = %g;  r1 = %g;  R = %g;  r2 = %g;  x2 = %g;  x3 = %g;  y3 = %g;  x0 = %g\n", r3, r1, R, r2, x2, x3, y3, x0)
   @printf("y00 = %g;  x01 = %g;  y01 = %g\n", y00, x01, y01)
   plot()
   circle(0, 0, R, :orange, beginangle=0, endangle=180)
   segment(-R, 0, R, 0, :orange)
   circle(0, r1, r1)
   circle2(x2, r2, r2, :green)
   circle2(x3, y3, r3, :blue)
   circle(0, -r3, r3, :blue)
   segment(x01, y01, 0, -y00)
   segment(-x01, y01, 0, -y00)
   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=:gray80, lw=0.5)
       point(x3, y3, "小円:r3,(x3,y3) ", :blue, :right, :vcenter)
       point(0, -r3, " C  小円:r3,(0,-r3)", :blue, :left, delta=-delta)
       point(x0, 0, "E (x0,0)", :black, :left, delta=-delta)
       point(x2, 0, "D (x2,0)", :black, :left, delta=-delta)
       point(0, r1, "大円:r1,(0,r1)\nA", :red, :right, :bottom, delta=delta)
       point(x2, r2, "中円:r2\n(x2,r2)\nB", :green, :center, :bottom, delta=delta)
       point(0, 0, "O ", :black, :right, :bottom, delta=delta)
       point(0, -y00, " -y00", :black, :left, :vcenter)
       point(x01, y01, " (x01,y01)", :black, :left, :vcenter)
       segment(x2, r2, 0, -r3, :magenta)
       segment(0, r1, 0, -r3, :magenta)
       segment(0, r1, x0, 0, :magenta)
       segment(x2, r2, x2, 0, :magenta)
   end
end;

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

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

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