裏 RjpWiki

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

算額(その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でシェアする

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

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