算額（その839）

算額（その839）

四十三　岩手県一関市真滝 熊野白山滝神社 弘化3年(1846)

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

大円内に水平な弦を引き，その上に中円2個，小円2個を乗せ，下に中円 2 個を上下に並べる。小円の直径を知って大円の直径を得る方法を述べよ。

問題文に「大円内に図のように圭（二等辺三角形）を設け」とあるのだが，図に圭はない。そのためかどうか，得られる解は，答えとも術とも違うものになる。

大円の半径と中心座標を R, (0, 0)
中円の半径と中心座標を r1, (0, r1 - R), (0, 3r1 - R), (r1, 5r1 - R)
小円の半径と中心座標を r2, (x2, 4r1 - R + r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms R::positive, r1::positive,
     r2::positive, x2::positive;
eq1 = r1^2 + (5r1 - R)^2 - (R - r1)^2
eq2 = x2^2 + (4r1 - R + r2)^2 - (R - r2)^2
eq3 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
solve([eq1, eq2, eq3], (R, r1, x2))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (225*r2/32, 9*r2/4, 21*r2/4)

大円の直径は小円の直径の 225/32 = 7.03125 倍である。

術では「大円径 = (√7 + 4)*小円径/0.64 = 10.38*小円径」としている。違いの原因がわからない。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (R, r1, x2) = r2 .* (225/32, 9/4, 21/4)
   @printf("R = %g;  r1 = %g;  x2 = %g\n", R, r1, x2)
   plot()
   circle(0, 0, R, :green)
   circle2(r1, 5r1 - R, r1)
   circle2(x2, 4r1 - R + r2, r2, :blue)
   circle(0, r1 - R, r1)
   circle(0, 3r1 - R, r1)
   y = 4r1 - R
   x = sqrt(R^2 - y^2)
   segment(-x, y, x, y, :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(0, r1 - R, "中円:r1\n(0,r1-R)", :red, :center, delta=-delta)
       point(0, 3r1 - R, "中円:r1\n(0,3r1-R)", :red, :center, delta=-delta)
       point(r1, 5r1 - R, "中円:r1\n(r1,5r1-R)", :red, :center, delta=-delta)
       point(x2, 4r1 - R + r2, "小円:r2,(x2,4r1-R+r2)", :blue, :center, delta=-delta)
       point(0, 4r1 - R, "4r1-R", :magenta, :center, :bottom, delta=delta)
       point(0, R, "R", :green, :center, :bottom, delta=delta)
       point(R, 0, " R", :green, :left, :bottom, delta=delta)
   end
end;

 

算額（その838）

算額（その838）

四十三　岩手県一関市真滝 熊野白山滝神社 弘化3年(1846)

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

長方形内に交差する四分円，甲円，乙円，等円を入れる。乙円の直径が与えられたとき，甲円の直径を求める方法を述べよ。

長方形の長辺，短辺を 2a, b
甲円の半径と中心座標を r1, (0, b - r1)
乙円の半径と中心座標を r2, (0, r2)
等円の半径と中心座標を r3, (0, 2r2 + r3), (x3, b - r3)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy

@syms a::positive, b::positive, r1::positive,
     r2::positive, r3::positive, x3::positive;
x3 = 2sqrt(r1*r3)
a = 2sqrt(b*r1)
eq1 = a^2 + r2^2 - (b - r2)^2 |> expand
eq2 = a^2 + (2r2 + r3)^2 - (b - r3)^2 |> expand
eq3 = (a - x3)^2 + (b - r3)^2 - (b + r3)^2 |> expand;

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")

   -b^2 + 4*b*r1 + 2*b*r2,  # eq1
   -b^2 + 4*b*r1 + 2*b*r3 + 4*r2^2 + 4*r2*r3,  # eq2
   -8*sqrt(b)*r1*sqrt(r3) + 4*b*r1 - 4*b*r3 + 4*r1*r3,  # eq3

たかが，三元連立方程式であるが SymPy の能力では，この連立方程式を解くことができない。
そこで，変数をへらしてゆく。

# solve([eq1, eq2, eq3], (b, r1, r3))

eq1 を b について解く。

res1 = solve(eq1, b)[1]
res1 |> println

   4*r1 + 2*r2

eq2, eq3 の b に res1 = 4r1 + 2r2 を代入する。

eq12 = eq2(b => res1)/8 |> simplify
eq13 = eq3(b => res1)/4 |> simplify;

eq12 |> println
eq13 |> println

   -r1*r2 + r1*r3 + r2*r3
   -2*sqrt(2)*r1*sqrt(r3)*sqrt(2*r1 + r2) + r1*r3 + 2*r1*(2*r1 + r2) - 2*r3*(2*r1 + r2)

eq12 を r3 について解く。

res2 = solve(eq12, r3)[1]
res2 |> println

   r1*r2/(r1 + r2)

eq13 を二乗し，ルートを外す。

eq14 = (eq13.args[1] + eq13.args[2] + eq13.args[3])^2 - (eq13.args[4])^2  |> expand |> simplify
eq14 |> println

   16*r1^4 + 16*r1^3*r2 - 40*r1^3*r3 + 4*r1^2*r2^2 - 36*r1^2*r2*r3 + 9*r1^2*r3^2 - 8*r1*r2^2*r3 + 12*r1*r2*r3^2 + 4*r2^2*r3^2

eq14 の r3 に res2 を代入する。

eq15 = eq14(r3 => res2) |> simplify
eq15 |> println

   r1^3*(16*r1^3 + 8*r1^2*r2 - 15*r1*r2^2 - 8*r2^3)/(r1^2 + 2*r1*r2 + r2^2)

eq15 に (r1^2 + 2*r1*r2 + r2^2)/r1^3 を掛けて単純にする。

eq16 = numerator(eq15)/r1^3 |> simplify
eq16 |> println

   16*r1^3 + 8*r1^2*r2 - 15*r1*r2^2 - 8*r2^3

特定の r2 のときの r1 は，eq16 を r1 について解けばよい。

eq17 = eq16(r2 => 1)
eq17 |> println

   16*r1^3 + 8*r1^2 - 15*r1 - 8

res3 = solve(eq17, r1)

   3-element Vector{Sym{PyCall.PyObject}}:
    -1/6 + (-1/2 - sqrt(3)*I/2)*(289/1728 + sqrt(237)*I/144)^(1/3) + 49/(144*(-1/2 - sqrt(3)*I/2)*(289/1728 + sqrt(237)*I/144)^(1/3))
    -1/6 + 49/(144*(-1/2 + sqrt(3)*I/2)*(289/1728 + sqrt(237)*I/144)^(1/3)) + (-1/2 + sqrt(3)*I/2)*(289/1728 + sqrt(237)*I/144)^(1/3)
                                              -1/6 + 49/(144*(289/1728 + sqrt(237)*I/144)^(1/3)) + (289/1728 + sqrt(237)*I/144)^(1/3)

SymPy では虚数解になってしまうが，evalf() で数値を見ると，虚部は極めて小さく実質的に実数解であることがわかる。しかも正の値になる解は r1 = 0.979095387936243 である。

res3[1].evalf() |> println
res3[2].evalf() |> println
res3[3].evalf() |> println
abs(res3[3].evalf()) |> println

   -0.549139398436065 + 0.e-20*I
   -0.929955989500177 + 0.e-23*I
   0.979095387936243 - 0.e-21*I
   0.979095387936243

pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(eq17, xlims=(-1, 1), ylims=(-15, 5))
hline!([0])

他のパラメータは，代入によって求めることができる。

r1 = 0.979095;  r3 = 0.494719;  b = 5.91638;  x3 = 1.39194;  a = 4.81361

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1
   r1 = 0.979095387936243
   r3 = r1*r2/(r1 + r2)
   b = 4*r1 + 2*r2
   x3 = 2sqrt(r1*r3)
   a = 2sqrt(b*r1)
   @printf("r1 = %g;  r3 = %g;  b = %g;  x3 = %g;  a = %g\n", r1, r3, b, x3, a)
   plot([a, a, -a, -a, a], [0, b, b, 0, 0], color=:black, lw=0.5)
   circle(a, 0, b, :green, beginangle=90, endangle=180)
   circle(-a, 0, b, :green, beginangle=0, endangle=90)
   circle(0, b - r1, r1)
   circle(0, r2, r2, :magenta)
   circle2(x3, b - r3, r3, :orange)
   circle(0, 2r2 + r3, r3, :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", :black, :left, :bottom, delta=delta)
       point(0, b, "b", :black, :center, :bottom, delta=delta)
       point(0, b - r1, "甲円:r1\n(0,b-r1)", :red, :center, delta=-delta)
       point(0, r2, "乙円:r2\n(0,r2)", :magenta, :center, delta=-delta)
       point(0, 2r2 + r3, " 等円:r3,(0,2r2+r3)", :black, :left, :vcenter)
       point(x3, b - r3, " 等円:r3,(x3,b-r3)", :black, :left, :vcenter)
   end
end;