裏 RjpWiki

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

算額(その1059)

2024年06月13日 | Julia

算額(その1059)

九十四 大船渡市立根町 五葉神社 文政5年(1822)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

外円の中に正方形と小円 2 個ずつを容れる。小円の直径が 1 寸のとき,正方形の一辺の長さはいかほどか。

外円の半径と中心座標を R, (0, 0)
小円の半径と中心座標を r, (0, R - r)
とおき,以下の方程式を解き,外円の半径 R を求める。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive
eq1 = dist2(0, 0, R/2, R/2, 0, R - r, r)
R = solve(eq1, R)[1]
R |> println

   r*(1 + sqrt(2))

外円の半径 R は,小円の半径 r の (1 + √2) 倍である。
正方形の一辺の長さは (R/2)\*√2 = r\*(√2/2 + 1) である。
したがって 小円の半径が 1/2 のとき 0.853553390593274 である。

r = 1/2
r*((1 + √2)/2*√2) |> println
r*(√2/2 + 1) |> println

   0.8535533905932737
   0.8535533905932737

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

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

算額(その1058)

2024年06月13日 | Julia

算額(その1058)

九十四 大船渡市立根町 五葉神社 文政5年(1822)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

正方形の中に正三角形と大円,小円を容れる。小円の直径が 1 寸のとき,大円の直径はいかほどか。

正方形の一辺の長さを a
正方形の辺上にある正三角形の頂点の座標を (0, b), (b, 0)
大円の半径と中心座標を r1, (r1, r1)
小円の半径と中心座標を r2, (a - r2, r2), (r2, a - r2)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive
eq1 = a + (a - b) - sqrt(a^2 + (a - b)^2) - 2r2
eq2 = 2b - sqrt(Sym(2))b - 2r1
eq3 = a^2 + (a - b)^2 - 2b^2;
res = solve([eq1, eq2, eq3], (r1, a, b))[1]

   (2*r2/((-2 + sqrt(6))*sqrt(2*sqrt(2) + 3)), r2*(sqrt(2) + 2 + sqrt(12*sqrt(2) + 18))/((-2 + sqrt(6))*sqrt(2*sqrt(2) + 3)), 2*r2*(sqrt(2) + 2)/((-2 + sqrt(6))*sqrt(2*sqrt(2) + 3)))

@syms d
ans_r2 = res[1]/r2
ans_r2 = apart(ans_r2, d) |> sympy.sqrtdenest |> simplify |> (x -> x*r2)
ans_r2 |> println

   r2*(-sqrt(6) - 2 + 2*sqrt(2) + 2*sqrt(3))

大円の半径 r1 は,小円の半径 r2 の (2(√2 + √3 - 1) - √6) = 1.8430389971007672 倍である。
小円の直径が 1 寸のとき,大円の直径は 1.8430389971007672 寸である。

2(√2 + √3 - 1) - √6

   1.8430389971007672

その他のパラメータは以下のとおりである。
   正三角形の一辺の長さは 4.44949
   r2 = 0.5;  r1 = 0.921519;  a = 4.29788;  b = 3.14626

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

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

算額(その1057)

2024年06月13日 | Julia

算額(その1057)

九十四 大船渡市立根町 五葉神社 文政5年(1822)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.
算額(その1057)

直角三角形の中に,大中小の正方形,甲円,乙円,丙円を容れる。甲円,丙円の直径が 9 寸,4 寸のとき,乙円の直径はいかほどか。


三角形内の正方形の一辺の長さについては,算額(その356)に記した。
https://blog.goo.ne.jp/r-de-r/e/21d75d77f3ea4dff26d8a4e76e74659d

本問の図形において,小正方形,中正方形,大正方形の一辺の長さは等比数列をなすということである。また,それぞれの一辺を斜辺とする直角三角形も相似であり,それぞれに内接する甲円,乙円,丙円も相似でその直径も等比数列をなす。甲円径:乙円径:丙円径 = 9寸:x寸:4寸なので,乙円の直径は x = sqrt(9*4) = 6 寸である。

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

算額(その1056)

2024年06月13日 | Julia

算額(その1056)

九十一 陸前高田市小友町只出 小友八幡神社 明治26年(1893)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

外円の中に甲円,乙円,丙円,丁円が入っている。丙円の直径が 8 分 5 厘のとき,甲円の直径はいかほどか。

後にわかるが,丁円の大きさは丙円の大きさに無関係である。

算額(その938)の類題(発展問題)である(甲円,乙円と大円,中円の位置が逆)。算額(その938)では中央の円は脇の 4 円と同じであったが,本問では中央の円(丁円)は脇の円(丙円)とは別で,任意の大きさになりうる。

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

include("julia-source.txt");

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

   (2*r3, 2*r3 + r4, 4*r3 + r4, sqrt(2*r3 + r4)*sqrt(4*r3 + r4))

甲円の半径 r1 は,丙円 r3 の 2 倍である。
丙円の直径が 0.85 寸のとき,甲円の直径は 1.7 寸である。

丁円の大きさは,甲円の大きさに影響を与えない。
「問」ではこの辺を暗示している。「只云乙丁若干丙径八分五厘甲径幾何」と,丁径は何でもよいといっている。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 0.85/2
   r4 = 0.5/2
   (r1, r2, R, x3) = (2*r3, 2*r3 + r4, 4*r3 + r4, sqrt(2*r3 + r4)*sqrt(4*r3 + r4))
   @printf("丙円,丁円の直径が %g, %g 寸のとき,甲円の直径は %g 寸である。\n", 2r3, 2r4, 2r1)
   @printf("r3 = %g;  r4 = %g;  r1 = %g;  r2 = %g;  R = %g;  x3 = %g\n", r3, r4, r1, r2, R, x3)
   plot()
   circle(0, 0, R, :blue)
   circle22(0, R - r2, r2)
   circle22(0, R - r1, r1, :green)
   circle(0, 0, r4, :magenta)
   circle4(x3, r3, r3, :brown)
   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 - r2, "乙円:r2,(0,R-r2)", :red, :center, delta=-delta/2)
       point(0, R - r1, "甲円:r1,(0,R-r1)", :green, :center, delta=-delta/2)
       point(x3, r3, "丙円:r3\n(x3,r3)", :brown, :center, delta=-delta/2)
       point(0, R, " R", :blue, :left, :bottom, delta=delta/2)
       point(0, 0, "丁円:r4,(0,0)", :magenta, :center, delta=-delta/2)
   end
end;

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

算額(その1055)

2024年06月13日 | Julia

算額(その1055)

九十一 陸前高田市小友町只出 小友八幡神社 明治26年(1893)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

直角三角形内に正方形と大円,甲円,乙円,丙円,丁円を容れる。甲円,乙円の直径がそれぞれ 1.17 寸,1.05 寸のとき,大円の直径はいかほどか。

術では「甲円と乙円の直径を加えて 1 で割れば大円の直径が得られる」などとデタラメを言っている。山村もなんのコメントも付けずオウム返しで解説しているだけ。

直角三角形の直角を挟む二辺の短い方を「鈎」,長い方を「股」とする
直角三角形の斜辺が正方形の頂点で分割されるが,短い方を「短弦」,長い方を「長弦」とする
大円の半径と中心座標を r1, (r1, r1)
甲円の半径と中心座標を r2, (2r1 + r2, r2)
乙円の半径と中心座標を r3, (r3, 2r1 + r3)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive,
     短弦::positive, 長弦::positive,
     鈎::positive, 股::positive
eq1 = 2r1 + (股 - 2r1) - 長弦 - 2r2
eq2 = (鈎 - 2r1) + 2r1 - 短弦 - 2r3
eq3 = (鈎 - 2r1)^2 + 4r1^2 - 短弦^2
eq4 = 4r1^2 + (股 - 2r1)^2 - 長弦^2
eq5 = 鈎^2 + 股^2 - (長弦 + 短弦)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (r1, 短弦, 長弦, 鈎, 股))[2]  # 2 0f 2

   (r2/2 + r3/2 + sqrt(r2^2 + r3^2)/2, r2 + sqrt(r2^2 + r3^2) + r3^2/r2 + r3*sqrt(r2^2 + r3^2)/r2, r2^2/r3 + r2*sqrt(r2^2 + r3^2)/r3 + r3 + sqrt(r2^2 + r3^2), r2 + 2*r3 + sqrt(r2^2 + r3^2) + r3^2/r2 + r3*sqrt(r2^2 + r3^2)/r2, (r2^2 + 2*r2*r3 + r2*sqrt(r2^2 + r3^2) + r3^2 + r3*sqrt(r2^2 + r3^2))/r3)

res[1] |> println
res[1](r2 => 1.17, r3 => 1.05) |> println

   r2/2 + r3/2 + sqrt(r2^2 + r3^2)/2
   1.89603435039443

大円の半径は (r2 + r3 + sqrt(r2^2 + r3^2))/2 である。術がいうような単純な和ではない。
甲円,乙円の直径がそれぞれ 1.17 寸,1.05 寸のとき,大円の直径は 1.89603435039443 寸である。

その他のパラメータは以下のとおりである。
r2 = 1.17;  r3 = 1.05;  r1 = 1.89603;  短弦 = 5.09521;  長弦 = 5.67752;  鈎 = 7.19521;  股 = 8.01752

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3) = (1.17, 1.05)
   (r1, 短弦, 長弦, 鈎, 股) = (r2/2 + r3/2 + sqrt(r2^2 + r3^2)/2, r2 + sqrt(r2^2 + r3^2) + r3^2/r2 + r3*sqrt(r2^2 + r3^2)/r2, r2^2/r3 + r2*sqrt(r2^2 + r3^2)/r3 + r3 + sqrt(r2^2 + r3^2), r2 + 2*r3 + sqrt(r2^2 + r3^2) + r3^2/r2 + r3*sqrt(r2^2 + r3^2)/r2, (r2^2 + 2*r2*r3 + r2*sqrt(r2^2 + r3^2) + r3^2 + r3*sqrt(r2^2 + r3^2))/r3)
   @printf("甲円,乙円の直径が %g, %g のとき,大円の直径は %g である。\n", 2r2, 2r3, 2r1)
   @printf("r2 = %g;  r3 = %g;  r1 = %g;  短弦 = %g;  長弦 = %g;  鈎 = %g;  股 = %g\n", r2, r3, r1, 短弦, 長弦, 鈎, 股)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
   plot!([0, 2r1, 2r1, 0, 0], [0, 0, 2r1, 2r1, 0], color=:magenta, lw=0.5)
   circle(r1, r1, r1, :green)
   circle(2r1 + r2, r2, r2, :orange)
   circle(r3, 2r1 + r3, r3, :purple)
   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, " 股", :blue, :left, :bottom, delta=delta)
       point(0, 鈎, " 鈎", :blue, :left, :bottom, delta=delta)
       point(r1, r1, "大円:r1,(r1,r1)", :green, :center, :bottom, delta=delta)
       point(2r1 + r2, r2, "甲円:r2\n(2r1+r2,r2)", :orange, :center, :bottom, delta=delta)
       point(r3, 2r1 + r3, "乙円:r3\n(r3,2r1+r3)", :purple, :center, :bottom, delta=delta)
   end
end;

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

2024/06/13

2024年06月13日 | 写真

これは何の花か分かりますか?

答えは、ゴボウの花


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

算額(その1054)

2024年06月13日 | Julia

算額(その1054)

八十九 陸前高田市小友町 常膳寺観音堂 天保13年(1842)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

全円の中に正三角形と等円 5 個を容れる。全円の直径が 1 寸のとき,等円の直径はいかほどか。

全円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (r, -r - R/2)
とおき,方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive
eq1 = r^2 + (r + R/2)^2 - (R - r)^2
res = solve(eq1, r)[1]
res |> println
2res(R => 1/2).evalf()|> println

   R*(-3 + 2*sqrt(3))/2
   0.232050807568877

等円の半径 r は,全円の半径 R の (2√3 - 3)/2 倍である。
全円の直径が 1 寸のとき,等円の直径は 0.2320508075688772 寸である。

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

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

算額(その1053)

2024年06月12日 | Julia

算額(その1053)

八十八 陸前高田市小友町 常膳寺観音堂 天保12年(1841)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

この算額は,「問」,「答」,「術」とも読み取ることができないそうだ。図形のみ認識できるようなので,いくつかは再現できそうである。

以下の図形から,「外円の中に正方形と大円,小円を容れる。小円の直径が 97 寸のとき,大円の直径はいかほどか」という問が可能であろう。

外円の半径と中心座標を R, (0, 0); R = 2r1
大円の半径と中心座標を r1, (r1, 0)
小円の半径と中心座標を r2, (0, y2)
とおき,連立方程式を解く。
ここでは,小円の半径が既知として解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive,
     r2::positive, y2::positive
R = 2r1
eq1 = r1^2 + y2^2 - (r1 + r2)^2
eq2 = dist2(R, 0, 0, R, 0, y2, r2)
res = solve([eq1, eq2], (r1, y2))[1];

大円の半径
res[1] |> simplify |> println

   r2*(1 + 2*sqrt(2) + sqrt(5 + 4*sqrt(2)))/4

小円の中心の y 座標
res[2] |> simplify |> println

   r2*(1 + sqrt(5 + 4*sqrt(2)))/2

大円の半径 r1 は,小円の半径 r2 の (1 + 2√2 + sqrt(5 + 4√2))/4 倍である。
小円の直径が 97 寸のとき,大円の直径は 97*(1 + 2√2 + sqrt(5 + 4√2))/4 = 172.00308906722228 = 172 寸有奇である。

その他のパラメータは以下のとおりである。
r2 = 48.5;  r1 = 86.0015;  y2 = 103.414;  R = 172.003

97*(1 + 2√2 + sqrt(5 + 4√2))/4

   172.00308906722228

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 97/2
   (r1, y2) = r2 .* ((1 + 2√2 + sqrt(5 + 4√2))/4, (1 + sqrt(5 + 4√2))/2)
   R = 2r1
   @printf("小円の直径が %g のとき,大円の直径は %g である。\n", 2r2, 2r1)
   @printf("r2 = %g;  r1 = %g;  y2 = %g;  R = %g\n", r2, r1, y2, R)
   plot([R, 0, -R, 0, R], [0, R, 0, -R, 0], color=:blue, lw=0.5)

   circle(0, 0, R, :magenta)
   circle2(r1, 0, r1)
   circle22(0, 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(R, 0, " R", :blue, :left, :bottom, delta=delta)
       point(0, R, " R", :green, :left, :bottom, delta=delta)
       point(r1, 0, "大円:r1,(1,0)", :red, :center, :bottom, delta=delta)
       point(0, y2, "小円:r2,(0,y2) ", :black, :center, delta=-delta)
   end
end;

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

算額(その1052)

2024年06月12日 | Julia

算額(その1052)

八十八 陸前高田市小友町 常膳寺観音堂 天保12年(1841)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

この算額は,「問」,「答」,「術」とも読み取ることができないそうだ。図形のみ認識できるようなので,いくつかは再現できそうである。

以下の図形から,「正三角形内に正方形を入れ,正方形の1つの頂点から正三角形の一つの頂点を結ぶ斜線を引き,分割された領域に大円,小円を 1 個ずつ容れる。正三角形の一辺の長さが5寸のとき,正方形の一辺の長さ,大円,小円の直径を求めよ。」という問が可能であろう。

正三角形の一辺の長さを 2a
正方形の一辺の長さを 2b
大円の半径と中心座標を r1, (b - r1, 2b - r1)
小円の半径と中心座標を r2, (b + r2, y2)
とおき,連立方程式を解く。
未知数は a, b, r1, r2, y2 の 5 個なので,条件式が 4 個あればどれか一つを既知として解を求めることが可能である。
ここでは,正三角形の一辺の長さを既知として解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r1::positive,
     y1::positive, r2::positive, y2::positive
eq1 = dist2(-b, 2b, a, 0, b - r1, 2b - r1, r1)
eq2 = dist2(-b, 2b, a, 0, b + r2, y2, r2)
eq3 = dist2(0, √Sym(3)a, a, 0, b + r2, y2, r2)
eq4 = (√Sym(3)a - 2b)/b - √Sym(3);
#solve([eq1, eq2, eq3, eq4], (b, r1, r2, y2))

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)
   (b, r1, r2, y2) = u
   return [
       4*b*(-2*a*b*r1 + a*r1^2 + 4*b^3 - 6*b^2*r1 + b*r1^2),  # eq1
       4*a^2*b^2 - 4*a^2*b*y2 - a^2*r2^2 + a^2*y2^2 - 8*a*b^3 - 8*a*b^2*r2 - 2*a*b*r2^2 + 4*a*b*r2*y2 + 2*a*b*y2^2 + 4*b^4 + 8*b^3*r2 + 4*b^3*y2 - b^2*r2^2 + 4*b^2*r2*y2 + b^2*y2^2,  # eq2
       3*a^2/4 - 3*a*b/2 - 3*a*r2/2 - sqrt(3)*a*y2/2 + 3*b^2/4 + 3*b*r2/2 + sqrt(3)*b*y2/2 - r2^2/4 + sqrt(3)*r2*y2/2 + y2^2/4,  # eq3
       2*b/(a - b) - sqrt(3),  # eq4
   ]
end;

a = 262/2
iniv = BigFloat[7, 3.1, 2.1, 6] .* a/15
res = nls(H, ini=iniv)

   ([60.79731158304585, 27.355505309267038, 18.002338812484446, 54.40898006283059], true)

「正三角形の一辺の長さが 262 寸のとき,小円の直径は 36.00寸有奇 である」というと,算額らしい。

正三角形の一辺の長さが 262 のとき,正方形の一辺の長さは 121.595,大円,小円の直径は 54.711, 36.0047 である。
a = 131;  b = 60.7973;  r1 = 27.3555;  r2 = 18.0023;  y2 = 54.409

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   #a = 15
   (b, r1, r2, y2) = res[1]
   @printf("正三角形の一辺の長さが %g のとき,正方形の一辺の長さは %g,大円,小円の直径は %g, %g である。\n", 2a, 2b, 2r1, 2r2)
   @printf("a = %g;  b = %g;  r1 = %g;  r2 = %g;  y2 = %g\n", a, b, r1, r2, y2)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:blue, lw=0.5)
   plot!([b, b, -b, -b, b], [0, 2b, 2b, 0, 0], color=:green, lw=0.5)
   segment(-b, 2b, a, 0, :orange)
   circle(b - r1, 2b - r1, r1, :red)
   circle(b + r2, y2, r2, :gray70)
   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)
       point(0, √3a, " √3a", :blue, :left, :bottom, delta=delta)
       point(b, 0, " b", :green, :left, :bottom, delta=delta)
       point(0, 2b, " 2b", :green, :left, :bottom, delta=delta)
       point(b - r1, 2b - r1, "大円:r1\n(b-r1,2b-r1)", :red, :center, :bottom, delta=delta)
       point(b + r2, y2, "小円:r2,(b+r2,y2) ", :black, :right, :vcenter)
   end
end;

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

算額(その1051)

2024年06月12日 | Julia

算額(その1051)

八十五 室根村折壁字大洞 入沢弥栄神社 明治16年(1883)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

この算額は,「問」,「答」,「術」のいずれも不適切極まりないものである。
算額の図から以下のように解釈して解を求めることにする。

この問題は,乙円の直径を求めることと,等円の直径を求めることが独立であり,それぞれ別々に解を求めるべきものである。

1. 外円内に甲円,乙円 3 個ずつ容れる。甲円の直径が 18寸5分のとき,乙円の直径はいかほどか。

算額(その27)
https://blog.goo.ne.jp/r-de-r/e/a729e33453ca3b7f28ed787cbb22faa9

大円の半径と中心座標を R, (0, 0); R = 2r1
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (x2, y2); x2 = (R - r2)*cosd(30), y2 = (R - r2)*sind(30)
とおき,方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive,
     x2::positive, y2::positive, r3::positive
R = 2r1
x2 = (R - r2)*cosd(Sym(30))
y2 = (R - r2)*sind(Sym(30))
eq1 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
res1 = solve(eq1, r2)[1]
res1 |> println
res1(r1 => 18.5) |> println

   2*r1/5
   7.40000000000000

乙円の半径 r2 は,甲円の半径 r1 の 2/5 倍である。
甲円の直径が 18.5 寸のとき,乙円の直径は 18.5*2/5 = 7.4 寸である。

術では「大円径三十七寸七分,乙円径五十分三十七」とある。いずれも不適切な解であるが,算額1052 のときと同じように,径を述べるときに実寸(三十七寸七分)と割合(五十分三十七)を混用する意味もわからない。村山の註も何を言っているのか。

2. 外円の外側に,外円に外接し,かつ隣同士概説する 28 個の等円を描くとき,等円の直径はいかほどか。

算額(その618)
https://blog.goo.ne.jp/r-de-r/e/ebb9afbef32e1f62bc9584a730d2ca83

等円の直径を聞いているのに,「等円径各二十五分十二」と書いている。どういうことだろうか。しかもここでも,割合(二十五分十二)を示している。

ここでは題意を汲んで,直径 18.5 寸の外円の外に並ぶ 28 個の等円の直径を求めることにする。

等円の半径と中心座標を r3; r3 = (R + r3)*sind(360/(28*2))
とおき,方程式を解く。

eq2 = r3 - (R + r3)*sind(Sym(360)//(Sym(28)*2))
res2 = solve(eq2, r3)[1]
res2 |> println
res2(r1 => 18.5).evalf() |> println

   2*r1*sin(pi/28)/(1 - sin(pi/28))
   4.66499988383834

等円の半径 r3 は,甲円の半径 r1 の 2sin(π/28)/(1 - sin(π/28)) 倍である。
甲円の直径が 18.5 寸のとき,等円の直径は 18.5*2sin(π/28)/(1 - sin(π/28)) = 4.664999883838343 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 18.5/2
   R = 2r1
   r2 = 2r1/5
   r3 = r1*2sin(π/28)/(1 - sin(π/28))
   @printf("甲円 = %g;  乙円 = %g;  等円 = %g;  大円 = %g\n", 2r1, 2r2, 2r3, 2R)
   plot()
   circle(0, 0, R, :green)
   circle(0, 0, R + 2r3, :green)
   rotate(0, R - r1, r1)
   x2 = (R - r2)*cosd(30)
   y2 = (R - r2)*sind(30)
   rotate(x2, y2, r2, :blue)
   rotate(0, R + r3, r3, :orange, angle=360/28)
   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)
       二十八宿 = [
          "角", "亢", "氐", "房", "心", "尾", "箕",  # 東方青龍
          "斗", "牛", "女", "虚", "危", "室", "壁",  # 北方玄武
          "奎", "婁", "胃", "昴", "畢", "觜", "参",  # 西方白虎
          "井", "鬼", "柳", "円", "張", "翼", "軫"   # 南方朱雀
       ]
       for i = 1:28
           x0 = (R + r3)*cosd(360/28*(mod(15 - i, 28)+1))
           y0 = (R + r3)*sind(360/28*(mod(15 - i, 28)+1))
           point(x0, y0, 二十八宿[i], :blue, :center, :vcenter, mark=false)
       end
       point(0, R, " R", :green, :left, delta=-delta/2)
       point(R, 0, "R ", :green, :right, :bottom, delta=delta/2)
       point(R+2r3, 0, " R+2r3", :green, :left, :bottom, delta=delta/2)
       point(0, R - r1, "甲円:r1,(0,R-r1)", :red, :center, delta=-delta/2)
       point(x2, y2, "乙円:r2\n(x2,y2)", :blue, :center, delta=-delta/2)
       plot!(xlims=(-R - 2r3 - 2delta, R + 2r3 + 8delta))
   end
end;

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

算額(その1050)

2024年06月12日 | Julia

算額(その1050)

八十五 室根村折壁字大洞 入沢弥栄神社 明治16年(1883)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

外円の中に,甲円,乙円,丙円,丁円を容れる。甲円の直径は十二寸五分,「甲乙円径五分三」(注1)のとき,丙円の直径はいかほどか。

注1:これは,「乙円の直径は甲円の直径の 3/5」であろう。ほかに解釈のしようがない。乙円の直径は七寸五分である。しかしこの時点で,算額のような図は描けないことが明らかである。下図の A が直径が 12.5 寸の甲円である。

甲円がどの程度互いに交差するのかは変化しうるが。乙円 B は直径が 7.5 寸である。乙円が外円に接するように描いているが,甲円と乙円は交わってしまう。甲円の交わりが大きくなると,丙円との交わりも大きくなってしまう。逆に甲円の交わりが小さくなり,交わりが 0 になってもまだ甲円と乙円は交わったままである。甲円の交わりがなくなって(中心が離れて)は,算額の図のようにはならない。
なぜこのようなことになるかといえば,それは乙円の大きさを固定してしまったからである。実際以下のように連立方程式を立ててみると,条件が 1 個多すぎる。それは乙円を既知とするからである。そこで,以下では乙円の直径は未知として解を求めることにする。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x31, 0), (x32, y3)
丁円の半径と中心座標を r4, (0, 0)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive,
     x2::positive, r3::positive, x31::positive,
     x32::positive,  y3::positive, r4::positive
eq1 = x2^2 + r2^2 - (R - r2)^2
eq2 = x32^2 + y3^2 - (R - r3)^2
eq3 = x2^2 + (R - r1 - r2)^2 - (r1 + r2)^2
eq4 = x31^2 + (R - r1)^2 - (r1 + r3)^2
eq5 = x32^2 + (y3 - R + r1)^2 - (r1 + r3)^2
eq6 = (x2 - x32)^2 + (y3 - r2)^2 - (r2 + r3)^2
eq7 = (x2 - x31)^2 + r2^2 - (r2 + r3)^2
eq8 = 2r1 - r4 - R;
## res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (R, r2, x2, r3, x31, x32, y3, r4))

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)
   #(R, x2, r3, x31, x32, y3, r4) = u
   (R, r2, x2, r3, x31, x32, y3, r4) = u
   return [
       r2^2 + x2^2 - (R - r2)^2,  # eq1
       x32^2 + y3^2 - (R - r3)^2,  # eq2
       x2^2 - (r1 + r2)^2 + (R - r1 - r2)^2,  # eq3
       x31^2 + (R - r1)^2 - (r1 + r3)^2,  # eq4
       x32^2 - (r1 + r3)^2 + (-R + r1 + y3)^2,  # eq5
       (-r2 + y3)^2 - (r2 + r3)^2 + (x2 - x32)^2,  # eq6
       r2^2 - (r2 + r3)^2 + (x2 - x31)^2,  # eq7
       -R + 2*r1 - r4,  # eq8
   ]
end;

r1 = 12.5/2
iniv = BigFloat[1.9, 0.45, 1.38, 0.24, 0.85, 1.22, 1.12, 0.099] .*r1
res = nls(H, ini=iniv)

   ([11.881752126800865, 2.815876063400433, 8.617479375809692, 1.5081750129387437, 5.335976715996894, 7.63184686468951, 7.02609532892722, 0.6182478731991338], true)

甲円の直径が 12.5 のとき,乙円の直径は 5.63175212680087(注2),丙円の直径は 3.01635002587749(注3)である。

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

   r1 = 6.25;  r2 = 2.81588;  r3 = 1.50818;  r4 = 0.618248;  R = 11.8818
   0.45054017014406933

注2:「問」では乙円の直径は甲円の直径の 3/5=0.6 と言っていたが,実際には約 0.45 である。先に示した図でも,乙円の直径が言っているより小さければ,甲円に外接し外円に内接するという条件を満たすことは自明であった。

注3:「答」では,「乙円径五十分十九」といっている。五十分の十九寸でないことは図を見ても明らかである。また,山村は術の解説で,「丙円径 = 19/50」としたあと,その算出過程を示しているが,最後に「27.8 - 20)÷22=3.9, 3.9=19.5/50,註 答に 19/50 とあり,小差がある」等と書いているが,「3.9=19.5/50」はひどい。19.5/50 は 0.39 である。0.39 寸でも 0.38 寸でもどちらも間違っている。小さすぎる。
では,別の解釈は丙円の直径はなにかの 19/50 だろうが,該当するようなものは見当たらない。「問」に続いて素直に読めば,「丙円は甲円の 19/50」と言いたいのだろうが,物差しで測ってみると「丙円径/甲円径≒1/6」のようだ(計算結果では「丙円径/甲円径≒1/4」)。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 12.5/2
   (R, r2, x2, r3, x31, x32, y3, r4) = res[1]
   @printf("甲円の直径が %g のとき,乙円の直径は %.15g,丙円の直径は %.15g である。\n", 2r1, 2r2, 2r3)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  r4 = %g;  R = %g\n", r1, r2, r3, r4, R)
   println(r2/r1)
   plot()
   circle(0, 0, R)
   circle22(0, R - r1, r1, :blue)
   circle4(x2, r2, r2, :green)
   circle4(x32, y3, r3, :magenta)
   circle2(x31, 0, r3, :magenta)
   circle(0, 0, r4, :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, " R", :red, :left, :bottom, delta=delta/2)
       point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(0, R - r1, "甲円:r1,(0,R-r1)", :blue, :center, delta=-delta)
       point(x2, r2, "乙円:r2\n(x2,r2)", :green, :center, delta=-delta)
       point(x31, 0, "丙円:r3\n(x31,r3)", :magenta, :center, :vcenter)
       point(x32, y3, "丙円:r3\n(x32,r3)", :magenta, :center, :vcenter)
       point(0, 0, "丁円:r4,(0,0)", :black, :center, delta=-2delta)
   end
end;

 

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

算額(その1049)

2024年06月12日 | Julia

算額(その1049)

八十五 室根村折壁字大洞 入沢弥栄神社 明治16年(1883)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

大円の内外に等円を 4 個描く。内外の 3 個ずつの等円は共通接線を持つ。大円の直径が 30 寸のとき,等円の直径はいかほどか。

大円の半径と中心座標を r1, (0, 0)
等円の半径と中心座標を r2, (0, y2), (r1 + r2, 0); y2 = r1 - r2
とおき,以下の方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, y2::positive,
     len::positive, θ::positive
len = sqrt((r1 + r2)^2 - r2^2)
sinθ = r2/(r1 + r2)
cosθ = sqrt(1 - (r2/(r1 + r2))^2)
y2 = r1 - r2
eq = y2*cosθ - r2
solve(eq, r2)[1] |> println

   r1*sqrt(-2 + sqrt(5))

等円の半径 r2 は 大円の半径の sqrt(√5 - 2) = 0.48586827175664576 倍である。
大円の直径が 30 寸のとき,等円の直径は 30*sqrt(√5 - 2) = 14.576048152699373 寸である。

r2 = 7.28802407634969;  y2 = 7.71197592365031

「術」では,等円径 = (大円径/2 + 0.15)*(√5 - √2) ≒ (大円径/2 + 0.15)*0.821 としているようだが,全くでたらめな解であるし,0.15 が何なのかの説明もない。山村もなんのコメントもしていない。
図を見ても(描いても),等円の直径は大円の半径に極めて近いことがわかる。大円径が 15 寸で,等円径が 12.438 寸は明らかにおかしいことになぜ気づかないのか。

大円径 = 30
等円径 = (大円径/2 + 0.15)*(√5 - √2)
等円径2 = (大円径/2 + 0.15)*0.821
(等円径, 等円径2)

   (12.451094389169425, 12.43815)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 30/2
   r2 = r1*sqrt(√5 - 2)
   y2 = r1 - r2
   @printf("大円の直径が %g のとき,等円の直径は %g である。\n", 2r1, 2r2)
   @printf("r2 = %.15g;  y2 = %.15g\n", r2, y2)
   plot()
   circle(0, 0, r1)
   circle2(r1 + r2, 0, r2, :blue)
   circle22(0, y2, r2, :blue)
   len = sqrt((r1 + r2)^2 - r2^2)
   sinθ = r2/(r1 + r2)
   cosθ = sqrt(1 - (r2/(r1 + r2))^2)
   (x0, y0) = (len*cosθ, len*sinθ)
   abline(x0, y0, y0/x0, -r1 - 2r2, r1 + 2r2, :green)
   abline(x0, -y0, -y0/x0, -r1 - 2r2, r1 + 2r2, :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, y0, " y0 = r1 - r2", :blue, :center, delta=-delta/2)
       point(0, y0, "等円:r2", :blue, :center, :bottom, delta=delta)
       point(0, r1, " r1", :red, :left, :bottom, delta=delta/2)
       point(r1 + r2, 0, "等円:r2\n(r1+r2,0)", :blue, :center, delta=-delta)
       
   end
end;

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

算額(その1048)

2024年06月12日 | Julia

算額(その1048)

八十二 藤沢町保呂羽保呂羽山 保呂羽神社 明治7年(1874)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

正方形の中に斜線(天斜,地斜と名付ける)を入れ,正三角形を作る。
天斜,地斜がそれぞれ 5 寸,7 寸のとき,正三角形の一辺の長さはいかほどか。

天斜,地斜の交点座標などを図のように定め,以下の連立方程式を解く。
正三角形の一辺の長さは sqrt((x1 - x3)^2 + (a - y1)^2) などである。

include("julia-source.txt");

using SymPy
@syms x1::positive, y1::negative,
     x2::positive, y2::positive,
     x3::positive, a::positive,
     天斜::positive, 地斜::positive
eq1 = (x1 - x3)^2 + (a - y1)^2 - (x3 - x2)^2 - (a - y2)^2 |> expand
eq2 = (x1 - x2)^2 + (y2 - y1)^2 - (x3 -x2)^2 - (a - y2)^2 |> expand
eq3 = a/(a - x3) - y1/(a - x1)
eq3 = a*(a - x1) - y1*(a - x3) |> expand
eq4 = (a - y2)/x2 - (a - y1)/x1
eq4 = (a - y2)*x1 - (a - y1)*x2 |> expand
eq5 = x1^2 + (a - y1)^2 - 天斜^2 |> expand
eq6 = (a - x3)^2 + a^2 - 地斜^2 |> expand;

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)
   (x1, y1, x2, y2, x3, a) = u
   return [
       -2*a*y1 + 2*a*y2 + x1^2 - 2*x1*x3 - x2^2 + 2*x2*x3 + y1^2 - y2^2,  # eq1
       -a^2 + 2*a*y2 + x1^2 - 2*x1*x2 + 2*x2*x3 - x3^2 + y1^2 - 2*y1*y2,  # eq2
       a^2 - a*x1 - a*y1 + x3*y1,  # eq3
       a*x1 - a*x2 - x1*y2 + x2*y1,  # eq4
       a^2 - 2*a*y1 + x1^2 + y1^2 - 天斜^2,  # eq5
       2*a^2 - 2*a*x3 + x3^2 - 地斜^2,  # eq6
   ]
end;

(天斜, 地斜) = [5, 7]

iniv = BigFloat[5, 6, 4, 6, 5, 7]
res = nls(H, ini=iniv)

   ([4.884718333076426, 5.602160370771752, 3.7901899873037372, 5.841354083664568, 4.544601991980781, 6.669646579819466], true)

天斜,地斜がそれぞれ 5, 7 のとき,正三角形の一辺の長さは 1.12035973329433 である。

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

   x1 = 4.88472;  y1 = 5.60216;  x2 = 3.79019;  y2 = 5.84135;  x3 = 4.5446;  a = 6.66965

(x1, y1, x2, y2, x3, a) = res[1]
sqrt((x1 - x3)^2 + (a - y1)^2)

   1.120359733294333

「術」は解析解を示している。

天斜 = 5
地斜 = 7
A = (√12天斜 + 地斜)*地斜 - 3天斜^2
一辺の長さ = ((天斜 + 地斜) - √A)/2

   1.1203597332943334

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x1, y1, x2, y2, x3, a) = res[1]
   @printf("天斜,地斜がそれぞれ %g, %g のとき,正三角形の一辺の長さは %.15g である。\n", 天斜, 地斜, sqrt((x1 - x3)^2 + (a - y1)^2))
   @printf("x1 = %g;  y1 = %g;  x2 = %g;  y2 = %g;  x3 = %g;  a = %g\n", x1, y1, x2, y2, x3, a)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:blue, lw=0.5)
   segment(0, a, x1, y1, :red)
   segment(a, 0, x3, a, :red)
   segment(x3, a, x2, y2, :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(0, a, " a", :blue, :left, :bottom, delta=delta/2)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(x1, y1, " (x1,y1)", :red, :left, :vcenter)
       point(x2, y2, "(x2,y2) ", :red, :right, :vcenter, delta=-delta)
       point(x3, a, "(x3,a)", :red, :center, :bottom, delta=delta/2)
   end
end;

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

算額(その1047)

2024年06月11日 | Julia

算額(その1047)

八十二 藤沢町保呂羽保呂羽山 保呂羽神社 明治7年(1874)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

直線の上に小円 2 個が互いに接して載っている。その上に,小円を 3 個含む,交差する大円 2 個が載っている。
黒積(図で灰色に塗った面積)を求めよ。



手計算すればよいのだが,大円とy 軸の交点の y 座標を求める。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::negative, x::positive, y::positive
r1 = 2r2
eq1 = (x - r2)^2 + (y - 4r2)^2 - 4r2^2
solve(eq1, y)[2](x => 0) |> println

   sqrt(3)*(-r2) + 4*r2

求める面積(黒積)は「台形OABC の面積から大円の面積の 1/12(OAD) と小円の面積の 1/4(CBD) を差し引いた面積の 2 倍」である。

# 台形OABC の面積
S1 = ((1 + (2 - √Sym(3))) + 3)/2*r2^2
S1 |> println
# 大円の面積の 1/12
S2 = PI*(2r2)^2/12
S2 |> println
# 小円の面積の 1/4 
S3 = PI*r2^2/4
S3 |> println

   r2^2*(3 - sqrt(3)/2)
   pi*r2^2/3
   pi*r2^2/4

小円の直径が 1 寸のときの黒積は -7*π/24 - sqrt(3)/4 + 3/2 = 0.150689440810758 である。

black = 2(S1 - S2 - S3)(r2 => 1//2)
black |> println
black.evalf() |> println

   -7*pi/24 - sqrt(3)/4 + 3/2
   0.150689440810758

function paintitblack(r2)
   θ = 240:0.1:270
   x = 2r2.*cosd.(θ) .+ r2
   y = 2r2.*sind.(θ) .+ 4r2
   θ2 = 90:0.1:180
   append!(x, r2.*cosd.(θ2) .+ r2)
   append!(y, r2.*sind.(θ2) .+ r2)
   append!(x, -reverse(x))
   append!(y, reverse(y))
   append!(x, [0, 0])
   append!(y, [r2, r2*(3 - √3)])
   plot!(x, y, seriestype=:shape, fillcolor=:gray80)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   r1 = 2r2
   plot(r2.*[1, 0, 0, 1, 1], r2.*[4, 4 - √3, 1, 1, 4], color=:black, lw=1)
   circle2(r2, 4r2, r1)
   circle(0, 4r2, r2, :blue)
   circle2(r2, r2, r2, :blue)
   circle2(2r2, 4r2, r2, :blue)
   paintitblack(r2)
   segment(-3r2, 0, 3r2, 0, :black, lw=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(r2, 4r2, "O ", :black, :right, :vcenter)
       point(0, r2*(4 - √3), "A ", :black, :center, :bottom, delta=delta)
       point(0, r2, "B ", :black, :right, :vcenter)
       point(r2, r2, " C", :black, :left, :vcenter)
       point(r2, 2r2, " D", :black, :left, :bottom, delta=delta)
       point(0, 4*r2 - sqrt(3)*r2)
       point(0, 4r2, "小円:r2,(0,4r2)", :blue, :center, :bottom, delta=delta)
       point(r2, r2, "小円:r2,(r2,r2)", :blue, :center, delta=-delta)
       point(r2, 4r2, " 大円:2r2,(r2,4r2)", :red, :left, :vcenter)
   end
end;

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

算額(その1046)

2024年06月11日 | Julia

算額(その1046)

八十二 藤沢町保呂羽保呂羽山 保呂羽神社 明治7年(1874)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

正三角形の中に全円,大円,小円と斜線を容れる。小円の直径が 1 寸のとき,大円の直径はいかほどか。

正三角形の一辺の長さを 2a
全円の半径と中心座標を r0, (0, r0); r0 = √3a/3
大円の半径と中心座標を r1, (x1, r1); x1 < 0
小円の半径と中心座標を r2, (x2, r2)
斜線と斜辺の交点座標を (b, √3(b + a)); b < 0
とおき,以下の連立方程式を解く。

一度に解けないので,逐次解いてゆく。

include("julia-source.txt");

using SymPy
@syms a::positive, b::negative, r0::positive,
     r1::positive, x1::negative, r2::positive, x2::positive,
     y2::positive
r0 = √Sym(3)a/3
eq1 = dist2(a, 0, 0, √Sym(3)a, x2, y2, r2)
eq2 = dist2(-a, 0, 0, √Sym(3)a, x1, r1, r1)
eq3 = dist2(a, 0, b, √Sym(3)*(b + a), x2, y2, r2)
eq4 = dist2(a, 0, b, √Sym(3)*(b + a), x1, r1, r1)
eq5 = x2^2 + (r0 - y2)^2 - (r0 + r2)^2
eq6 = x1^2 + (r0 - r1)^2 - (r0 + r1)^2;

eq2, eq6 から a, x1 を求める。

res = solve([eq2, eq6], (a, x1))[2]

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

eq1, eq3, eq5, eq6 に a, x1 を代入する。

eq1 = 4eq1(a => res[1]);
eq3 = numerator(eq3(a => res[1]));
eq4 = numerator(eq4(a => res[1], x1 => res[2]));
eq5 = eq5(a => res[1]);

eq1 |> println
eq3 |> println
eq4 |> println
eq5 |> println

   81*r1^2 - 18*sqrt(3)*r1*x2 - 18*r1*y2 - 4*r2^2 + 3*x2^2 + 2*sqrt(3)*x2*y2 + y2^2
   81*b^2*r1^2 - 18*sqrt(3)*b^2*r1*x2 + 18*b^2*r1*y2 - 4*b^2*r2^2 + 3*b^2*x2^2 - 2*sqrt(3)*b^2*x2*y2 + b^2*y2^2 + 486*sqrt(3)*b*r1^3 - 324*b*r1^2*x2 - 12*sqrt(3)*b*r1*r2^2 + 18*sqrt(3)*b*r1*x2^2 - 6*sqrt(3)*b*r1*y2^2 + 2187*r1^4 - 486*sqrt(3)*r1^3*x2 - 486*r1^3*y2 - 108*r1^2*r2^2 + 81*r1^2*x2^2 + 54*sqrt(3)*r1^2*x2*y2 + 27*r1^2*y2^2
   252*b^2*r1^2 + 1332*sqrt(3)*b*r1^3 + 5184*r1^4
   x2^2 - (3*r1 + r2)^2 + (3*r1 - y2)^2

eq4 を解いて b を求める。

ans_b = solve(eq4, b)[2]
ans_b |> println

   -16*sqrt(3)*r1/7

eq3, eq5 の b に ans_b を代入する。

eq3 = 49eq3(b => ans_b)/r1^2 |> simplify
eq5 = eq5(b => ans_b) |> expand |> simplify;

eq1 |> println
eq3 |> println
eq5 |> println

   81*r1^2 - 18*sqrt(3)*r1*x2 - 18*r1*y2 - 4*r2^2 + 3*x2^2 + 2*sqrt(3)*x2*y2 + y2^2
   6075*r1^2 - 1350*sqrt(3)*r1*x2 - 9990*r1*y2 - 4332*r2^2 + 225*x2^2 + 1110*sqrt(3)*x2*y2 + 4107*y2^2
   -6*r1*r2 - 6*r1*y2 - r2^2 + x2^2 + y2^2

eq1, eq3, eq5 を解いて r1, x2, y2 を求める。

res2 = solve([eq1, eq3, eq5], (r1, x2, y2))[3];  # 3 of 3

res2[1] |> println
res2[2] |> println
res2[3] |> println

   2*r2*(sqrt(5) + 3)/9
   5*sqrt(3)*r2/6 + 2*sqrt(15)*r2/3
   3*r2/2

res2[1](r2 => 1/2).evalf() |> println
res2[2](r2 => 1/2).evalf() |> println
res2[3](r2 => 1/2).evalf() |> println

   0.581785330833310
   2.01268228522284
   0.750000000000000

大円の半径 r1 は,小円の半径 r2 の 2(√5 + 3)/9 倍である。
小円の直径が 1 寸のとき大円の直径は 2(√5 + 3)/9 = 1.16357066166662 寸である。
「術」の「小円径×(√20 + 6)/9」 と一致した。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   r1 = r2*2(√5 + 3)/9
   x2 = r2*(5√3/6 + 2√15/3)
   y2 = r2*3/2
   b = r1*(-16√3/7)
   a = r1*3√3
   x1 = -2√3*r1
   r0 = a/√3
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:blue, lw=0.5)
   circle(0, √3a/3, √3a/3)
   circle(x1, r1, r1, :green)
   circle(x2, y2, r2, :magenta)
   segment(a, 0, b, √3(b + 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, r0, "全円:r0\n(0,r0)", :red, :center, :bottom, delta=delta/2)
       point(x1, r1, "大円:r1\n(x1,r1)", :green, :center, :bottom, delta=delta/2)
       point(x2, y2, "小円:r2\n(x2,y2)", :magenta, :center, :bottom, delta=delta/2)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(0, √3a, " √3a", :blue, :left, :vcenter)
   end
end;

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

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

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