裏 RjpWiki

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

算額(その1079)

2024年06月19日 | Julia

算額(その1079)

九十七 岩手県大船渡市猪川町 雨宝堂 現雨宝山竜宝院 文政7年(1824)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

キーワード:円9個,外円,弦2本

全円の中に 2 本の水平な弦を引き,大円 1 個,甲円 5 個,乙円 2 個を容れる。乙円の直径が 1 寸のとき,全円の直径はいかほどか。

全円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r0, (0, R - r0)
甲円の半径と中心座標を r1, (0, R - r1), (r1, R - 3r1), (r1, R - 2r0 - r1)
乙円の半径と中心座標を r2, (x2, R - 2r0 + r2)
とおき,以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms R::positive, r0::positive,
     r1::positive, r2::positive, x2::positive
eq1 = x2^2 + (R - 2r0 + r2)^2 - (R - r2)^2
eq2 = r1^2 + (R - 2r0 - r1)^2 - (R - r1)^2
eq3 = r1^2 + (3r1 - r0)^2 - (r0 - r1)^2
eq4 = x2^2 + (r0 - r2)^2 - (r0 + r2)^2
res = solve([eq1, eq2, eq3, eq4], (R, r1, r0, x2))[1]

   (14641*r2/3240, 121*r2/90, 121*r2/40, 11*sqrt(10)*r2/10)

全円の半径 R は 乙円の半径 r2 の 14641/3240 = 4.518827160493827 倍である。
術では「4 と 1681/3240」と小学校以来使ったことのない帯分数で表している。
したがって,乙円の直径が 1 寸のとき,全円の直径は 14641/3240 = 4.518827160493827 寸である。

その他のパラメータは以下のとおりである。
r2 = 0.672222;  R = 2.25941;  r1 = 0.672222;  r0 = 1.5125;  x2 = 1.73925

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (R, r1, r0, x2) = r2 .* (14641/3240, 121/90, 121/40, 11√10/10)
   @printf("乙円の直径が %g のとき,全円の直径は %g である。\n", 2r2, 2R)
   @printf("r2 = %g;  R = %g;  r1 = %g;  r0 = %g;  x2 = %g\n", r1, R, r1, r0, x2)
   plot()
   circle(0, 0, R, :green)
   circle(0, R - r0, r0, :magenta)
   circle(0, R - r1, r1)
   y01 = R - 2r1
   x01 = sqrt(R^2 - y01^2)
   segment(-x01, y01, x01, y01)
   circle2(r1, R - 3r1, r1)
   y02 = R - 2r0
   x02 = sqrt(R^2 - y02^2)
   segment(-x02, y02, x02, y02)
   circle2(r1, y02 - r1, r1)
   circle2(x2, y02 + r2, r2, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, R, " R", :green, :left, :bottom, delta=delta/2)
       point(0, R - r1, "甲円:r1,(0,R-r1)", :red, :center, delta=-delta/2)
       point(r1, R - 3r1, "甲円:r1\n(r1,R-3r1)", :red, :center, delta=-delta/2)
       point(r1, y02 - r1, "甲円:r1\n(r1,y02-r1)", :red, :center, delta=-delta/2)
       point(x2, y02 + r2, "乙円:r2\n(x2,y02+r2)", :blue, :center, :bottom, delta=delta/2)
       point(0, y01, "y01", :black, :center, delta=-delta/2)
       point(0, y02, "y02", :black, :center, delta=-delta/2)
   end
end;

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

算額(その1078)

2024年06月19日 | Julia

算額(その1078)

百 大船渡市猪川町 雨宝堂 現雨宝山竜宝院 文政10年(1827)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html
キーワード:円2個,正三角形,正方形

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

正方形の一辺の長さを a
正三角形の頂点の座標を (0, b), (b, 0), (a/2, √3a/2)
日円の半径と中心座標を r1, (x1, y1); y1 = x1
月円の半径と中心座標を r2, (r2, y2)
とおき,以下の連立方程式を解く。

日円と月円の方程式は独立なので,別々に解く。

1. 日円について

include("julia-source.txt")

using SymPy

@syms a::positive, b::positive,
     r1::positive, x1::positive, y1::positive,
     r2::positive, y2::positive
b = a*(√Sym(3) - 1)  # a - 2(a - √Sym(3)a/2)
y1 = x1
eq1 = dist2(0, b, a, a, x1, y1, r1)
eq2 = dist2(a/2, √Sym(3)a/2, a, 0, x1, y1, r1)
res1 = solve([eq1, eq2], (r1, x1))[2];
r1 = apart(res1[1]) |> simplify
r1 |> println
x1 = apart(res1[2]) |> simplify
x1 |> println

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

日円の半径は正方形の一辺の長さの (√3 - 1)*sqrt(3 - √6 + √3 - √2)/4 倍である。

2. 月円について

r2 = y2*(2 - √Sym(3))
eq3 = dist2(0, b, b, 0, r2, y2, r2)
y2 = solve(eq3, y2)[2]
y2 = y2 |> sympy.sqrtdenest |> simplify
r2 = y2*(2 - √Sym(3)) |> simplify
r2 |> println
y2 |> println

   a*(-7*sqrt(2) - 5*sqrt(3) + 9 + 4*sqrt(6))/2
   a*(-2*sqrt(2) - sqrt(3) + sqrt(6) + 3)/2

月円の半径は正方形の一辺の長さの (-7√2 - 5√3 + 9 + 4√6)/2 倍である。

3. 月円と日円の比について

月円の直径が 1 のときの日円の直径は r2/r1 = 1.43185165257814 である。

ついでに,月円の直径が 1 寸になるときの正方形の一辺の長さを求めてみる。

@syms d
日円の直径 = apart(r1/r2, d) |> simplify
日円の直径 |> println
日円の直径.evalf() |> println

   -1/2 + sqrt(2)/2 + sqrt(6)/2
   1.43185165257814

正方形(小正三角形)の一辺 の長さは r2 = 1/2 を解いて,a = √2/4 + √6/4 + 3/2 + √3 = 4.197976633857945 である。

eq = r2 - 1//2
res = solve(eq, a)[1]
ans_a = apart(res[1], d)
ans_a |> println
ans_a.evalf() |> println

   sqrt(2)/4 + sqrt(6)/4 + 3/2 + sqrt(3)
   4.19797663385795

4. 図を描いて確認する

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = sqrt(2)/4 + sqrt(6)/4 + 3/2 + sqrt(3)
   b = a*(√3 - 1)
   r1 = a*(√3 - 1)*sqrt(3 - √6 + √3 - √2)/4
   x1 = a*(-√2 + 2 + √6)/4
   r2 = a*(-7√2 - 5√3 + 9 + 4√6)/2
   y2 = a*(-2√2 - √3 + √6 + 3)/2
   @printf("月円の直径が %g のとき,日円の直径は %g である。\n", 2r2, 2r1)
   @printf("a = %g;  r1 = %g;  x1 = %g;  r2 = %g;  y2 = %g\n", a, r1, x1, r2, y2)
   plot([0, a, a, 0, 0], [0, 0, a, a], color=:blue, lw=0.5)
   plot!([0, a, a/2, 0], [0, 0, √3a/2, 0], color=:green, lw=0.5)
   plot!([b, a, 0, b], [0, a, b, 0], color=:magenta, lw=0.5)
   circle(x1, x1, r1)
   circle(r2, y2, r2, :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(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(0, a, "a ", :blue, :right, :vcenter)
       point(x1, x1, "日円:r1,(x1,x1)", :red, :center, delta=-delta/2)
       point(r2, y2, "月円:r2,(r2,y2)", :orange, :center, delta=-delta/2)
       point(0, b, "b ", :blue, :right, :vcenter)
       point(b, 0, "b", :blue, :center, delta=-delta)
       point(a/2, √3a/2, "(a/2,√3a/2)", :green, :center, :bottom, delta=delta, deltax=-3delta)
       plot!(xlims=(-5delta, a + 5delta), ylims=(-5delta, a + 5delta))
   end
end;

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

算額(その1077)

2024年06月19日 | Julia

算額(その1077)

百 大船渡市猪川町 雨宝堂 現雨宝山竜宝院 文政10年(1827)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html
キーワード:円8個,累円,直線上

直線上に甲円を中心として,左右に乙円,丙円,丁円 が互いに外接している。その上にそれら全ての円に外接する天円が載っている。乙円,丁円の直径がそれぞれ 4 寸,16 寸のとき,天円の直径はいかほどか。

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

include("julia-source.txt")

using SymPy

@syms R::positive, r1::positive, x1, r2::positive, x2::positive,
     r3::positive, x3::positive, r4::positive, x4::positive
(r2, r4) = (4, 16) .// 2
x1 = 0
eq1 = x2^2 + (2r1 + R - r2)^2 - (R + r2)^2 |> expand
eq2 = x3^2 + (2r1 + R - r3)^2 - (R + r3)^2 |> expand
eq3 = x4^2 + (2r1 + R - r4)^2 - (R + r4)^2 |> expand
eq4 = (x2 - x1)^2 + (r2 - r1)^2 - (r2 + r1)^2 |> expand
eq5 = (x3 - x2)^2 + (r3 - r2)^2 - (r3 + r2)^2 |> expand
eq6 = (x4 - x3)^2 + (r4 - r3)^2 - (r4 + r3)^2 |> expand
eq4 = x2^2 - 4r1*r2
eq5 = (x3 - x2)^2 - 4r2*r3
eq6 = (x4 - x3)^2 - 4r3*r4
res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (R, r1, x2, r3, x3, x4))[1]

   (49/4, 7/4, sqrt(14), 28/9, 7*sqrt(14)/3, 5*sqrt(14))

乙円,丁円の直径がそれぞれ 4 寸,16 寸のとき,天円の直径は 49/2 = 24.5 寸である。

「答」では,32 寸としているが,それでは図が描けない。

条件として与えられたのが「乙円径 4 寸,丁円径 16 寸」,答えは「天円径 32 寸」
これを図にしようと,まず直径 32 の天円を描き,それに接するように直径 4 の甲円を描くが,その左に天円に外接する甲円が描けるように天円と乙円の中心座標を調整する。

次に,直径 16 の丁円を天円に外接するように描く。これはたいした苦も無く描ける。次いで,その左に天円と外接するように丙円を描く。これも容易に描ける。
が,しかし。乙円と丙円は外接しない。
要するに,「問」と「答」を満たすような解はないということだ。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r4) = [4, 16] ./ 2
   (R, r1, x2, r3, x3, x4) = (49/4, 7/4, sqrt(14), 28/9, 7*sqrt(14)/3, 5*sqrt(14))
   @printf("乙円,丁円の直径が %g, %g のとき,天円の直径は %g である。\n", 2r2, 2r4, 2R)
   @printf("(r1 => %.15g, r2 => %.15g, r3 = %.15g, r4 = %.15g)\n", r1, r2, r3, r4)
   plot()
   circle(0, 2r1 + R, R, :blue)
   circle(0, r1, r1)
   circle2(x2, r2, r2, :green)
   circle2(x3, r3, r3, :magenta)
   circle2(x4, r4, 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, 2r1 + R, "天円:R,(0,2r1+R)", :blue, :center, delta=-delta)
       point(x4, r4, " 丁円:r4,(x4,r4)", :orange, :center, delta=-delta)
       point(x3, r3, " 丙円:r3,(x3,r3)", :magenta, :center, delta=-delta, deltax=10delta)
       point(x2, r2, "乙円:r2,(x2,r2)", :green, :center, delta=10delta)
       point(0, r1, "甲円:r1,(0,r1) ", :black, :right, :vcenter)
   end
end;

インチキな図を書くプログラム

function draw2(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r4) = [4, 16] ./ 2
   (R, r1, x2, r3, x3, x4) = (16, 1.78, 3.73, 3.4, 10.6,20.96)    
   @printf("乙円,丁円の直径が %g, %g のとき,天円の直径は %g である。\n", 2r2, 2r4, 2R)
   @printf("(r1 => %.15g, r2 => %.15g, r3 = %.15g, r4 = %.15g)\n", r1, r2, r3, r4)
   plot()
   circle(0, 2r1 + R, R, :blue)
   circle(0, r1, r1)
   circle2(x2, r2, r2, :green)
   circle2(x3, r3, r3, :magenta)
   circle2(x4, r4, 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, 2r1+R, " 天円径=32", :blue, :center, delta=-delta, mark=false)
       point(x4, r4, " 丁円径=16", :orange, :center, delta=-delta, mark=false)
       point(x3, r3, " 丙円", :magenta, :center, :vcenter, mark=false)
       point(x2, r2, "乙円径=4", :green, :center, delta=10delta, mark=false)
       point(0, r1, "甲円", :red, :center, :vcenter, mark=false)
   end
end;

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

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

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