裏 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でシェアする

算額(その1014)

2024年05月30日 | Julia

算額(その1014)

一〇四 桶川町加納 氷川天満神社 明治43年(1910)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

一辺の長さが 15 寸の正三角形に,円 2 個を入れる。円の直径と正三角形の高さはいかほどか。

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

include("julia-source.txt");

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

   a*(-1 + sqrt(3))/4
   2.74519052838329
   12.9903810567666

円の半径は正三角形の一辺の長さの (√3 - 1)/4 倍である。
正三角形の一辺の長さが 15 寸のとき,円の直径は 15*(√3 - 1)/2 = 5.490381056766579 寸である。
また正三角形の高さは一辺の長さの √3/2 倍である。すなわち 15*√3/2 = 12.990381056766578 寸である。

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

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

算額(その1013)

2024年05月30日 | Julia

算額(その1013)

七八 加須市外野 棘脱地蔵堂 明治9年(1876)

埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

正方形内に甲円 1 個と乙円 3 個が入っている。乙円の直径が 10 寸のとき,甲円の直径はいかほどか。

正方形の一辺の長さを 2a; 2a = 6r2
甲円の半径と中心座標を r1, (0, 2a - r1)
乙円の半径と中心座標を r2, (2r2, r2)
とおき,以下の方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive
a = 3r2
eq1 = 4r2^2 + (2a - r1 - r2)^2 - (r1 + r2)^2
res = solve(eq1, r1)[1]
res |> println
res(r2 => 10/2).evalf() |> println

   7*r2/3
   11.6666666666667

甲円の半径は乙円の半径の 7/3 倍である。
乙円の直径が 10 寸のとき,甲円の直径は 10*7/3 = 23.333333333333332 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 10/2
   r1 = 7r2/3
   a = 3r2
   @printf("乙円の直径が %g のとき,甲円の直径は %g\n", 2r2, 2r1)
   plot([a, a, -a, -a, a], [0, 2a, 2a, 0, 0], color=:blue, lw=0.5)
   circle(0, 2a - r1, r1, :green)
   circle2(2r2, r2, r2)
   circle(0, r2, r2)
   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(0, 2a - r1, "甲円:r1,(0,2a-r1)", :green, :center, delta=-delta/2)
       point(0, r2, "乙円:r2\n(0,r2)", :red, :center, delta=-delta/2)
       point(2r2, r2, "乙円:r2\n(2r2,r2)", :red, :center, delta=-delta/2)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(a, 2a, "(a,2a)", :blue, :right, :bottom, delta=delta/2)
   end
end;

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

算額(その1012)

2024年05月29日 | Julia

算額(その1012)

八七 加須市大越六間 天神社 明治14年(1881)

埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

5 個の等円が隣同士外接し,かつ,内側にある正五角形の辺にも外接している。等円の直径が 49.0165 寸のとき,正五角形の一辺の長さを求めよ。

等円の中心は,原点 O を中心とする半径 OA の円周上にある。OA = r1 + r2*cos(36°)
等円の半径 AD = r1
AB/2 = r1 = OA*sin(36°)
正五角形の頂点は,原点 O を中心とする半径 OC の円周上にある。 OC = r2 = 4r1*(-sqrt(5 - √5) + 2√2)/((1 + √5)*sqrt(5 - √5))
∠COD は 36°
正五角形の辺と等円の接点は,原点 O を中心とする半径 OD の円周上にある。
正五角形の一辺の長さは 2CD = 2OC*sin(36°)

等円の直径 2r1 = 49.0165 のとき,五角形の一辺の長さは 24.97515419514368 である。

術では五角形の一辺の長さは 25 寸とあるが,正確ではない。
正確に25 寸であるためには等円の直径が 49.0652626376286 でなければならない。
問で「等円径四十九寸令壱リ六毛五糸」とあるのが誤記で,「等円径四十九寸令分六釐五毛」なのかもしれない。

以下により,r1, r2 を求め,正五角形の一辺の長さ 2r2*sin(36°) = 24.97515419514368 を求める。

include("julia-source.txt");

using SymPy

@syms r1::poitive, r2::poitive

eq1 = (r1 + r2*cosd(Sym(36)))*sind(Sym(36)) - r1;
solve(eq1, r2)[1] |> println

   4*r1*(-sqrt(5 - sqrt(5)) + 2*sqrt(2))/((1 + sqrt(5))*sqrt(5 - sqrt(5)))

r1 = 49.0165/2
r2 = 4*r1*(-sqrt(5 - sqrt(5)) + 2*sqrt(2))/((1 + sqrt(5))*sqrt(5 - sqrt(5)))
2r2*sind(36)

   24.97515419514368

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 49.0165/2
   r2 = 4r1*(-sqrt(5 - √5) + 2√2)/((1 + √5)*sqrt(5 - √5))
   l = 2r2*sind(36)
   println("i = $l")
   θ = 54:72:414
   x = cosd.(θ)
   y = sind.(θ)
   plot()
   r0 = r1 + r2*cosd(36)
   for i = 1:5
       segment(r2*x[i], r2*y[i], r2*x[i + 1], r2*y[i + 1], :blue, lw=0.5)
       circle(r0*cosd(θ[i] + 36), r0*sind(θ[i] + 36), r1)
   end
   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)
       circle(0, 0, r0, :gray70)
       #circle(r0*cosd(θ[1] + 36), r0*sind(θ[1] + 36), r1, :magenta)
       point(r0*cosd(θ[1] + 36), r0*sind(θ[1] + 36), " A", :green, :left, :bottom, delta=delta/2)
       segment(0, 0, r0*cosd(θ[1] + 36), r0*sind(θ[1] + 36), :green)
       point(r0*cosd(θ[5] + 36), r0*sind(θ[5] + 36), " B")
       segment(r0*cosd(θ[1] + 36), r0*sind(θ[1] + 36), r0*cosd(θ[5] + 36), r0*sind(θ[5] + 36), :green)
       point(0, 0, " O", :green, :left, delta=-delta/2)
       circle(0, 0, r2, :gray80)
       point(r2*cosd(θ[1]), r2*sind(θ[1]), "C", :green, :left, :bottom, delta=0.1delta, deltax=0.3delta)
       segment(0, 0, r2*cosd(θ[1]), r2*sind(θ[1]), :green)
       point(0, r2*sind(θ[1]), " D", :green, :left, delta=-delta/2)
       circle(0, 0, r2*sind(θ[1]), :gray80)
   end
end;

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

算額(その1011)

2024年05月29日 | Julia

算額(その1011)

八六 加須市多聞寺 愛宕神社 明治13年(1880)

埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

長方形の中に甲円 4 個,乙円 2 個,丙円 8 個と菱形が入っている。甲円の直径が 1 寸のとき,乙円の直径はいかほどか。

長方形の長辺と短辺の長さを 2a, 2b; b = 2r1
菱形の短い方の対角線の長さを 2c; 長い方の対角線は 2b
甲円の半径と中心座標を r1, (a - r1, r1)
乙円の半径と中心座標を r2, (0, y2)
丙円の半径と中心座標を r3, (x3, b - r3), (r3, r3)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::poitive, b::poitive, c::poitive,
     r1::poitive, r2::poitive, y2::poitive,
     r3::poitive, x3::poitive
b = 2r1
eq1 = (a - r1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2 |> expand
eq2 = r3^2 + (y2 - r3)^2 - (r2 + r3)^2 |>expand
eq3 = dist2(0, b, c, 0, x3, b - r3, r3)/4r1
eq4 = dist2(0, b, c, 0, a - r1, r1, r1)/(4r1^2)
eq5 = dist2(0, b, c, 0, 0, y2, r2)
eq6 = dist2(0, b, c, 0, r3, r3, r3)/(4c*r1);

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (a, c, r2, y2, r3, x3) = u
   return [
       a^2 - 2*a*r1 - 2*a*x3 + r1^2 - 4*r1*r3 + 2*r1*x3 + x3^2,  # eq1
       -r2^2 - 2*r2*r3 + r3^2 - 2*r3*y2 + y2^2,  # eq2
       -c*r3*x3 - r1*r3^2 + r1*x3^2,  # eq3
       a^2 - a*c - 2*a*r1 + c*r1,  # eq4
       4*c^2*r1^2 - 4*c^2*r1*y2 - c^2*r2^2 + c^2*y2^2 - 4*r1^2*r2^2,  # eq5
       c*r1 - c*r3 - 2*r1*r3 + r3^2,  # eq6
   ]
end;
r1 = 1/2
iniv = BigFloat[1.17, 0.29, 0.17, 0.39, 0.13, 0.17]
res = nls(H, ini=iniv)

   ([1.1666666666666667, 0.2916666666666667, 0.1701388888888889, 0.3923611111111111, 0.125, 0.16666666666666666], true)

甲円の直径が 1 寸のとき,乙円の直径は 0.340278 寸(3 分 4 厘有奇)である。
注1:「答」の乙円径□□四釐有奇」の欠損文字は「三分」であろう。
注2:「術」「甲円の直径を 45 倍し 144 で割る」というのは数値解を求めたということであろうが,45/144=0.3125 は答えとも不一致である。17/50 = 0.34 は良い近似であろう。

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

   a = 1.16667;  b = 1;  c = 0.291667;  r2 = 0.170139;  y2 = 0.392361;  r3 = 0.125;  x3 = 0.166667

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, c, r2, y2, r3, x3) = res[1]
   b = 2r1
   @printf("甲円の直径が %g のとき,乙円の直径は %g である。\n", 2r1, 2r2)
   @printf("a = %g;  b = %g;  c = %g;  r2 = %g;  y2 = %g;  r3 = %g;  x3 = %g\n", a, b, c, r2, y2, r3, x3)
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:black, lw=0.5)
   plot!([0, c, 0, -c, 0], [-b, 0, b, 0, -b], color=:orange, lw=0.5)
   circle4(a - r1, b - r1, r1)
   circle22(0, y2, r2, :green)
   circle4(x3, b - r3, r3, :blue)
   circle4(r3, r3, r3, :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(a - r1, r1, "甲円:r1,(a-r1,r1)", :red, :center, :bottom, delta=delta/2)
       point(0, y2, " 乙円:r2,(0,y2)", :green, :left, :vcenter)
       point(r3, b - r3, " 丙円:r3,(x3,b-r3)", :blue, :left, :vcenter)
       point(r3, r3, " 丙円:r3,(r3,r3)", :blue, :left, :vcenter)
       point(a, b, "(a,b)", :black, :right, :bottom, delta=delta/2)
   end
end;

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

算額(その1010)

2024年05月28日 | Julia

算額(その1010)

七八 加須市大字外野 棘脱地蔵堂 明治9年(1876)

埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

直線の上に甲円が 3 個あり,互いに接している。その上に両端の 2 円に接するように乙円を描く。乙円と中央の甲円の交差した部分と甲円と直線の間の 3 箇所に丙円を描く。甲円の直径が 10 寸のとき,乙円の直径はいかほどか。

注:算額ではそれぞれの円に名前を付けていないのであろうか(少なくとも『埼玉の算額』の図では円の名前はない)。通常,大きい順に甲・乙・丙と付けていくが,この算額では一番大きいのは乙円である。また,小さい 3 個の円は図には描かれているが名前もついておらず「問」の文言中にも言及がない。以下ではこれらを補完したうえで,解を示す。

甲円の半径と中心座標を r1, (2r1, r1)
乙円の半径と中心座標を r2, (0, y2)
丙円の半径と中心座標を r3, (r1, r3)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::poitive, y2::poitive,
     r2::poitive, r3::poitive
eq1 = y2 - r2 + 2r3 - 2r1
eq2 = 4r1^2 + (y2 - r1)^2 - (r1 + r2)^2
eq3 = r1^2 + (r1 - r3)^2 - (r1 + r3)^2
res = solve([eq1, eq2, eq3], (r2, r3, y2))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (13*r1/4, r1/4, 19*r1/4)

乙円の半径 r2 は,甲円の半径 r1 の 13/4 倍である。
甲円の直径が 10 寸のとき,乙円の直径は 10*13/4 = 32.5 寸である。

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

   r1 = 5;  r2 = 16.25;  r3 = 1.25;  y2 = 23.75

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 10/2
   (r2, r3, y2) = (13*r1/4, r1/4, 19*r1/4)
   @printf("甲円の直径が %g のとき,乙円の直径は %g である。\n", 2r1, 2r2)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  y2 = %g\n", r1, r2, r3, y2)
   plot()
   circle(0, y2, r2)
   circle2(2r1, r1, r1, :blue)
   circle(0, r1, r1, :blue)
   circle2(r1, r3, r3, :green)
   circle(0, y2 - r2 + r3, 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, y2, "乙円:r2,(0,y2)", :red, :center, delta=-delta/2)
       point(2r1, r1, "甲円:r1\n(2r1,r1)", :blue, :center, delta=-delta/2)
       point(r1, r3, "丙円:r3,(r1,r3)", :green, :center, delta=-6delta/2)
       plot!(ylims=(-4delta, y2 + r2 + 3delta))
       hline!([0], color=:black, lw=0.5)
   end
end;

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

算額(その1009)

2024年05月28日 | Julia

算額(その1009)

七八 加須市大字外野 棘脱地蔵堂 明治9年(1876)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.
キーワード円5個,正方形,直線上

直線の上に載っている 2 個の大円が交差する隙間に中円と正方形を入れる。また,直線との隙間に小円を入れる。中円の直径が最大になるときの小円の直径を求めよ。


大円の半径と中心座標を r1, (x1, 0)
正方形の対角線の長さを 2a
中円の半径と中心座標を r2, (0, a + r2)
小円の半径と中心座標を r3, (0, r3 - r1)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::poitive,  r1::poitive, x1::poitive,
     r2::poitive, r3::poitive
r1 = 35//20
eq1 = x1 - r1 + a
eq2 = x1^2 + (a + r2)^2 - (r1 - r2)^2
eq3 = x1^2 + (r3 - r1)^2 - (r1 + r3)^2
res = solve([eq1, eq2, eq3], (r2, a, x1))[1]

   ((8*sqrt(7)*r3^(3/2) - 7*sqrt(7)*sqrt(r3) + 14*r3)/(2*(4*r3 - 7)), -sqrt(7)*sqrt(r3) + 7/4, sqrt(7)*sqrt(r3))

f = res[1];

using Plots
pyplot(size=(400, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(f, xlims=(0.1,0.21), xlabel="r3", ylabel="r2")



r3 = 1.5 前後のとき,r2 は最大値を取る。
導関数を求め,導関数 = 0 を解けばよい。

r3 = 21/8 - 7*sqrt(2)/4 = 0.150126265847084 のとき,r2 は最大値 0.300252531694167 になる。

res = solve(diff(f))[1]
res |> println
res.evalf() |> println

   21/8 - 7*sqrt(2)/4
   0.150126265847084

f(r3 => 21/8 - 7*sqrt(2)/4).evalf() |> println

   0.300252531694167

小円の直径が 2*0.150126265847084 = 0.300252531694168 のとき,中円の直径は最大値 2*0.300252531694167 = 0.600505063388334 になる。
中円の直径は小円の直径のちょうど 2 倍である。

2*0.150126265847084 |> println
2*0.300252531694167 |> println

   0.300252531694168
   0.600505063388334

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

   r3 = 0.150126;  r2 = 0.300253;  a = 0.724874;  x1 = 1.02513

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 21/8 - 7*sqrt(2)/4
   (r2, a, x1) = ((8*sqrt(7)*r3^(3/2) - 7*sqrt(7)*sqrt(r3) + 14*r3)/(2*(4*r3 - 7)), -sqrt(7)*sqrt(r3) + 7/4, sqrt(7)*sqrt(r3))
   @printf("r3 = %g;  r2 = %g;  a = %g;  x1 = %g\n", r3, r2, a, x1)
   @printf("中円の直径は,小円の直径の %g 倍である。\n", r2/r3)
   plot()
   circle2(x1, 0, r1)
   circle(0, a + r2, r2, :blue)
   circle(0, - a - r2, r2, :blue)
   plot!([0, a, 0, -a, 0], [-a, 0, a, 0, -a], color=:green, lw=0.5) 
   circle(0, r3 - r1, r3, :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(x1, 0, "大円:r1,(x1,r1)", :red, :left, delta=-delta/2)
       point(0, a + r2, "     中円:r2,(0,a+r2)", :blue, :left, :vcenter)
       point(0, r3 - r1, "小円:r3,(0,r3-r1)", :magenta, :center, delta=-3delta)
       point(a, 0, "a ", :green, :right, :vcenter)
       plot!(ylims=(-r1 - 7delta, r1 + 3delta))
   end
end;

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

算額(その1008)

2024年05月28日 | Julia

算額(その1008)

七八 加須市大字外野 棘脱地蔵堂 明治9年(1876)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

正三角形の中に正方形を 3 個入れる。正三角形の一辺の長さが 10 寸のとき,正方形の一辺の長さはいかほどか。

正三角形の一辺の長さを 2a
正方形の一辺の長さを 2b
とおき,以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::poitive,  b::poitive
eq1 = 2b/(a - b) - 1/sqrt(Sym(3))
res = solve(eq1, b)[1] |> simplify
res |> println
res(a => 10/2).evalf() |> println

   a*(-1 + 2*sqrt(3))/11
   1.12004618869898

正方形の一辺の長さは正三角形の一辺の長さの (2√3 - 1)/11 倍である。
正三角形の一辺の長さが 10 寸のとき,正方形の一辺の長さは 10(2√3 - 1)/11 = 2.2400923773979584 寸である。

「答」は「等方面五厘有奇」と,とんでもないことを書いている。
「術」では,「置十二開平方内一个減余乗面以除十一个見寸合問」と,「12 の平方根(2√3)から 1 引いて正三角形の一辺の長さを掛けて 11 で割る」と正解を書いているのに。

function rect2(b, sign=1)
   segment(sign*b, 2b, sign*(b + 2b*cosd(30)), 2b + 2b*sind(30), :red)
   segment(0, 2b + 2b*√3/2, sign*b, 2b, :red)
   segment(0, 2b + 2b*√3/2, sign*(2b*cosd(30)), 2b + 2b*√3/2 + 2b*sind(30), :red)
   segment(sign*(b + 2b*cosd(30)), 2b + 2b*sind(30), sign*(2b*cosd(30)), 2b + 2b*√3/2 + 2b*sind(30), :red)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 10/2
   b = a*(2sqrt(3) - 1)/11
   @printf("正三角形の一辺の長さが %g のとき,正方形の一辺の長さは %g である。\n", 2a, 2b)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:blue, lw=0.5)
   rect(-b, 0, b, 2b, :red)
   rect2(b)
   rect2(b, -1)
   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(a, 0, "a", :blue, :left, :bottom, delta=delta/2)
       point(0, √3a, " √3a", :blue, :left, :vcenter)
       point(0, 2b, " 2b", :red, :left, :bottom, delta=delta/2)
       point(b, 0, " b", :red, :left, :vcenter)
   end
end;

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

算額(その1007)

2024年05月28日 | Julia

算額(その1007)

七〇 加須市大字外野 棘脱地蔵堂 明治6年(1873)

埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

外円の中に直径が外円と同じ 1/3 円 4 個,甲円 2 個,乙円 4 個が入っている。乙円の直径が 1 寸のとき,甲円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
1/3 円の半径と中心座標を R, (-R\*cos(pi/3), R\*sin(pi/3))
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, ((R - r2)cos(pi/6), (R - r2)sin(pi/6))
とおき以下の方程式を解く。

二元連立方程式のように見えるが,それぞれが独立な方程式で,逆に二元連立方程式として解くように SymPy を使うと(有限の時間内には)解が求まらない。
eq1 から r2 を未知数とする R,引き続いて eq2 から R を未知数とする r1 が求まる。

include("julia-source.txt");

using SymPy

@syms x0::poitive,  y0::poitive,  x2::poitive, y2::positive,
     R::positive, r1::positive, r2::positive
(s30, s120) = (Sym(30), Sym(120))
(x0, y0) = R .* (cosd(s120), sind(s120))
(x2, y2) = (R - r2) .* (cosd(s30), sind(s30))
eq1 = (x2 - x0)^2 + (y2 - y0)^2 - (R + r2)^2
eq2 = x0^2 + (y0 - R + r1)^2 - (R - r1)^2;
# res = solve([eq1, eq2], (r1, R))

ans_R = solve(eq1, R)[1]
ans_R |> println

   4*r2

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

   R*(3 - sqrt(3))/3

甲円の直径は乙円の直径の R*(3 - sqrt(3))/3 = 4r2(3 - √3)/3 = 4(3 - √3)/3 = (12- √48)/3 倍である。

「術」は「置十八個開平方内十ニ個減以除五個乗乙円径」とあるが,十八は四十八,五は三の脱字,誤記,誤判読であろう。
つまり「置四十八個開平方内十ニ個減以除三個乗乙円径」が正しい。

また,算額の問の中の脱字「只云乙円径□□甲円径問幾何」は「一寸」であることもわかる(算額にはよくある設定)。

以上をまとめると,「乙円の直径が 1 寸のとき,甲円の直径は 1.690598923241497 寸である」。

r2 = 1/2
R = 4r2
r1 = R*(3 - √3)/3
2r1

   1.690598923241497

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   R = 4r2
   r1 = R*(3 - √3)/3
   @printf("乙円の直径が %g のとき,甲円の直径は %g である。\n", 2r2, 2r1)
   @printf("r2 = %g;  r1 = %g;  R = %g\n", r2, r1, R)
   (x0, y0) = R .* (cosd(120), sind(120))
   (x2, y2) = (R - r2) .* (cosd(30), sind(30))
   plot()
   circle(0, 0, R)
   circle(R*cosd(60), R*sind(60), R, beginangle=180, endangle=300)
   circle(-R*cosd(60), R*sind(60), R, beginangle=240, endangle=360)
   circle(R*cosd(60), -R*sind(60), R, beginangle=60, endangle=180)
   circle(-R*cosd(60), -R*sind(60), R, beginangle=0, endangle=120)
   circle22(0, R - r1, r1, :blue)
   circle4(x2, y2, r2, :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(x0, y0, "(x0,y0)", :red, :right, :bottom, delta=delta/2)
       point(0, R - r1, "甲円:r1,(0,R-r1)", :blue, :center, delta=-delta/2)
       point(x2, y2, "乙円:r2,(x2,y2)", :green, :center, delta=-delta/2)
   end
end;

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

算額(その1006)

2024年05月28日 | Julia

算額(その1006)

六三 羽生市須影 八幡神社 慶応元年(1865)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.
キーワード:円3個,等脚台形,対角線

等脚台形の中に対角線を引き,区画された領域に等円 3 個を入れる。下底(下頭),上底(上頭)がそれぞれ 2 寸,1 寸 5 分のとき,高さはいかほどか。

下底,上底,高さをそれぞれ 2a, 2b, h
等円の半径と中心座標を r, (r, r), (0, h - r)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a, b, h, r
eq1 = dist2(a, 0, -b, h, r, r, r)
eq2 = dist2(a, 0, -b, h, 0, h - r, r)
res = solve([eq1, eq2], (r, h))[6]  # 6 of 7

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

r = a*(a - 2*b)*(a + b - sqrt((a^2*(a - 2*b)^2*(a^2 + 2*a*b + b^2) + 4*b^2*(a^2 - b^2)^2)/(a^2*(a - 2*b)^2)))/(2*(a^2 - b^2))
h = 2*b*(-a^2 + b^2)/(a*(a - 2*b));

r はもう少し簡約化できる...と思って,まずルートの中を simplify() すると,以下のようになるが,結果を見るともっと簡約化できそう。

r = a*(a - 2b)*(a + b - sqrt((a + b)^2*(a^2 - 2a*b + 2b^2)^2/(a^2*(a - 2b)^2)))/(2(a^2 - b^2));

自動的にはこれ以上簡約化できないが,ルートの中は全部2乗になっているので,ルートを外せる。

r = a*(a - 2b)*(a + b + (a + b)*(a^2 - 2a*b + 2b^2)/(a*(a - 2b)))/(2(a^2 - b^2));

一旦展開してから simplify() すると,驚きの結果になる r = a - b。

r = r |> expand |> simplify
r |> println

   a - b

下底,上底がそれぞれ 2 寸,1 寸 5 分のとき,高さは 1.3125 寸である。

r(a => 2/2, b => 1.5/2) |> println
h(a => 2/2, b => 1.5/2) |> println

   0.250000000000000
   1.31250000000000

ちなみに,整数解としては,a = 9;  b = 6;  r = 3;  h = 20 のようなのがある。

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

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

算額(その1005)

2024年05月28日 | Julia

算額(その1005)

六一 越谷市下間久里 第六天 文久2年(1862)

埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

外円の中に円弧(1/3円),甲円,乙円それぞれ 3 個を入れる。乙円の直径が 1 寸 2 分のとき,外円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
1/3円の半径と中心座標を R, (x0, y0)
甲円の半径と中心座標を r1, (x1, y1)
乙円の半径と中心座標を r2, (x2, y2)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::poitive,  x0::poitive, y0::poitive,
     r1::poitive, x1::poitive, y1::poitive,
     r2::poitive, x2::poitive, y2::poitive
(s30, s150) = (Sym(30), Sym(150))
(x0, y0) = R .* (cosd(s150), sind(s150))
(x1, y1) = (R - r1) .* (cosd(s30), sind(s30))
(x2, y2) = (R - 2r1 - r2) .* (cosd(s30), sind(s30))
eq1 = (x1 - x0)^2 + (y1 - y0)^2 - (R + r1)^2
eq2 = (x2 - x0)^2 + (y2 - y0)^2 - (R + r2)^2
res = solve([eq1, eq2], (R, r1))[3]  # 3 of 3

   (85*r2/6, 17*r2/3)

外円の半径 R は乙円の半径 r2 の 85/6 倍である(甲円の半径 r1 は 17/3 倍)。
乙円の直径が 1 寸 2 分のとき,外円,甲円の直径はそれぞれ 17 寸,6 寸 8 分である。

1.2 .* (85/6, 17/3)

   (17.0, 6.8)

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

   x0 = -7.36122;  y0 = 4.25;  x1 = 4.41673;  y1 = 2.55;  x2 = 0.952628;  y2 = 0.55

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1.2/2
   (s30, s150) = (30, 150)
   (R, r1) = (85*r2/6, 17*r2/3)
   (x0, y0) = R .* (cosd(s150), sind(s150))
   (x1, y1) = (R - r1) .* (cosd(s30), sind(s30))
   (x2, y2) = (R - 2r1 - r2) .* (cosd(s30), sind(s30))
   @printf("乙円の直径が %g のとき,外円,甲円の直径はそれぞれ %g, %g である。\n", 2r2, 2R, 2r1)
   @printf("x0 = %g;  y0 = %g;  x1 = %g;  y1 = %g;  x2 = %g;  y2 = %g\n", x0, y0, x1, y1, x2, y2)
   plot()
   circle(0, 0, R, :magenta)
   rotate(x1, y1, r1)
   rotate(x2, y2, r2, :blue)
   circle(R*cosd(s30), R*sind(s30), R, beginangle=150, endangle=270, :green)
   circle(R*cosd(s150), R*sind(s150), R, beginangle=270, endangle=390, :green)
   circle(0, -R, R, beginangle=30, endangle=150, :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(x1, y1, "甲円:r1,(x1,y1)", :red, :center, delta=-delta/2)
       point(x2, y2, "   乙円:r2,(x2,y2)", :blue, :left, :vcenter)
       point(x0, y0, " 1/3円:R,(x0,y0)", :green, :left, :vcenter)
   end
end;

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

算額(その1004)

2024年05月27日 | Julia

算額(その1004)

二十 武州不動岡村 不動堂 文政元年(1818)

埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

楕円の中に甲円,乙円が 2 個ずつ入っている。楕円の長径,短径がそれぞれ 25 寸,15 寸のとき,甲円の直径はいかほどか。

楕円の長半径,短半径,中心座標を a, b, (0, 0)
甲円の半径と中心座標を r1, (x1, 0)
乙円の半径と中心座標を r2, (0, R - r2)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::poitive, b::poitive,
     r1::poitive, x1::poitive, r2::poitive
r2 = b/2
eq1 = x1^2 + r2^2 - (r1 + r2)^2
eq2 = (a^2 - b^2)*(b^2 - r1^2)/b^2 - x1^2  # 「算法助術」公式84
res = solve([eq1, eq2], (r1, x1))[3]

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

r1 の式は簡略化できる。

楕円の長径,短径がそれぞれ a, b のとき甲円の直径は (b - b^3/a^2) である。

res[1] |> simplify |> println

   b - b^3/a^2

楕円の長径,短径がそれぞれ 25 寸,15 寸のとき甲円の直径は 9.6 寸である。

(a, b) = (25, 15) ./ 2
(b - b^3/a^2, b*sqrt((a - b)*(a + b)*(2*a^2 - b^2))/a^2)

   (4.8, 7.683749084919418)

なお,例としては「楕円の長径,短径がそれぞれ 8, 4 のとき,甲円の直径は 3 である(乙円の直径は 2)」などがきれいである(上の図はこのときのもの)。

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

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

算額(その1003)

2024年05月27日 | Julia

算額(その1003)

十八 大里郡岡部村岡 稲荷社(久留里社算題集) 文化14年(1817)

埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

楕円の中に各辺が同じ長さの六角形を入れる。長径と短径が 14 寸,7 寸のとき,各辺の長さはいかほどか。

楕円の長半径,短半径,中心座標を a, b, (0, 0)
六角形の頂点の一つの座標を (x, y)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::poitive, b::poitive,
     x::poitive, y::poitive
eq1 = x^2/a^2 + y^2/b^2 - 1  # (x, y) が楕円の周上にある
eq2 = sqrt(x^2 + (b - y)^2) - 2y  # 2辺が同じ長さである
res = solve([eq1, eq2], (x, y))[2]

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

x 座標は 2*a*sqrt(b^4*(a^2 + 2*b^2)/(a^2 + 3*b^2)^2)/b
y 座標は b*(a^2 + b^2)/(a^2 + 3*b^2))
長半径,短半径,がそれぞれ 14/2, 7/2 のとき,x = 4.898979485566356, y = 2.5
また,六角形の一辺の長さは 5 である。

(a, b) = (14/2, 7/2)
x = 2*a*sqrt(b^4*(a^2 + 2*b^2)/(a^2 + 3*b^2)^2)/b
y = b*(a^2 + b^2)/(a^2 + 3*b^2)
(x, y) |> println
2y |> println

   (4.898979485566356, 2.5)
   5.0

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (14, 7) ./ 2
   x = 2*a*sqrt(b^4*(a^2 + 2*b^2)/(a^2 + 3*b^2)^2)/b
   y = b*(a^2 + b^2)/(a^2 + 3*b^2)
   @printf("(x, y) = (%g, %g)\n", x, y)
   println("sqrt(x^2 + (b - y)^2) = $(sqrt(x^2 + (b - y)^2))\n2y = $(2y)")
   plot([0, -x, -x, 0, x, x, 0], [b, y, -y, -b, -y, y, b], color=:blue, lw=0.5)
   ellipse(0, 0, a, b, color=:red)
   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(a, 0, " a", :red, :left, :bottom, delta=delta)
       point(0, b, " b", :red, :left, :bottom, delta=delta)
       point(x, y, " (x,y)", :blue, :left, :bottom, delta=delta)
   end
end;

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

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

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