算額（その47）

算額（その47）

岩手県東磐井郡藤沢町 藤勢寺 元治 2 年
http://www.wasan.jp/iwate/fujise.html

七十九 藤沢町藤沢道場 藤勢寺 元治2年
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

前者のリンクは，画像が不鮮明で，問，答，術とも正確に読めていなかった。
後者のリンク内容に基づき，第二版を作成する。

1 辺が 1 寸の正方形に2種類の円が収まっている。それぞれの径を求めよ。
注：算額の図では円が2種類あるように見えるが，問では「小円6個」と書いてある。

正方形の一辺の長さを a
小円の半径と中心座標を r, (x1, r), (r, x1), (x2, x2)
菱形の頂点座標を (x0, x0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, r::positive,
     x0::positive, x1::positive, x2::positive;
eq1 = dist2(a, 0, 0, a, x2, x2, r)
eq2 = dist2(a, 0, x0, x0, x2, x2, r)
eq3 = dist2(a, 0, x0, x0, x1, r, r)
eq4 = dist2(0, 0, a, a, x1, r, r)
res = solve([eq1, eq2, eq3, eq4], (r, x0, x1, x2))[6]  # 6 of 8

   (a*(-4553*sqrt(-10 + 12*sqrt(2)) - 4605*sqrt(-5 + 6*sqrt(2)) - 3537 + 17080*sqrt(2))/(2*(-4663*sqrt(-5 + 6*sqrt(2)) - 351*sqrt(2)*sqrt(-5 + 6*sqrt(2)) + 3845 + 4092*sqrt(2))), a*(-104*sqrt(2) - 3 + 25*sqrt(-10 + 12*sqrt(2)) + 45*sqrt(-5 + 6*sqrt(2)))/(12*(-6*sqrt(2) + 1 + 4*sqrt(-5 + 6*sqrt(2)))), (-4*sqrt(2)*a*sqrt(-5 + 6*sqrt(2)) - 5*a*sqrt(-5 + 6*sqrt(2)) + 5*sqrt(2)*a + 13*a)/(8 - 4*sqrt(-5 + 6*sqrt(2))), a*(-sqrt(2)/8 + 1/4 + sqrt(-5/32 + 3*sqrt(2)/16)))

a = 10
r = a*(-4553*sqrt(-10 + 12√2) - 4605*sqrt(-5 + 6√2) - 3537 + 17080√2)/(2*(-4663*sqrt(-5 + 6√2) - 351√2*sqrt(-5 + 6√2) + 3845 + 4092√2))
x0 = a*(-104√2 - 3 + 25*sqrt(-10 + 12√2) + 45*sqrt(-5 + 6√2))/(12*(-6√2 + 1 + 4*sqrt(-5 + 6√2)))
x1 = (-4√2a*sqrt(-5 + 6√2) - 5*a*sqrt(-5 + 6√2) + 5√2a + 13*a)/(8 - 4*sqrt(-5 + 6√2))
x2 = a*(-sqrt(2)/8 + 1/4 + sqrt(-5/32 + 3√2/16))
(r, x0, x1, x2)

   正方形の一辺の長さが 10 のとき，小円の直径は 2.73661 である。
   a = 10;  r = 1.36831;  x0 = 2.98964;  x1 = 3.30338;  x2 = 4.03246

r は，a の 6次式を解いて求めることができる。
山村は「術の式（の係数）に誤りがある」として，訂正の上で解を求めたが，「許容範囲内の誤差がある」と言っているが，そもそも山村が訂正した係数にも誤りがあるし，さらに言えば，係数は整数ではない。
どのような 6 次式なのか，手作業で方程式を解いてみる。
方針としては，eq1, eq2, eq4 の連立方程式から(a, r を含む） x0, x1, x2 （の式）を求めて，eq3 に代入して a, r を含む 式を求め，r について解く。 

res = solve([eq1, eq2, eq4], (x0, x1, x2))[3]  # 3 of 4

   ((a^3 - 4*a^2*(a/2 - sqrt(2)*r/2) + 2*a*r^2)/(-2*a^2 + 4*r^2), r*(1 + sqrt(2)), a/2 - sqrt(2)*r/2)

eq = eq3(x0 => res[1], x1 => res[2]) |> expand |> simplify
eq = numerator(eq)/a^2

    a^6/100 - 3*sqrt(2)*a^5*r/50 - a^5*r/25 + 3*sqrt(2)*a^4*r^2/25 + 6*a^4*r^2/25 - 2*sqrt(2)*a^3*r^3/25 - 8*sqrt(2)*a^2*r^4/25 - 11*a^2*r^4/25 + 2*sqrt(2)*a*r^5/25 + 4*a*r^5/25 + 4*r^6/25 + 4*sqrt(2)*r^6/25

a に定数を代入して eq = 0 を解けば r が求まる。
a = 10 としたとき，
38.62741699796952r^6+ 273.1370849898476r^5 -8925.483399593904r^4 -11313.70849898476r^3 +409705.6274847714r^2 -1.248528137423857e6r+ 1000000 = 0
である。係数は整数ではないので，「術」や山村のように整数係数を使用しても，正しい解は得られない。グラフを描けば，近似解すら得られないことがわかる。
SymPy でも解を得られないので，ニュートン・ラフソン法で数値解を求める。

function newton(fun, x1; delta = 1e-5, epsilon = 1e-14, maxrotation = 100)
   fun2(x) = (fun(x + delta) - fun(x)) / delta
   for i = 1:maxrotation
       x2 = x1 - fun(x1) / fun2(x1)
       if abs((x2 - x1) / x2) < epsilon
           return x2
       end
       x1 = x2
   end
   error("収束しませんでした")
end;

関数定義
fun(r) = 16*r^6 + 16*sqrt(2)*r^6 + 80*sqrt(2)*r^5 + 160*r^5 - 3200*sqrt(2)*r^4 - 4400*r^4 - 8000*sqrt(2)*r^3 + 120000*sqrt(2)*r^2 + 240000*r^2 - 600000*sqrt(2)*r - 400000*r + 1000000;

r = newton(fun , 1.368)
a = 10
x0 = a*(a^2 - 2*sqrt(2)*a*r - 2*r^2)/(2*(a^2 - 2*r^2))
x1 = r*(1 + sqrt(2))
x2 = a/2 - sqrt(2)*r/2
(r, x0, x1, x2)

   (1.368306828856772, 2.989643583187269, 3.30338490371374, 4.032460962571516)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 10
   r = a*(-4553*sqrt(-10 + 12√2) - 4605*sqrt(-5 + 6√2) - 3537 + 17080√2)/(2*(-4663*sqrt(-5 + 6√2) - 351√2*sqrt(-5 + 6√2) + 3845 + 4092√2))
   x0 = a*(-104√2 - 3 + 25*sqrt(-10 + 12√2) + 45*sqrt(-5 + 6√2))/(12*(-6√2 + 1 + 4*sqrt(-5 + 6√2)))
   x1 = (-4√2a*sqrt(-5 + 6√2) - 5*a*sqrt(-5 + 6√2) + 5√2a + 13*a)/(8 - 4*sqrt(-5 + 6√2))
   x2 = a*(-sqrt(2)/8 + 1/4 + sqrt(-5/32 + 3√2/16))
   @printf("正方形の一辺の長さが %g のとき，小円の直径は %g である。\n", a, 2r)
   @printf("a = %g;  r = %g;  x0 = %g;  x1 = %g;  x2 = %g\n", a, r, x0, x1, x2)
   #plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:blue, lw=0.5)
   plot([0, x0, a, a - x0, 0, a, a, 0, 0, a], [a, x0, 0, a - x0, a, 0, a, a, 0, 0], color=:blue, lw=0.5)
   segment(0, 0, x0, x0, :blue)
   segment(a, a, a - x0, a - x0, :blue)
   circle(x1, r, r)
   circle(r, x1, r)
   circle(x2, x2, r)
   circle(a - x1, a - r, r)
   circle(a - r, a - x1, r)
   circle(a - x2, a - x2, r)
   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(x0, x0, "(x0,x0)", :black, :left, :bottom, delta=delta/2)
       point(x1, r, "小円:r,(x1,r)", :red, :center, delta=-delta/2)
       point(x2, x2, "小円:r,(x2,x2)", :red, :center, delta=-delta/2)
   end
end;