算額（その131）

算額（その131）

新潟県長岡市与板町本与板 与板八幡宮 文化5年(1808)
http://www.wasan.jp/niigata/yoitahatiman3.html

短辺（平）の長さが 595 寸の長方形の中に，甲円，乙円，丙円，丁円，戊円が入っている。甲円の径はいかほどか。

長辺（長）の長さを w, 甲円，乙円，丙円，丁円，戊円の半径を r1, r2, r3, r4, r5，戊円の中心の x 座標を x5 として以下のように記号を定め，方程式を立てて，解く。

using SymPy

@syms r1::positive, r2::positive, r3::positive, r4::positive, r5::positive, x5::positive, w::positive;

eq1 = (w - r1 - r2)^2 + (595 - r1 - r2)^2 - (r1 + r2)^2
eq2 = (w - r1 - r3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (r4 - r1)^2 + (595 - r1 - r4)^2 - (r1 + r4)^2
eq4 = (w - r1 - x5)^2 + (595 - r1 - r5)^2 - (r1 + r5)^2
eq5 = (r2 - r3)^2 + (r2 - 595 + r3)^2 - (r2 + r3)^2
eq6 = (r2 - x5)^2 + (r2 - r5)^2 - (r2 + r5)^2
eq7 = (x5 - w + r4)^2 + (r4 - r5)^2 - (r4 + r5)^2;

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r1, r2, r3, r4, r5, x5, w))

これも solve() では解けないので，nlsolve() で解く。

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

   -(r1 + r2)^2 + (-r1 - r2 + 595)^2 + (-r1 - r2 + w)^2,
   (r1 - r3)^2 - (r1 + r3)^2 + (-r1 - r3 + w)^2,
   (-r1 + r4)^2 - (r1 + r4)^2 + (-r1 - r4 + 595)^2,
   -(r1 + r5)^2 + (-r1 - r5 + 595)^2 + (-r1 + w - x5)^2,
   (r2 - r3)^2 - (r2 + r3)^2 + (r2 + r3 - 595)^2,
   (r2 - r5)^2 - (r2 + r5)^2 + (r2 - x5)^2,
   (r4 - r5)^2 - (r4 + r5)^2 + (r4 - w + x5)^2,

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)
(r1, r2, r3, r4, r5, x5, w) = u
return [
   -(r1 + r2)^2 + (-r1 - r2 + 595)^2 + (-r1 - r2 + w)^2,
   (r1 - r3)^2 - (r1 + r3)^2 + (-r1 - r3 + w)^2,
   (-r1 + r4)^2 - (r1 + r4)^2 + (-r1 - r4 + 595)^2,
   -(r1 + r5)^2 + (-r1 - r5 + 595)^2 + (-r1 + w - x5)^2,
   (r2 - r3)^2 - (r2 + r3)^2 + (r2 + r3 - 595)^2,
   (r2 - r5)^2 - (r2 + r5)^2 + (r2 - x5)^2,
   (r4 - r5)^2 - (r4 + r5)^2 + (r4 - w + x5)^2,
   ]
end;

iniv = [200.0, 165, 135, 95, 90, 400, 700]

res = nls(H, ini=iniv)
println(res)

   ([213.00034502265584, 164.1726947717987, 134.08788967981687, 96.00257736188509, 88.286549252726, 404.9567532384933, 685.0868546118494], false)

甲円の半径 = 213.000,  乙円の半径 = 164.173,  丙円の半径 = 134.088,  丁円の半径 = 96.003,  戊円の半径 = 88.287
長 = 685.087,  x = 404.957
となるので，算額の問である「甲円の径］は 213*2 = 426寸である。

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(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, r4, r5, x5, w) = res[1]
   @printf("甲円の半径 = %.3f,  乙円の半径 = %.3f,  丙円の半径 = %.3f,  丁円の半径 = %.3f,  戊円の半径 = %.3f\n", r1, r2, r3, r4, r5)
   @printf("長 = %.3f,  x = %.3f\n", w, x5)
   plot([0, w, w, 0, 0], [0, 0, 595, 595, 0])
   circle(w - r1, 595 - r1, r1)
   circle(r2, r2, r2)
   circle(r3, 595 - r3, r3)
   circle(w - r4, r4, r4)
   circle(x5, r5, r5)
   if more == true
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

   甲円の半径 = 213.000,  乙円の半径 = 164.173,  丙円の半径 = 134.088,  丁円の半径 = 96.003,  戊円の半径 = 88.287
   長 = 685.087,  x = 404.957

 

 

算額（その130）

算額（その130）

算額（その51）と上下反転しているが同じである。51では未知数を絞り，代数解を求めている。

群馬県桐生市梅田町 日枝神社 明治12年(1879)10月
http://www.wasan.jp/gunma/hieda.html

円内に4種類の円が入っている。甲円の径が与えられたとき，外円の径を求めよ。

外円の半径を 1 とし，図形が点対称なので -30度〜90度の範囲に限定する。また以下のように記号を定め，方程式を立てて，解く。
甲円の中心座標，半径　0, 1-r1, r1
乙円の中心座標，半径　x2, y2, r2
丙円の中心座標，半径　x3, y3, r3
丁円の中心座標，半径　0, 0, r4

using SymPy

@syms r1::positive, x2::positive, y2::positive, r2::positive, x3::positive, y3::positive, r3::positive, r4::positive;

eq1 = x2^2 + (1 - r1 - y2)^2 - (r1 + r2)^2  # 甲乙
eq2 = x3^2 + (1 - r1 - y3)^2 - (r1 + r3)^2  # 甲丙
eq3 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2  # 乙丙
eq4 = x3^2 + y3^2 - (r3 + r4)^2  # 丙丁
eq5 = x2^2 + y2^2 - ( 1 - r2)^2  # 外乙　（内接）
eq6 = r4 + 2r1 - 1  # 半径の関係
eq7 = r4 + 2r3 + 2r2 - 1  # 半径の関係
eq8 = y3/x3 - y2/x2;  # 乙円，丙円，丁円の中心が一直線上にある

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (r1, x2, y2, r2, x3, y3, r3, r4))

solve() では解けないので，nlsolve() を用いる。

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)
(r1, x2, y2, r2, x3, y3, r3, r4) = u
return [
   x2^2 - (r1 + r2)^2 + (-r1 - y2 + 1)^2,
   x3^2 - (r1 + r3)^2 + (-r1 - y3 + 1)^2,
   -(r2 + r3)^2 + (x2 - x3)^2 + (y2 - y3)^2,
   x3^2 + y3^2 - (r3 + r4)^2,
   x2^2 + y2^2 - (1 - r2)^2,
   2*r1 + r4 - 1,
   2*r2 + 2*r3 + r4 - 1,
   y3/x3 - y2/x2,    
   ]
end;

iniv = [0.40, 0.6, 0.36, 0.27, 0.28, 0.16, 0.12, 0.20]

res = nls(H, ini=iniv)
println(res)

   ([0.39728455784707767, 0.6275824608249739, 0.36622451524418415, 0.27337758037549836, 0.28444851157211404, 0.16598937154089602, 0.12390697747157697, 0.20543088430584824], true)

外円の径は，甲円の径の 1/0.39728455784707767 ≒ 2.5 倍である（正確には 1/(-7/9 + 4*sqrt(7)/9) = 2.511857892036909）。

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 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, color=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")
   (r1, x2, y2, r2, x3, y3, r3, r4) = res[1]
   @printf("甲円の半径 = %.3f,  乙円の半径 = %.3f,  丙円の半径 = %.3f, 丁円の半径 = %.3f\n", r1, r2, r3, r4)
   plot()
   circle(0, 0, 1, :black)
   circle(0, 1 - r1, r1)
   circle((1 - r1)√3/2, -(1 - r1)/2, r1)
   circle(-(1 - r1)√3/2, -(1 - r1)/2, r1)
   circle(x2, y2, r2, :blue)
   circle(-x2, y2, r2, :blue)
   circle(0, r2 - 1, r2, :blue)
   circle(x3, y3, r3, :green)
   circle(-x3, y3, r3, :green)
   circle(0, -r4 - r3, r3, :green)
   circle(0, 0, r4)
   if more == true
       point(0, 1-r1, "甲:(0,1-r1,r1)", :red)
       point(x2, y2, "乙:(x2,y2,r2)", :blue)
       point(x3, y3, "丙:(x3,y3,r3)", :green, :left, :bottom)
       point(0, 0, "丁:(0,0,r4)", :green, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

   甲円の半径 = 0.397,  乙円の半径 = 0.273,  丙円の半径 = 0.124, 丁円の半径 = 0.205