算額（その60）

算額（その60）

石川県松任市若宮 若宮八幡宮 文政13年吉祥日
http://www.wasan.jp/isikawa/wakamiyahatiman1.html

石川県鳳珠郡穴水町 美麻奈比古神社 文政13年8月
http://www.wasan.jp/isikawa/mimanahiko.html

群馬県高崎市 祖師堂 明治16年
http://www.wasan.jp/gunma/sosido.html

外円の中に 8 個の小円と 1 個の大円がある。それぞれの径を求めよ。

外円の半径を 1 とし，図のように記号を定め，方程式を解く。

using SymPy
@syms r1::positive, x2::positive, y2::positive, r2::positive;

x2 = (1 - r2) / sqrt(Sym(2))
y2 = x2
eq1 = r1 + 2r2 - 1 |> expand
eq2 = x2^2 + y2^2 - (1 - r2)^2 |> expand
eq3 = x2^2 + (1 - r2 - y2)^2 - 4r2^2 |> expand;

res = solve([eq1, eq3], (r1, r2))

   2-element Vector{Tuple{Sym, Sym}}:
    (-4*sqrt(2) - 2*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 7, -3 + sqrt(2)*sqrt(10 - 7*sqrt(2)) + 2*sqrt(2))
    (-4*sqrt(2) + 2*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 7, -3 - sqrt(2)*sqrt(10 - 7*sqrt(2)) + 2*sqrt(2))

解が 2 組得られるが，最初のものが適切な解である。

using Plots

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=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   (r1, r2) = (-4*sqrt(2) - 2*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 7, -3 + sqrt(2)*sqrt(10 - 7*sqrt(2)) + 2*sqrt(2))
   x2 = y2 = (1 - r2)/sqrt(2)
   println("r1 = $r1, x2 = $x2, y2 = $y2, r2 = $r2")
   circle(0, 0, 1, :black)
   circle(0, 0, r1, :brown)
   circle(x2, y2, r2)
   circle(x2, -y2, r2)
   circle(-x2, y2, r2)
   circle(-x2, -y2, r2)
   circle(0, 1-r2, r2, :green)
   circle(0, r2-1, r2, :green)
   circle(1-r2, 0, r2, :magenta)
   circle(r2-1, 0, r2, :magenta)
   if more
       point(0, 0, "0 ", :black, :right)
       point(0, 1-r2, "1-r2 ", :green, :right)
       point(0, r1, "r1 ", :brown, :right)
       point(1-r2, 0, "1-r2 ", :green, :right)
       point(r1, 0, "r1 ", :brown, :right)
       point(x2, y2, "(x2,y2)", :red, :center)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   end
end;

   r1 = 0.4464626921716892, x2 = 0.5114017891839755, y2 = 0.5114017891839755, r2 = 0.2767686539141554

算額（その59）

算額（その59）

静岡県掛川市大坂 観音堂 明治28年
http://www.wasan.jp/sizuoka/kanondo.html

千葉胤秀: 算法新書，文政13年(1830)
一関市博物館>>和算に挑戦>>平成22年度出題問題&解答例>>平成22年度出題問題(3)[上級問題]
https://www.city.ichinoseki.iwate.jp/museum/wasan/h22/normal.html

直角三角形の中に 5 つの円が入っている。甲円，乙円の径をそれぞれ 2，1 としたとき，中円の径を求めよ。

以下のように記号を定め方程式を解く。

using SymPy
@syms x1::positive, x2::positive, x3::positive, x4::positive;
@syms y2::positive, y4::positive, y5::positive;
@syms w::positive, h::positive;

eq1 = w^2 + h^2 - (w + h - 2x1)^2 |> expand
eq2 = x1/(w - x1) - y2/(w - x2) |> expand
eq3 = x1/(w - x1) - 2/(w - x3) |> expand
eq4 = x1/(h - x1) - x4/(h - y4) |> expand
eq5 = x1/(h - x1) - 1/(h - y5) |> expand
eq6 = (x1 - x2)^2 + (x1 - y2)^2 - (x1 + y2)^2 |> expand
eq7 = (x2 - x3)^2 + (y2 - 2)^2 - (y2 + 2)^2 |> expand
eq8 = (x1 - x4)^2 + (y4 - x1)^2 - (x1 + x4)^2 |> expand
eq9 = (x4 - 1)^2 + (y5 - y4)^2 - (1 + x4)^2 |> expand;

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9])

この問題でも 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)
 x1, x2, x3, x4, y2, y4, y5, w, h = u
 return [
       -2*h*w + 4*h*x1 + 4*w*x1 - 4*x1^2,
       x1/(w - x1) - y2/(w - x2),
       x1/(w - x1) - 2/(w - x3),
       x1/(h - x1) - x4/(h - y4),
       x1/(h - x1) - 1/(h - y5),
       x1^2 - 2*x1*x2 - 4*x1*y2 + x2^2,
       x2^2 - 2*x2*x3 + x3^2 - 8*y2,
       x1^2 - 4*x1*x4 - 2*x1*y4 + y4^2,
       -4*x4 + y4^2 - 2*y4*y5 + y5^2
       ];
end;

iniv = [3.33, 8.33, 10.83, 1.25, 1.67, 7.08, 8.33, 14.58, 9.58];
res = nls(H, ini=iniv)

   ([7.13576989434318, 17.51985221378425, 23.017320107299383, 2.6712861872781755, 3.777769155028713, 15.867707492460653, 19.136521459194263, 29.202001485041404, 21.09238873771053], false)

収束フラッグは false になっているが，十分な精度で解が求まっている。

中円の径は 7.13576989434318 となった。算額では，7.13576989446576 である。

using Plots

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=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot(xlims=(-2,32), ylims=(-2,24))
   (x1, x2, x3, x4, y2, y4, y5, w, h) = [7.13576989434318, 17.51985221378425, 23.017320107299383, 2.6712861872781755, 3.777769155028713, 15.867707492460653, 19.136521459194263, 29.202001485041404, 21.09238873771053]
   y1 = x1
   println("x1 = $x1, x2 = $x2, x3 = $x3, x4 = $x4")
   println("y2 = $y2, y4 = $y4, y5 = $y5")
   println("w = $w, h = $h")
   plot!([0, w, 0, 0], [h, 0, 0, h], color=:black, lw=0.5)
   circle(x1, y1, x1, :black)
   circle(x2, y2, y2)
   circle(x3, 2, 2, :green)
   circle(x4, y4, x4, :magenta)
   circle(1, y5, 1, :brown)
   if more
       point(x1, x1, "(x1,x1)", :black, :center)
       point(x2, y2, "(x2,y2)", :red, :center)
       point(x3, 2, "(x3,2)", :green, :center)
       point(x4, y4, "(x4,y4)", :magenta, :center)
       point(1, y5, "(1,y5)", :brown)
       point(w, 0, " w", :black)
       point(0, h, "h ", :black, :right)
   end
end;

   x1 = 7.13576989434318, x2 = 17.51985221378425, x3 = 23.017320107299383, x4 = 2.6712861872781755
   y2 = 3.777769155028713, y4 = 15.867707492460653, y5 = 19.136521459194263
   w = 29.202001485041404, h = 21.09238873771053