算額（その287）

算額（その287）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(251)
長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

鉤股弦で，
股×弦 - 鉤×弦 = 425 平方寸
弦 - 鉤 = 18 寸
のとき，鉤，股，弦を求めよ。

以下の連立方程式を解く。
(鉤^2 + 股^2) - 弦^2 を使わないのは，虚数解を除くため。

include("julia-source.txt");

using SymPy
@syms 鉤::positive, 股::positive, 弦::positive;

eq1 = (股*弦 - 鉤*弦) - 425
eq2 = (弦 - 鉤) - 18
eq3 = (鉤^2 + 股^2) - 弦^2
eq3 = sqrt(鉤^2 + 股^2) - 弦
res = solve([eq1, eq2, eq3], (鉤, 股, 弦));

解は 2 組ある。

(鉤, 股, 弦) = res[1];
鉤.evalf(), 股.evalf(), 弦.evalf()

   (7.00000000000000, 24.0000000000000, 25.0000000000000)

(鉤, 股, 弦) = res[2];
鉤.evalf(), 股.evalf(), 弦.evalf()

   (25.3934598546653, 35.1875625010877, 43.3934598546653)

---

算額への招待 第4章 で術および現代解法が紹介されているが，方程式の係数に誤りがあるため，a=7, 25.3934598546653 の片方しか求めていない。

@syms a, b, c, A, B
A = 425
B = 18
eq1 = a^2 + b^2 - c^2
eq2 = b*c - a*c - A
eq3 = c - a - B;

solve(eq3, c)[1] |> println  # c

   a + 18

eq12 = eq2(c => a + 18)
eq12 |> println

   -a*(a + 18) + b*(a + 18) - 425

solve(eq12, b)[1] |> println  # b

   (a^2 + 18*a + 425)/(a + 18)

eq11 = eq1(c => a + 18, b => (a^2 + 18*a + 425)/(a + 18))
eq111 = eq11 |> simplify
eq111|> println

   (a^4 - 446*a^2 - 8028*a + 75649)/(a^2 + 36*a + 324)

factor(numerator(eq111)) |> println

   (a - 7)*(a^3 + 7*a^2 - 397*a - 10807)

pyplot(size=(200, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(eq111, xlims=(0, 30))
hline!([0])

res = solve(eq111)
res[1].evalf() |> println
res[2].evalf() |> println
res[3].evalf() |> println
res[4].evalf() |> println

   7.00000000000000
   -16.1967299273327 - 12.7768525871672*I
   -16.1967299273327 + 12.7768525871672*I
   25.3934598546653

件の現代解法では定数項が + 10807 となっており，もう一つの解が見逃されている。

A^2 - B^4 + (2A*B - 4B^3)a + (2A - 4B^2)a^2 + a^4 |> factor |> println

   (a - 7)*(a^3 + 7*a^2 - 397*a - 10807)

SymPy の solve() は，ここに展開したような式変形と求解を間違いなくやってくれている。

 

算額（その286）

算額（その286）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(252)
長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

等脚台形内に，甲，乙，丙の 5 個の円がある。それぞれ台形に内接したり，互いに外接しあっている。
甲円の直径が 1 寸のとき，乙円の直径を求めよ。

台形の上底と下底の長さを 2b, 2a，高さを c とする。
甲円の半径と中心座標 r1, (0, r1)
乙円の半径と中心座標 r2, (0, r2 + 2r1)
上の丙円の半径と中心座標 r3, (0, r3 + 2r1 + 2r2)
右の丙円の半径と中心座標 r3, (x3, r3)
c = 2(r1 + r2 + r3) である。

以下の連立方程式を解く。なお，円の中心から直線までの平方距離は分数式になるので分子 = 0 を解けば数式は若干短くなる。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive,
     r3::positive, x3::positive;

r1 = 1//2
c = 2(r1 + r2 + r3)
eq1 = x3^2 + (r3 - r1)^2 - (r3 + r1)^2
eq2 = x3^2 + (r2 + 2r1 - r3)^2 - (r2 + r3)^2
eq3 = distance(a, 0, b, c, 0, c - r3) - r3^2
eq3 = -2*a*b*r3 - 2*b^2*r2 - b^2 + 2*r2*r3^2 + 2*r3^3 + r3^2
eq4 = distance(a, 0, b, c, 0, r2 + 2r1) - r2^2
eq4 = 4*a^2*r2*r3 + 4*a^2*r3^2 + 4*a*b*r2^2 + 4*a*b*r2*r3 + 2*a*b*r2 + 4*a*b*r3 + 2*b^2*r2 + b^2 - 4*r2^4 - 8*r2^3*r3 - 4*r2^3 - 4*r2^2*r3^2 - 4*r2^2*r3 - r2^2
eq5 = distance(a, 0, b, c, x3, r3) - r3^2
eq5 = 2*a^2*r2 + a^2 + 2*a*b*r3 - 4*a*r2*x3 - 2*a*r3*x3 - 2*a*x3 - 2*b*r3*x3 - 2*r2*r3^2 + 2*r2*x3^2 - 2*r3^3 - r3^2 + 2*r3*x3^2 + x3^2;
# res = solve([eq1, eq2, eq3, eq4, eq5])#, (a, r1, r2, r3, x3, r6, x6, y6))

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

   x3^2 + (r3 - 1/2)^2 - (r3 + 1/2)^2,  # eq1
   x3^2 - (r2 + r3)^2 + (r2 - r3 + 1)^2,  # eq2
   -2*a*b*r3 - 2*b^2*r2 - b^2 + 2*r2*r3^2 + 2*r3^3 + r3^2,  # eq3
   4*a^2*r2*r3 + 4*a^2*r3^2 + 4*a*b*r2^2 + 4*a*b*r2*r3 + 2*a*b*r2 + 4*a*b*r3 + 2*b^2*r2 + b^2 - 4*r2^4 - 8*r2^3*r3 - 4*r2^3 - 4*r2^2*r3^2 - 4*r2^2*r3 - r2^2,  # eq4
   2*a^2*r2 + a^2 + 2*a*b*r3 - 4*a*r2*x3 - 2*a*r3*x3 - 2*a*x3 - 2*b*r3*x3 - 2*r2*r3^2 + 2*r2*x3^2 - 2*r3^3 - r3^2 + 2*r3*x3^2 + x3^2,  # eq5

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=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r2, r3, x3, a, b) = u
   return [
       x3^2 + (r3 - 1/2)^2 - (r3 + 1/2)^2,  # eq1
       x3^2 - (r2 + r3)^2 + (r2 - r3 + 1)^2,  # eq2
       -2*a*b*r3 - 2*b^2*r2 - b^2 + 2*r2*r3^2 + 2*r3^3 + r3^2,  # eq3
       4*a^2*r2*r3 + 4*a^2*r3^2 + 4*a*b*r2^2 + 4*a*b*r2*r3 + 2*a*b*r2 + 4*a*b*r3 + 2*b^2*r2 + b^2 - 4*r2^4 - 8*r2^3*r3 - 4*r2^3 - 4*r2^2*r3^2 - 4*r2^2*r3 - r2^2,  # eq4
       2*a^2*r2 + a^2 + 2*a*b*r3 - 4*a*r2*x3 - 2*a*r3*x3 - 2*a*x3 - 2*b*r3*x3 - 2*r2*r3^2 + 2*r2*x3^2 - 2*r3^3 - r3^2 + 2*r3*x3^2 + x3^2,  # eq5
      ]
end;

iniv = [1.3, 0.7, 1.2, 2.1, 0.5]
res = nls(H, ini=iniv);
println(res);

   ([1.3090169943749475, 0.6909830056250524, 1.1755705045849463, 2.1266270208801, 0.5020285397155682], true)

乙円の半径 = 1.3090169943749475 なので，直径はほぼ 2.618 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   (r2, r3, x3, a, b) = res[1]
   c = 2(r1 + r2 + r3)
   plot([a, b, -b, -a, a], [0,c, c, 0, 0], color=:black, lw=0.5)
   circle(0, r1, r1, :blue)
   circle(0, r2 + 2r1, r2, :orange)
   circle(0, c - r3, r3)
   circle(x3, r3, r3)
   circle(-x3, r3, r3)
   if more
       str = @sprintf("r2 ≒ %.5f", 2r2)
       point(0, r1, " r1\n 甲:r1", :blue)
       point(0, r2 + 2r1, " r2+2r1\n 乙:r2\n $str", :orange)
       point(0, c - r3, " r3+2r2+2r1\n 丙:r3", :red)
       point(x3, r3, " 丙:r3\n (x3,r3)", :red)
       point(b, c, "(b,2(r1+r2+r3))", :black, :left, :bottom)
       point(a, 0, "a ", :black, :right, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;