算額あれこれ

算額問題をコンピュータで解きます

算額(その1074)

17 岩手県江刺市大通り 中善観音 文政10年(1827)

安富有恒:和算—岩手の現存算額のすべて,青磁社,東京都,1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード:円2個,二等辺三角形,正方形,斜線
#Julia #SymPy #算額 #和算 #数学


正方形内に二等辺三角形と甲円,乙円,界斜を容れる。乙円の直径が 1 寸のとき,界斜の長さはいかほどか。

注:界斜は乙円と甲円の両方に接している。

正方形の一辺の長さを \(2a\)
界斜と正方形の一辺との交点座標を \( (2a,\ b)\)
界斜を延長して x 軸との交点座標を \( (c,\ 0)\)
甲円の半径と中心座標を \(r_1,\ (a,\ r_1)\)
乙円の半径と中心座標を \(r_2,\ (x_2,\ 2a - r_2)\)
とおき,以下の連立方程式を解く。

円の中心から直線までの距離を求める際に,同じ直線でも定義に使う点の座標によって連立方程式が複雑になることがある。この算額でいえば界斜を \( (0,\ 2a) - (2a,\ b)\) で定義するのと \( (0,\ 2a) - (c,\ 0)\) で定義する場合に当たる。連立方程式が少しでも複雑になると SymPy で解けなくなってしまうので,工夫が必要になる。

図において ⊿ABX と ⊿OCX は相似で,界斜の長さは \(\text{BC}\cdot 2a/C\) になる。

include("julia-source.txt");  # julia-source.txt ソース

using SymPy

@syms a::positive, b::positive, c::positive,
     r1::positive, r2::positive, x2::positive
#eq1 = dist2(0, 2a, 2a, b, a, r1, r1)
#eq2 = dist2(0, 2a, 2a, b, x2, 2a - r2, r2)
eq1 = dist2(0, 2a, c, 0, a, r1, r1)/4a
eq2 = dist2(0, 2a, c, 0, x2, 2a - r2, r2)/4a
eq3 = dist2(0, 0, a, 2a, a, r1, r1)
eq4 = dist2(0, 0, a, 2a, x2, 2a - r2, r2)
eq5 = r1*x2 - r2*(c - a)
res = solve([eq2, eq3, eq4, eq5], (a, r1, x2, c))[2];

@syms d
res[1] |> factor |> println
apart(res[2]) |> factor |>  println
apart(res[3]) |> factor |> println
res[4] |> factor |> println

   r2*(-sqrt(10*sqrt(5) + 23) + 5 + 3*sqrt(5) + sqrt(5)*sqrt(10*sqrt(5) + 23))/4
   -r2*(-3*sqrt(10*sqrt(5) + 23) - 5 - sqrt(5) + sqrt(5)*sqrt(10*sqrt(5) + 23))/4
   r2*(-sqrt(10*sqrt(5) + 23) + sqrt(5) + 3 + sqrt(5)*sqrt(10*sqrt(5) + 23))/4
   r2*(7 + 4*sqrt(5) + sqrt(5)*sqrt(10*sqrt(5) + 23))/2

以上をまとめて,\(a,\ r_1,\ x_2,\ c\) は以下のようにして求めることができる。

r2 = 1/2
t = sqrt(10√5 + 23)
a = r2*(5 - t + √5(3 + t))/4
r1 = r2*(5 + 3t + √5(1 - t))/4
x2 = r2*(3 - t + √5(1 + t))/4
c = r2*(7 + √5(4 + t))/2
(r2, a, r1, x2, c)

   r2 = 0.5;  a = 2.50415;  r1 = 1.54765;  x2 = 1.69513;  c = 7.75107

乙円の直径が 1 寸のとき,界斜の長さは \( (2a/c)\sqrt{4a^2 +c^2} = 5.962810866213448\) 寸である。

sqrt(4a^2 +c^2)*(2a/c)

   5.962810866213448

術は以下の通りで,式の見かけは違うが同じ解である。

using SymPy
@syms 乙円径
天 = √Sym(5) + 3
(sqrt(70天 - 20) + 天)*乙円径/4

    \(\displaystyle \frac{\left( \sqrt{70 \cdot 天 - 20} + 天 \right) \cdot 乙円径}{4}\)

( (sqrt(70天 - 20) + 天)*乙円径/4).evalf() |>  println

    5.96281086621345*乙円径

描画関数プログラムのソースを見る

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 0.5
   t = sqrt(10√5 + 23)
   a = r2*(5 - t + √5(3 + t))/4
   r1 = r2*(5 + 3t + √5(1 - t))/4
   x2 = r2*(3 - t + √5(1 + t))/4
   c = r2*(7 + √5(4 + t))/2
   len = sqrt(4a^2 + c^2)*2a/c
   @printf("乙円の直径が %g のとき,界斜は %g である。\n", 2r2, len)
   @printf("r2 = %g;  a = %g;  r1 = %g;  x2 = %g;  c = %g\n", r2, a, r1, x2, c)
   plot([0, 2a, 2a, 0, 0], [0, 0, 2a, 2a, 0], color=:green, lw=0.5)
   plot!([0, a, 2a], [0, 2a, 0], color=:blue, lw=0.5)
   circle(a, r1, r1)
   circle(x2, 2a - r2, r2, :magenta)
   segment(0, 2a, c, 0)
   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(2a, 0, " 2a", :black, :left, :bottom, delta=delta/2, deltax=-0.5delta)
       point(2a, 2a*(c - 2a)/c, " (2a,b)", :black, :left, :vcenter)
       point(c, 0, " C", :red, :left, :bottom, delta=delta/2)
       point(a, 2a, " A", :red, :left, :bottom, delta=delta/2)
       point(0, 2a, "B ", :red, :right, :bottom, delta=delta/2)
       point(0, 0, "O ", :red, :right, :bottom, delta=delta/2)
       (x0, y0) = intersection(0, 2a, c, 0, 0, 0, a, 2a)
       point(float(x0), float(y0), "  X", :red, :left, :vcenter) 
       point(a, r1, "甲円:r1,(a,r1)", :red, :center, delta=-delta/2)
       point(x2, 2a - r2, "乙円:r2\n(x2,2a-r2)", :magenta, :right, delta=-4delta, deltax=-6delta)
       plot!(xlims=(-5delta, c + delta))
   end
end;


以下のアイコンをクリックして応援してください