算額（その65）

算額（その65）

長野県木曽郡上松町寝覚 臨川寺弁財天社 文政13年
http://www.wasan.jp/nagano/rinsenji.html

3 個の甲円が交差している部分に内接する乙円の径を求めよ。

甲円の半径を 1 として，図のように記号を定め，方程式を解く。
x, y は簡単に決められる。方程式は，甲円と乙円が内接することに基づく。

using SymPy
@syms x::positive, y::positive, r::positive;
x = sqrt(Sym(3))/4
y = 1//4
eq = x^2 + (1 - y)^2 - (1 - r)^2;

res = solve(eq) |> println

   Sym[1 - sqrt(3)/2, sqrt(3)/2 + 1]

解は 2 つ得られるが，1 - sqrt(3)/2 = 0.1339745962155614 が適切な解である。

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()
   r = 1 - sqrt(3)/2
   println("r = $r")
   circle(0, 1, 1, :green)
   circle(cos(pi/6), -sin(pi/6), 1, :green)
   circle(-cos(pi/6), -sin(pi/6), 1, :green)
   circle(x, y, r)
   circle(-x, y, r)
   circle(0, -1/2, r)
   if more
       point(x, y, "(x,y)", :red, :center)
       point(0, 1, "1 ", :green, :right)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

   r = 0.1339745962155614

算額（その64）

算額（その64）

長野県飯山市滝澤家所蔵
http://www.wasan.jp/nagano/takizawa.html

大円の中に正方形，中円，小円が収まっている。大円の径を 5 寸としたとき，中円の中心を結ぶ正方形の一辺の長さと小円の径を求めよ。

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

using SymPy
@syms r1::positive, r2::positive;
eq1 = 2r2 - 5/sqrt(Sym(2));
eq2 = 2r2^2 - (r1 + r2)^2;
res = solve([eq1, eq2], (r1, r2));
res[1][1] |> println  # r1
res[1][2] |> println  # r2

   5/2 - 5*sqrt(2)/4
   5*sqrt(2)/4

正方形の一辺は 2r2。算額では 3寸5分3厘としている。

2 * 5*sqrt(2)/4

   3.5355339059327378

小円の径は 5/2 - 5*sqrt(2)/4。算額では 7分8厘としている。

5/2 - 5*sqrt(2)/4 |> N

   0.7322330470336311

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) = (5/2 - 5*sqrt(2)/4, 5*sqrt(2)/4)
   println("r1 = $r1, r2 = $r2")
   circle(0, 0, 5, :black)
   circle(0, 0, r1)
   circle(r2, r2, r2, :brown)
   circle(r2, -r2, r2, :brown)
   circle(-r2, r2, r2, :brown)
   circle(-r2, -r2, r2, :brown)
   x = 5/sqrt(2)
   plot!([-x, x, x, -x, -x], [-x, -x, x, x, -x], color=:black, lw=0.5)
   plot!([-r2, r2, r2, -r2, -r2], [-r2, -r2, r2, r2, -r2], color=:blue, lw=0.5, linestyle=:dot)
   if more
       point(0, 0, "0 ", :black, :right)
       point(r1, 0, "r1 ", :black, :right)
       point(r2, 0, "r2 ", :brown, :right)
       point(5/√2, 0, "5/√2 ", :brown, :right)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

   r1 = 0.7322330470336311, r2 = 1.7677669529663689

 

算額（その63）

算額（その63）

飯山市常磐柳新田 瀧澤治夫宅所蔵 天保11年(1840)7月
http://www.wasan.jp/nagano/takizawa2.html

半円の中に，3 種類の円が収まっている。各円の径を求めよ。

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

using SymPy
@syms x1::positive, r1::positive, x2::positive, r2::positive;
eq1 = x1^2 + r1^2 - (1 - r1)^2;
eq2 = x2^2 + r2^2 - (1 - r2)^2;
eq3 = x1^2 + (1//2 - r1)^2 - (1//2 + r1)^2;
eq4 = (x2 - x1)^2 + (r1 - r2)^2 - (r1 + r2)^2;

res = solve([eq1, eq2, eq3, eq4], (x1, r1, x2, r2))

   1-element Vector{NTuple{4, Sym}}:
    (sqrt(2)/2, 1/4, 2*sqrt(2)/3, 1/18)

算額では，甲円は1尺2寸（12寸），乙円は8寸，丙円は結果として 3 寸としている。しかし，算額の図はこの比率になっていない（図のほうが正しい）。

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()
   (x1, r1, x2, r2) = (sqrt(2)/2, 1/4, 2*sqrt(2)/3, 1/18)
   println("x1 = $x1, r1 = $r1\nx2 = $x2, r2 = $r2")
   circle(0, 0, 1, :black)
   circle(0, 1/2, 1/2)
   circle(x1, r1, r1, :brown)
   circle(-x1, r1, r1, :brown)
   circle(x2, r2, r2, :green)
   circle(-x2, r2, r2, :green)
   hline!([0], color=:black, lw=0.5)
   if more
       point(0, 0, "0 ", :black, :right)
       point(0, 1/2, "1/2 ", :black, :right)
       point(x1, r1, "(x1,r1)", :brown, :center)
       point(x2, r2, "\n    (x2,r2)", :green, :center)
       vline!([0], color=:black, lw=0.5)
   end
end;

   x1 = 0.7071067811865476, r1 = 0.25
   x2 = 0.9428090415820635, r2 = 0.05555555555555555

 

算額（その62）

算額（その62）

長野市篠ノ井 布制神社 文化3年(1806)
http://www.wasan.jp/nagano/fuse.html

玄が12寸，甲円の径が2寸5分のとき丙円の径を求めよ。

玄を 240，甲円の半径を 25 とし，図のように記号を定め，方程式を解く。

using SymPy
@syms x1::positive, r1::positive, x2::positive, r2::positive, y0::negative, r0::positive;
eq1 = x1^2 + (25 - r1)^2 - (25 + r1)^2;
eq2 = x2^2 + (25 - r2)^2 - (25 + r2)^2;
eq3 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2;
eq4 = y0^2 + 120^2 - r0^2;
eq5 = x1^2 + (r1 - y0)^2 - (r0 - r1)^2;
eq6 = 50 - y0 - r0;

res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (x1, r1, x2, r2, y0, r0))

   2-element Vector{NTuple{6, Sym}}:
    (600/13, 3600/169, 24, 144/25, -119, 169)
    (600/13, 3600/169, 600, 3600, -119, 169)

丙円の半径は算額では 5分7厘5毛となっているが，正しくは 5分7厘6毛（5.76）となる。

また，算額の各円の見かけの大きさは，計算結果とはかなり異なっている。

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()
   (x1, r1, x2, r2, y0, r0) = (600/13, 3600/169, 24, 144/25, -119, 169)
   println("x1 = $x1, r1 = $r1\nx2 = $x2, r2 = $r2\ny0 = $y0, r0 = $r0")
   circle(0, y0, r0, :black)
   circle(0, 25, 25)
   circle(x1, r1, r1, :brown)
   circle(-x1, r1, r1, :brown)
   circle(x2, r2, r2, :green)
   plot!([-120, 120], [0, 0], color=:black, lw=0.5)
   if more
       point(0, 0, "0 ", :black, :right)
       point(0, 25, "25 ", :black, :right)
       point(120, 0, " 120", :black)
       point(x1, r1, "(x1,r1)", :brown, :center)
       point(x2, r2, "\n    (x2,r2)", :green, :center)
       point(0, y0, "y0 ", :blue, :right)
       vline!([0], color=:black, lw=0.5)
   end
end;

   x1 = 46.15384615384615, r1 = 21.301775147928993
   x2 = 24, r2 = 5.76
   y0 = -119, r0 = 169

算額（その61）

算額（その61）

長野市鬼無里 智光山文珠堂 明治11年(1878)8月
http://www.wasan.jp/nagano/monjudo1.html

半円の中に 5 個の円がある。それぞれの径を求めよ。

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

using SymPy
@syms x1::positive, r1::positive, x2::negative, r2::positive,
   x3::positive, r3::positive, x4::negative, r4::positive,
   x5::positive, r5::positive;
eq1 = x1^2 + r1^2 - (1 - r1)^2
eq2 = x3^2 + r3^2 - (1 - r3)^2
eq3 = x2^2 + r2^2 - (1 - r2)^2
eq4 = (x1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq5 = (x1 - x5)^2 + (r1 - r5)^2 - (r1 + r5)^2
eq6 = (x5 - x3)^2 + (r5 - r3)^2 - (r5 + r3)^2
eq7 = (x2 - x1)^2 + (r2 - r1)^2 - (r2 + r1)^2
eq8 = (x4 - x1)^2 + (r4 - r1)^2 - (r4 + r1)^2
eq9 = (x2 - x4)^2 + (r2 - r4)^2 - (r2 + r4)^2
eq10 = 1/r1 + 1/r3 + 2sqrt(1/r1/r3) - 1/r5;  # デカルトの円定理

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

この問題の場合も solve() では解けないので nlsolve() を用いる。

eq1 |> expand |> println
eq2 |> expand |> println
eq3 |> expand |> println
eq4 |> expand |> println
eq5 |> expand |> println
eq6 |> expand |> println
eq7 |> expand |> println
eq8 |> expand |> println
eq9 |> expand |> println
eq10 |> expand |> println

   2*r1 + x1^2 - 1
   2*r3 + x3^2 - 1
   2*r2 + x2^2 - 1
   -4*r1*r3 + x1^2 - 2*x1*x3 + x3^2
   -4*r1*r5 + x1^2 - 2*x1*x5 + x5^2
   -4*r3*r5 + x3^2 - 2*x3*x5 + x5^2
   -4*r1*r2 + x1^2 - 2*x1*x2 + x2^2
   -4*r1*r4 + x1^2 - 2*x1*x4 + x4^2
   -4*r2*r4 + x2^2 - 2*x2*x4 + x4^2
   -1/r5 + 1/r3 + 1/r1 + 2/(sqrt(r1)*sqrt(r3))

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, r1, x2, r2, x3, r3, x4, r4, x5, r5 = u
return [
       2*r1 + x1^2 - 1
       2*r3 + x3^2 - 1
       2*r2 + x2^2 - 1
       -4*r1*r3 + x1^2 - 2*x1*x3 + x3^2
       -4*r1*r5 + x1^2 - 2*x1*x5 + x5^2
       -4*r3*r5 + x3^2 - 2*x3*x5 + x5^2
       -4*r1*r2 + x1^2 - 2*x1*x2 + x2^2
       -4*r1*r4 + x1^2 - 2*x1*x4 + x4^2
       -4*r2*r4 + x2^2 - 2*x2*x4 + x4^2
       -1/r5 + 1/r3 + 1/r1 + 2/(sqrt(r1)*sqrt(r3))
       ];
end;

iniv = [0.01, 0.5, -0.7, 0.22, 0.8, 0.19, -0.5, 0.05, 0.4, 0.05];
res = nls(H, ini=iniv)

   ([0.21236918980148334, 0.47744966361153107, -0.5821590777643445, 0.3305454040882841, 0.7994277491932246, 0.18045763690992736, -0.22131238148404836, 0.0984814314508113, 0.5759211572513911, 0.06920626565999469], true)

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(ylims=(-0.07, 1.07))
   x1, r1, x2, r2, x3, r3, x4, r4, x5, r5 = [0.21236918980148334, 0.47744966361153107, -0.5821590777643445, 0.3305454040882841, 0.7994277491932246, 0.18045763690992736, -0.22131238148404836, 0.0984814314508113, 0.5759211572513911, 0.06920626565999469]
   println("x1 = $x1, r1 = $r1\nx2 = $x2, r2 = $r2\nx3 = $x3, r3 = $r3\nx4 = $x4, r4 = $r4\nx5 = $x5, r5 = $r5")
   circle(0, 0, 1, :black)
   circle(x1, r1, r1, :brown)
   circle(x2, r2, r2, :green)
   circle(x3, r3, r3, :blue)
   circle(x4, r4, r4, :red)
   circle(x5, r5, r5, :magenta)
   hline!([0], color=:black, lw=0.5)
   if more
       point(0, 0, " 0", :black)
       point(x1, r1, "(x1,r1)", :brown, :center)
       point(x2, r2, "(x2,r2)", :green, :center)
       point(x3, r3, "(r3,r3)", :blue, :center)
       point(x4, r4, "(x4,r4)", :red, :center)
       point(x5, r5, "(x5,r5)", :magenta, :center)
       vline!([0], color=:black, lw=0.5)
   end
end;

   x1 = 0.21236918980148334, r1 = 0.47744966361153107
   x2 = -0.5821590777643445, r2 = 0.3305454040882841
   x3 = 0.7994277491932246, r3 = 0.18045763690992736
   x4 = -0.22131238148404836, r4 = 0.0984814314508113
   x5 = 0.5759211572513911, r5 = 0.06920626565999469