裏 RjpWiki

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

算額(その847)

2024年04月09日 | Julia

算額(その847)

百四 岩手県大船渡市田茂山 根城八幡宮 天保12年(1841)

山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市. http://www.wasan.jp/yamamura/yamamura.html

正方形内に左上隅に縦に 2 個,右上隅に 1 個の等円を入れ,1 つの頂点から 2 個の等円に接する斜線を引く。等円の直径が 4 寸のとき,斜線の長さを求めよ。

正方形の一辺の長さを a,上辺と斜線の交点座標を (b, a)
等円の半径と中心座標を r, (r, a - r), (a, a - 3r), (a - r, a - r)
斜線の長さを L
とおき,以下の連立方程式を解く。
左の三角形の面積を,3 辺を底辺とする三角形の高さ(等円の半径),右の台形の面積を,4 辺を底辺とする三角形の高さ(等円の半径)で面積を分解する条件式である。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r::positive, L::positive;
L = sqrt(a^2 +  b^2)
eq1 = (a + L)r + b*3r - a*b
eq2 = (L + (a - b) + a)r + a*(a - r) - (2a^2 - a*b) |> simplify;

たかだか2元連立方程式だが,r が定数のままだと SymPy の能力では解けないようだ。
r = 4//2
を加えると解ける。
整数方程式だろうと予測して,総当たりで解を求めるのもよかろう。

using SymPy
@syms a::positive, b::positive, r::positive, L::positive;
r = 4//2
L = sqrt(a^2 +  b^2)
eq1 = (a + L)r + b*3r - a*b
eq2 = (L + (a - b) + a)r + a*(a - r) - (2a^2 - a*b) |> simplify
res = solve([eq1, eq2], (a, b))[1]

   (12, 9)

a = 12, b = 9 とわかったので,斜線の長さは sqrt(a^2 +  b^2) = 15 である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 4//2
   (a, b) = (12, 9)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:blue, lw=0.5)
   circle(r, a - r, r)
   circle(r, a - 3r, r)
   circle(a - r, a - r, r)
   segment(0, 0, b, a, :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(r, a - r, "等円:r,(r,a-r)", :red, :center, delta=-delta/2)
       point(r, a - 3r, "等円:r,(r,a-3r)", :red, :center, delta=-delta/2)
       point(a - r, a - r, "等円:r,(a-r,a-r)", :red, :center, delta=-delta/2)
       point(b, a, "(b,a)", :magenta, :center, :bottom, delta=delta/2)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(0, a, " a", :blue, :left, :bottom, delta=delta/2)
       
   end
end;

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

算額(その846)

2024年04月09日 | Julia

算額(その846)

百四 岩手県大船渡市田茂山 根城八幡宮 天保12年(1841)

山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市. http://www.wasan.jp/yamamura/yamamura.html

天円の中に水平な弦を引き,その上に地円 3 個,人円 1個,下に地円 1 個を入れる。人円の直径が与えられたとき,地円の直径を求めよ。

人円の中心座標は原点になる。
天円の半径と中心座標を R, (0, 0)
地円の半径と中心座標を r1, (0, R - r1), (x1, r1 - r2), (0, r1 - R)
人円の半径と中心座標を r2, (0, 0)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, x1::positive, r2::positive;

eq1 = x1^2 + (r1 - r2)^2 - (R - r1)^2
eq2 = x1^2 + 4r2^2 - 4r1^2
eq3 = x1^2 + (r1 - r2)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3], (R, r1, x1))[1]

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

地円の半径 r1 は人円の半径 r2 の (1 + √5)/2 倍である。
例えば,人円の直径が 1 のとき,地円の直径は 1.618033988749895 である。
ちなみに天円の直径は 人円の直径の 2 + √5 倍の 4.23606797749979 である。

(1 + √5)/2, 2 + √5

   (1.618033988749895, 4.23606797749979)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (R, r1, x1) = (2*r2 + sqrt(5)*r2, r2*(1 + sqrt(5))/2, 2*r2*sqrt(1/2 + sqrt(5)/2))
   @printf("人円の直径 = %g;  地円の直径 = %g;  天円の直径 = %g\n", 2r2, 2r1, 2R)
   plot()
   circle(0, 0, R)
   circle2(x1, r1 - r2, r1, :blue)
   circle22(0, R - r1, r1, :blue)
   circle(0, 0, r2, :green)
   y = 2r1 - R
   x = sqrt(R^2 - y^2)
   segment(-x, y, x, y, :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(R, 0, " R", :red, :left, :bottom, delta=delta/2)
       point(0, R - r1, " 地円:r1,(0,R-r1)", :blue, :center, delta=-delta/2)
       point(0, r1 - R, " 地円:r1,(0,r1-R)", :blue, :center, delta=-delta/2)
       point(x1, r1 - r2, " 地円:r1,(x1,r1-r2)", :blue, :center, delta=-delta/2)
       point(0, 0, "人円:r2,(0,0)", :green, :center, delta=-delta/2)
       point(0, -r2, "-r2", :black, :center, delta=-delta/2)
   end
end;

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

算額(その845)

2024年04月09日 | Julia

算額(その845)

九十九 雨宝堂 現在岩手県江刺市中善観音堂 文政10年(1827)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

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

岩手県大船渡市 稲子沢雨宝堂 文政10年(1827)
一関市博物館>>和算に挑戦>>平成24年度出題問題&回答例

https://www.city.ichinoseki.iwate.jp/museum/wasan/h24/index.html

正方形内に四分円,甲円,乙円,丙円が入っている。丙円の直径が 1 寸のとき,乙円の直径を求めよ。

四分円の半径(正方形の一辺の長さ)と中心座標を r0, (0, 0), (r0, r0)
甲円の半径と中心座標を r1, (x1, r0 - x1)
乙円の半径と中心座標を r2, (r2, r2)
丙円の半径と中心座標を r3, (x3, x3)
とおき,以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms r0::positive, r1::positive, x1::positive, r2::positive,
     r3::positive, x3::positive;
eq1 = x1^2 + (r0 - x1)^2 - (r0 - r1)^2
eq2 = 2(r0 - r2)^2 - (r0 + r2)^2
eq3 = 2(r0 - x3)^2 - (r0 - r3)^2
eq4 = (x3 - x1)^2 + (r0 - x1 - x3)^2 - (r1 + r3)^2
eq5 = 2(r0 - 2x1)^2 - 4r1^2
res = solve([eq1, eq2, eq3, eq4, eq5], (r0, r1, x1, r2, x3))
res[5]

   (r3*(4*sqrt(2) + 7), r3*(4*sqrt(2) + 7)/4, -sqrt(2)*r3*sqrt(56*sqrt(2) + 81)/8 + r3*(4*sqrt(2) + 7)/2, -2*sqrt(2)*r3*sqrt(56*sqrt(2) + 81) + 3*r3*(4*sqrt(2) + 7), r3*(sqrt(2) + 3))

8 組の解が得られるが,5 番目のものが適解である。
r0, r1, x1, r2, x3 は順に 丙円の半径 r3 の 4√2 + 7 倍, √2 + 7/4 倍, 2√2 + 7/2 - √2sqrt(56√2 + 81)/8 倍, 12√2 + 21 - 2√2sqrt(56√2 + 81) 倍, √2 + 3 倍である。

簡約化すると,乙円の直径は丙円の直径の 5 - 2√2 倍である。
丙円の直径が 1 寸のとき,乙円の直径は 2.1715728752538084 = 2 寸 1 分 7 厘 1 毛 5 糸有奇である。

# r2
res[5][4] |> sympy.sqrtdenest |> simplify |> println

   r3*(5 - 2*sqrt(2))

5 - 2√2

   2.1715728752538097

# x1
res[5][3] |> sympy.sqrtdenest |> simplify |> println

   r3*(9*sqrt(2) + 20)/8

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

r0 = 6.32843;  r1 = 1.58211;  x1 = 2.0455;  r2 = 1.08579;  x3 = 2.20711

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 1/2
   (r0, r1, x1, r2, x3) = r3 .* (4√2 + 7, √2 + 7/4, (9√2 + 20)/8, 5 - 2√2, √2 + 3)
   @printf("乙円の直径 = %g\n", 2r2)
   @printf("r0 = %g;  r1 = %g;  x1 = %g;  r2 = %g;  x3 = %g\n", r0, r1, x1, r2, x3)
   plot([0, r0, r0, 0, 0], [0, 0, r0, r0, 0], color=:black, lw=0.5)
   circle(0, 0, r0, :blue, beginangle=0, endangle=90)
   circle(r0, r0, r0, :blue, beginangle=180, endangle=270)
   circle(r2, r2, r2)
   circle(r0 - r2, r0 - r2, r2)
   circle(x1, r0 - x1, r1, :magenta)
   circle(r0 - x1, x1, r1, :magenta)
   circle(x3, x3, r3, :green)
   circle(r0 - x3, r0 - x3, 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(x1, r0 - x1, "甲円:r1,(x1,r0-x1)", :magenta, :center, delta=-delta/2)
       point(r2, r2, "乙円:r2,(r2,r2)", :red, :center, delta=-delta/2)
       point(x3, x3, " 丙円:r3,(x3,x3)", :green, :left, :vcenter)
       point(0, r0, " r0", :black, :left, :bottom, delta=delta/2)
       point(r0, 0, " r0", :black, :left, :bottom, delta=delta/2)
   end
end;

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

算額(その844)

2024年04月09日 | Julia

算額(その844)

四十四 岩手県一関市真滝 熊野白山滝神社 文久元年(1861)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

一関市博物館>>和算に挑戦>>平成20年度出題問題&回答例
https://www.city.ichinoseki.iwate.jp/museum/wasan/h20/index.html

圭(二等辺三角形)の中に斜線を入る。(1)斜線の長さと(2)底辺と,(3)斜辺と斜線の交点と頂点までの距離の 3 つの長さが 1 寸である。上部の三角形に内接する円を描いたとき図に色付した部分の面積(黒積)を求めよ。
注:図では「圭」と言いながら,不等辺三角形に見える図が描かれている。また,どこが黒積なのかという記述もない。

まず,θ を求める。
AD = BD = BC = a,∠BAC = θ とする。
△ABC は △BCD と相似なので,∠ACB = ∠ABC = ∠ABD + ∠CBD = 2θ, ∠ACB + ∠ABC + ∠BAC = 5θ = 180° ゆえ,θ = 36° である。
次に,円の半径 r を求める。

include("julia-source.txt")

using SymPy

@syms L::positive, a::positive, θ::positive, r::positive, h::positive;
θ = Sym(36)
h = a*sind(θ)
L = sqrt(a^2 - h^2)
eq4 = tand(θ/2) - r/L;
r = solve(eq4, r)[1];
println("r = $r")

   r = a*sqrt(50 - 10*sqrt(5))/20

灰色の部分の面積 S は,△BEF の面積から,半径 r の円の面積の 72/360 を引いたものものなので,以下のようになる。

S = L*r/2 - PI*r^2*((90 - θ/2)/360)
S |> println

   -pi*a^2*(50 - 10*sqrt(5))/2000 + a*sqrt(50 - 10*sqrt(5))*sqrt(-a^2*(5/8 - sqrt(5)/8) + a^2)/40

求める面積(黒積)は S の 2 倍なので,a = 1 を代入して 0.125831216718928 を得る。

算額では,「答曰黒積▢分二厘毛八糸」とあるが,この問題でも誤記(?)がある。「黒積一分二厘毛八四」であろう。

2*S(a => 1).evalf() |> println

   0.125831216718928

なお,tan(18°),tan(36°) は結構きれいな形で表すことができるので,覚えておいて損はない。

tand(Sym(18)) |> println
tand(Sym(36)) |> println

   sqrt(1 - 2*sqrt(5)/5)
   sqrt(5 - 2*sqrt(5))

function draw(a, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, θ) = (1, 36)
   r = a*sqrt(50 - 10√5)/20
   h = (a/2)tand(2θ)
   plot([a, 0, -a, a] ./ 2, [0, h, 0, 0], color=:black, lw=0.5)
   segment(a/2, 0, a/2 - a*cosd(θ), a*sind(θ))
   h2 = a*((12 + 4√5)sqrt(25 - 5√5) + 5(√2 + √10)*sqrt(10 + 2√5))/(20(3 + √5)^1.5)
   (Ex, Ey) = (a/4, h/2)
   plot!([Ex, 0, a/2, Ex], [Ey, h2, 0, Ey], color=:black, lw=0.5, seriestype=:shape, fillcolor=:gray80)
   circlef(0, h2, r, :white)
   circle(0, h2, r)
   plot!([Ex, 0, a/2, Ex], [Ey, h2, 0, Ey], color=:black, lw=0.5)
   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, h, " A:(0,h)", :black, :left, :vcenter)
       point(a/2, 0, " B:(a/2,0)", :black, :left, :vcenter)
       point(-a/2, 0, " C:(-a/2,0)", :black, :left, :bottom, delta=delta)
       point(a/2 - a*cosd(θ), a*sind(θ), "D ", :black, :right, :vcenter)
       point(Ex, Ey, " E", :black, :left, :vcenter)
       point(0, h2, "F:(0,h2) ", :red, :right, :vcenter)
       
   end
end;

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

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

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