算額（その197）

算額（その197）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（91）
長野県下高井郡木島平村穂高 一川谷大元神社 文化8年(1811)

問2 大円の中に 3 個の小円が互いに外接して入っている。小円の径が 5 寸のとき，大円の径はいかほどか。

大円半径を r とおき方程式を解く。

using SymPy

@syms r::positive;

eq1 = (√3(r - 5)/2)^2 + (3(r - 5)/2)^2 - (2*5)^2 
res = solve(eq1)[1]
res.evalf() |> println

   10.7735026918963

小円の径が 5 寸のとき，大円の径は 10寸7分7厘3毛5糸である。
算額では 17.77352 寸としているが，末尾桁は 0 または 17.773502 寸の写し間違いであろう。17.773502 寸としても，正しい末尾は 3 である。

なお，eq1 を簡約化すると 3r^2 - 30r - 25 = 0 である。
術では「4 を 3 で割り，平方根を求めて 1 を加え，小円の径を掛ける」とあるが，これは二次方程式の正の解を計算することを表している。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:blue, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = res[1]
   println("r = $r")
   plot()
   circle(0, 0, r)
   circle(0, r - 5, 5, :blue)
   circle(√3(r - 5)/2, (5 - r)/2, 5, :blue)
   circle(-√3(r - 5)/2, (5 - r)/2, 5, :blue)
   if more == true
       point(0, r - 5, " r-5")
       point(0, r, " r", :red)
       point(√3(r - 5)/2, (5 - r)/2, "√3(r-5)/2,(5-r)/2", :blue, :center, :top)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その196）

算額（その196）

中村信弥「改訂増補 長野県の算額」（91,165）
http://www.wasan.jp/zoho/zoho.html
長野県下高井郡木島平村穂高 一川谷大元神社 文化8年(1811)
長野県飯山市下木島 鳥出神社 天保14年(1843)

問1 2つの大円が交わっている部分に中円，それ以外の部分に小円が入っている。大円，中円の径が与えられたとき，小円の径を求めよ。

大円，中円，小円の半径を r1, r2, r3 とおく。
小円の中心の y 座標を y3 とおく。
以下の連立方程式を r3, y3 について解く。

using SymPy

@syms r1::positive, r2::positive, r3::positive, y3::positive;

eq1 = r3^2 + (r1 - r2 + y3)^2 - (r1 + r3)^2 |> expand
eq2 = r3^2 + (r2 - r1 + y3)^2 - (r1 - r3)^2 |> expand
solve([eq1, eq2], (r3, y3))

   1-element Vector{Tuple{Sym, Sym}}:
    (sqrt(r2)*(r1 - r2)*sqrt(2*r1 - r2)/r1, sqrt(r2)*sqrt(2*r1 - r2))

式を変形して，小円の径は「（1 - 中円の径/大円の径）× 平方根(2×大円の径×中円の径 - 中円の径の二乗）」である。当然であるが，文中の「径」は「直径」と読み替えても「半径」と読み替えてもよい。

中円の径が大円の径に比して小さいと，「小円が中円より大きい」ということになる（その限界値を求めるのも一興である）。

r1 = 10;  r2 = 4.56310987307924;  r3 = 4.5631098730792345;  y3 = 8.392867552141613

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(r1, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r3, y3) = (sqrt(r2)*(r1 - r2)*sqrt(2*r1 - r2)/r1, sqrt(r2)*sqrt(2*r1 - r2))
   println("r1 = $r1;  r2 = $r2;  r3 = $(r3);  y3 = $y3")
   plot()
   circle(0, 0, r2)
   circle(0, r1 - r2, r1, :green)
   circle(0, r2 - r1, r1, :green)
   circle(r3, y3, r3, :blue)
   circle(-r3, y3, r3, :blue)
   circle(r3, -y3, r3, :blue)
   circle(-r3, -y3, r3, :blue)
   if more == true
       point(0, r1 - r2, "r1-r2")
       point(0, r2 - r1, "r2-r1")
       point(r3, y3, "(r3,y3)", :blue)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その195）

算額（その195）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（91）
長野県下高井郡木島平村穂高 一川谷大元神社 文化8年(1811)

長方形内に2本の斜線を引き，4区分されたそれぞれの領域に大円，甲円，乙円，丙円を入れる。甲円，丙円の径がそれぞれ 126寸，66寸のとき，乙円の径はいかほどか。

長方形の長辺，短辺の長さをそれぞれ x, y とおく。その他，
大円の中心座標，半径を (x - r0, r0, r0)
甲円の中心座標，半径を (r1, r1, r1) ただし r1 = 63
乙円の中心座標，半径を (x2, y - r2, r2)
丙円の中心座標，半径を (x3, r3, r3) ただし r3 = 33
として，各円の中心から2本の斜線までの距離とそれぞれの円の半径の関係についての連立方程式を nlsolve() で解く。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
   p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
   l = sympy.Line(p1, p2)
   l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

@syms a::positive, b::positive, c::positive, x::positive, y::positive,
     x2::positive, r0::positive, r2::positive, x3::positive;

r1 = 126 // 2
r3 = 66 // 2
r0 = y / 2
eq1 = distance(a, 0, c, y, x- r0, r0) - r0^2
eq2 = distance(0, y, b, 0, x- r0, r0) - r0^2
eq3 = distance(a, 0, c, y, x2, y - r2) - r2^2
eq4 = distance(0, y, b, 0, x2, y - r2) - r2^2
eq5 = distance(a, 0, c, y, r1, r1) - r1^2
eq6 = distance(0, y, b, 0, r1, r1) - r1^2
eq7 = distance(a, 0, c, y, x3, r3) - r3^2
eq8 = distance(0, y, b, 0, x3, r3) - r3^2;

println(eq1, ",")
println(eq2, ",")
println(eq3, ",")
println(eq4, ",")
println(eq5, ",")
println(eq6, ",")
println(eq7, ",")
println(eq8, ",")

   -y^2/4 + (y/2 - y*(a^2 - 2*a*c + c^2 + y^2 + (a - c)*(a + c - 2*x + y))/(2*(a^2 - 2*a*c + c^2 + y^2)))^2 + (x - y/2 - (y^2*(a + c - 2*x + y)/2 + (2*x - y)*(a^2 - 2*a*c + c^2 + y^2)/2)/(a^2 - 2*a*c + c^2 + y^2))^2,
   -y^2/4 + (y/2 - y*(2*b^2 - 2*b*x + b*y + y^2)/(2*(b^2 + y^2)))^2 + (-b*(2*b*x - b*y + y^2)/(2*(b^2 + y^2)) + x - y/2)^2,
   -r2^2 + (x2 - (x2*(a^2 - 2*a*c + c^2 + y^2) + y*(a*r2 - c*r2 + c*y - x2*y))/(a^2 - 2*a*c + c^2 + y^2))^2 + (-r2 + y - y*(a^2 - a*c - a*x2 + c*x2 - r2*y + y^2)/(a^2 - 2*a*c + c^2 + y^2))^2,
   -r2^2 + (-b*(b*x2 + r2*y)/(b^2 + y^2) + x2)^2 + (-r2 + y - y*(b^2 - b*x2 - r2*y + y^2)/(b^2 + y^2))^2,
   (63 - (63*a^2 - 126*a*c + a*y^2 - 63*a*y + 63*c^2 + 63*c*y)/(a^2 - 2*a*c + c^2 + y^2))^2 + (-y*(a^2 - a*c - 63*a + 63*c + 63*y)/(a^2 - 2*a*c + c^2 + y^2) + 63)^2 - 3969,
   (-b*(63*b + y^2 - 63*y)/(b^2 + y^2) + 63)^2 + (-y*(b^2 - 63*b + 63*y)/(b^2 + y^2) + 63)^2 - 3969,
   (x3 - (x3*(a^2 - 2*a*c + c^2 + y^2) + y*(a*y - 33*a + 33*c - x3*y))/(a^2 - 2*a*c + c^2 + y^2))^2 + (-y*(a^2 - a*c - a*x3 + c*x3 + 33*y)/(a^2 - 2*a*c + c^2 + y^2) + 33)^2 - 1089,
   (-b*(b*x3 + y^2 - 33*y)/(b^2 + y^2) + x3)^2 + (-y*(b^2 - b*x3 + 33*y)/(b^2 + y^2) + 33)^2 - 1089,

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)
   (a, b, c, x, y, x2, r2, x3) = u
   return [
       -y^2/4 + (y/2 - y*(2*a^2 - 2*a*c - 2*a*x + a*y + 2*c*x - c*y + y^2)/(2*(a^2 - 2*a*c + c^2 + y^2)))^2 + (x - y/2 - (a^2*x - a^2*y/2 - 2*a*c*x + a*c*y + a*y^2/2 + c^2*x - c^2*y/2 + c*y^2/2)/(a^2 - 2*a*c + c^2 + y^2))^2,
       -y^2/4 + (y/2 - y*(2*b^2 - 2*b*x + b*y + y^2)/(2*(b^2 + y^2)))^2 + (-b*(2*b*x - b*y + y^2)/(2*(b^2 + y^2)) + x - y/2)^2,
       -r2^2 + (x2 - (x2*(a^2 - 2*a*c + c^2 + y^2) + y*(a*r2 - c*r2 + c*y - x2*y))/(a^2 - 2*a*c + c^2 + y^2))^2 + (-r2 + y - y*(a^2 - a*c - a*x2 + c*x2 - r2*y + y^2)/(a^2 - 2*a*c + c^2 + y^2))^2,
       -r2^2 + (-b*(b*x2 + r2*y)/(b^2 + y^2) + x2)^2 + (-r2 + y - y*(b^2 - b*x2 - r2*y + y^2)/(b^2 + y^2))^2,
       (63 - (63*a^2 - 126*a*c + a*y^2 - 63*a*y + 63*c^2 + 63*c*y)/(a^2 - 2*a*c + c^2 + y^2))^2 + (-y*(a^2 - a*c - 63*a + 63*c + 63*y)/(a^2 - 2*a*c + c^2 + y^2) + 63)^2 - 3969,
       (-b*(63*b + y^2 - 63*y)/(b^2 + y^2) + 63)^2 + (-y*(b^2 - 63*b + 63*y)/(b^2 + y^2) + 63)^2 - 3969,
       (x3 - (x3*(a^2 - 2*a*c + c^2 + y^2) + y*(a*y - 33*a + 33*c - x3*y))/(a^2 - 2*a*c + c^2 + y^2))^2 + (-y*(a^2 - a*c - a*x3 + c*x3 + 33*y)/(a^2 - 2*a*c + c^2 + y^2) + 33)^2 - 1089,
       (-b*(b*x3 + y^2 - 33*y)/(b^2 + y^2) + x3)^2 + (-y*(b^2 - b*x3 + 33*y)/(b^2 + y^2) + 33)^2 - 1089,     
   ]
end;

iniv = [big"85.0", 315, 310, 420, 170, 175, 45, 185]
res = nls(H, ini=iniv)
println(res)

   (BigFloat[84.0000000000000000000001735702499203810234362830123112659631373822954941885675, 314.9999999999999999999947583595807515206167302592172476055552893595369151008652, 307.9999999999999999999947727498833028947963179438839311145144054791595619767401, 419.9999999999999999999951927515441018951317502337587509125424356893062721470855, 168.0000000000000000000006020767375679450188933432128189504750790778747497805281, 175.999999999999999999997350763673708796543835009537017244584622091582997453893, 44.00000000000000000000017089021421936614679777263092841262631968462778651775099, 182.9999999999999999999975156955872708477881618570615211838481194365997674781264], true)

   a = 84.00000, b = 315.00000, c = 308.00000
   x = 420.00000, y = 168.00000, x2 = 176.00000, r2 = 44.00000, x3 = 183.00000
   乙円の直径 = 88.00000

乙円の半径は 44 なので，直径は 88 寸である。
図に描いてみると，算額の図は正確なものではない（与えられた条件を満たすものではない）ことがわかる。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, c, x, y, x2, r2, x3) = res[1]
   r0 = y/2
   @printf("a = %.5f, b = %.5f, c = %.5f\nx = %.5f, y = %.5f, x2 = %.5f, r2 = %.5f, x3 = %.5f\n",
       a, b, c, x, y, x2, r2, x3)
   @printf("乙円の直径 = %.5f\n", 2r2)
   plot([0, x, x, 0, 0], [0, 0, y, y, 0], color=:black, lw=0.5)
   circle(x - r0, r0, r0)
   circle(r1, r1, r1, :blue)
   circle(x2, y - r2, r2, :green)
   circle(x3, r3, r3, :orange)
   segment(a, 0, c, y)
   segment(0, y, b, 0)
   if more == true
       point(x - r0, r0, "大:(x-r0,r0,r0)", :red, :center)
       point(r1, r1, "甲:(r1,r1,r1)", :blue, :center)
       point(x2, y - r2, "乙:(x2,y-r2,r2)", :green, :center)
       point(x3, r3, "丙:(x3,r3,r3)", :orange, :center)
       point(0, y, " y", :black, :left, :bottom)
       point(x, 0, " x", :black, :left, :bottom)
       point(a, 0, "a ", :black, :right, :bottom)
       point(b, 0, " b", :black, :left, :bottom)
       point(c, y, " c", :black, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;