算額（その43）

算額（その43）

埼玉県加須市騎西町 玉敷神社

足柄上群大井町 三嶋神社 明治15年
http://www.wasan.jp/kanagawa/misima.html

2024/01/12 追記：日本橋 福徳神社（芽吹稲荷）に埼玉県かぞ算額文化保存会が奉納した 2 枚目の算額の問1として加須市騎西町 玉敷神社の算額問題が掲載されている。https://mebuki.jp/2125/

図のように，外円の中に甲，乙，丙の円が入っている。それぞれの円の径を求めよ。

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

using SymPy

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

eq2 = (r2 + r3)^2 + ((r2 + r3) + r3)^2 - ((r2 + r3) + r2)^2
eq3 = (r2 + r3)^2 + ((r2 + r3) + r3)^2 - (1 - r3)^2;

res = solve([eq2, eq3], (r2, r3))

   1-element Vector{Tuple{Sym, Sym}}:
    (1/3, 1/6)

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")
   (r2, r3) = (1/3, 1/6)
   r1 = r2 + r3
   x2 = r1 + r3
   println("r1 = $r1;  r2 = $r2;  r3 = $r3;  x2 = $x2")
   plot()
   circle(0, 0, 1, :black)
   circle(0, r1, r1, :red)
   circle(0, -r1, r1, :red)
   circle(x2, 0, r2, :blue)
   circle(-x2, 0, r2, :blue)
   circle(x2, r2+r3, r3, :brown)
   circle(x2, -(r2+r3), r3, :brown)
   circle(-x2, r2+r3, r3, :brown)
   circle(-x2, -(r2+r3), r3, :brown)
   if more
       point(0, r1, "甲円 r1", :red)
       point(x2, 0, "乙円 r2", :blue)
       point(r1+r3, r2+r3, "丙円 r3", :red)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
end;

算額（その42）

算額（その42）

山形県鶴岡市羽黒町 出羽三山神社 文政6年
http://www.wasan.jp/yamagata/haguro.html

外円内に大円 3 個と小円 4 個が収まっている。それぞれの径を求めよ。

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

using SymPy

@syms r1::real, r2::real, x1::positive, y1::negative, x2::positive, y2::negative;

x1 = (1 - r1)sqrt(Sym(3))/2
y1 = (r1 - 1)/2    

x2 = (r2 - r1)sqrt(Sym(3))/2
y2 = (r1 - r2)/2

eq1 = x1^2 + (1 - r2 - y1)^2 - (r1 + r2)^2
eq2 = x2^2 + (1 - r1 - y2)^2 - (r1 + r2)^2

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

   3-element Vector{Tuple{Sym, Sym}}:
    (1, 0)
    (-2*sqrt(7)/3 - 2*(-2 + sqrt(7))^2/3 + 7/3, -2 + sqrt(7))
    (-2*(-sqrt(7) - 2)^2/3 + 2*sqrt(7)/3 + 7/3, -sqrt(7) - 2)

2 番目の解が適。

(-2*sqrt(7)/3 - 2*(-2 + sqrt(7))^2/3 + 7/3, -2 + sqrt(7))

簡約化すると

(2*sqrt(7) - 5, sqrt(7) - 2)

   (0.291502622129181, 0.6457513110645907)

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")
   (r1, r2) = (2*sqrt(7) - 5, sqrt(7) - 2)
   x1 = (1 - r1)sqrt(3)/2
   y1 = (r1 - 1)/2   
   println("r1 = $r1;  r2 = $r2;  x1 = $x1;  y1 = $y1")
   plot()
   circle(0, 0, 1, :black)
   circle(0, 0, r1, :blue)
   circle(0, 1-r1, r1, :blue)
   circle((1-r1)sqrt(3)/2, -(1-r1)/2, r1, :blue)
   circle(-(1-r1)sqrt(3)/2, -(1-r1)/2, r1, :blue)
   x2 = (r2 - r1)sqrt(Sym(3))/2
   y2 = (r1 - r2)/2
   circle(x2, y2, r2, :magenta)
   circle(-x2, y2, r2, :magenta)
   circle(0, 1-r2, r2, :magenta)
   if more
       point(0, 1-r1, "1-r1 ", :blue)
       point(0, 1-r2, "1-r2 ", :magenta)
       point((1-r2)sqrt(3)/2, -(1-r2)/2, "  (x2,y2)", :magenta)
       point(x1, y1, "(x1,y1)", :blue, :center)
       point(0, 0, "0 ", :magenta)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
end;

 

算額（その41）

算額（その41）

大阪府茨木市 井於神社 弘化3年
http://www.wasan.jp/osaka/iyo.html

3辺が 3,4,5の直角三角形に，図のように，大円，中円，小円が入っている。それぞれの円の径を求めよ。

図のように，それぞれの円の中心座標と半径を (0, r1, r1), (x2, r2, r2), (x3, r3, r3), 大円と斜辺の接点の座標を (x1, y1) とする。

using SymPy

@syms r1::positive, r2::positive, r3::positive;
@syms x2::positive, x3::positive, x0::positive, x1::positive, y1::positive, l::positive;

まず，AB = AE = 3 - r1, BC = CD = 4 - r1, AB + BC = 5  から 7 - 2r1 = 5。よって r1 = 1 である。

eq1 = (3 - r1) + (4 - r1) - 5
solve(eq1, r1)[1] |> println

   1

次に，以下の方程式を解く。

r1 = 1;
eq2 = (x3 - x2)^2 + (r2 - r3)^2 - (r2 + r3)^2;
eq3 = (4 - x3)/(4 - x2) - r3/r2;
eq4 = (4 - r1)/(4 - x2) - r1/r2;
eq6 = (x3 - r1)^2 + (r1 - r3)^2 - (r1 + 2r2 + r3)^2;
eq8 = 4y1 - 3(4 - x1);
eq9 = (x1 - r1)^2 + (y1 - r1)^2 - r1^2;

res = solve([eq2, eq3, eq4, eq6, eq8, eq9], (r2, r3, x2, x3, x1, y1))

   2-element Vector{NTuple{6, Sym}}:
    (-179/79 + 62*sqrt(10)/79, -81/79 + 36*sqrt(10)/79, 853/79 - 186*sqrt(10)/79, 559/79 - 108*sqrt(10)/79, 8/5, 9/5)
    (11/9 - 2*sqrt(10)/9, 161/81 - 44*sqrt(10)/81, 1/3 + 2*sqrt(10)/3, -53/27 + 44*sqrt(10)/27, 8/5, 9/5)

2 組の解が得られるが，最初の解は r2 < r3 となり不適切解なので，2番目の解が適切。

name = ("r2", "r3", "x2", "x3", "x1", "y1")
j = 2
for i in 1:6
   println("$(name[i]) = $(res[j][i]) = $(res[j][i].evalf())")
end

   r2 = 11/9 - 2*sqrt(10)/9 = 0.519493853295916
   r3 = 161/81 - 44*sqrt(10)/81 = 0.269873863612238
   x2 = 1/3 + 2*sqrt(10)/3 = 2.44151844011225
   x3 = -53/27 + 44*sqrt(10)/27 = 3.19037840916328
   x1 = 8/5 = 1.60000000000000
   y1 = 9/5 = 1.80000000000000

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")
   r1 = 1
   (r2, r3, x2, x3, x1, y1) = (11/9 - 2*sqrt(10)/9, 161/81 - 44*sqrt(10)/81, 1/3 + 2*sqrt(10)/3, -53/27 + 44*sqrt(10)/27, 8/5, 9/5)
   println("r1 = $r1;  r2 = $r2;  r3 = $r3;  x2 = $x2;  x3 = $x3")
   plot()
   circle(r1, r1, r1, :red)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :brown)
   plot!([0, 4, 0, 0], [0, 0, 3, 0], color=:black, linewidth=0.25)
   if more
       point(r1, r1, "大円 r1", :red)
       point(x2, r2, "中円 r2", :blue)
       point(x3, r3, "小円 r3", :brown)
       point(0, 3, "  A", :red)
       point(x1, y1, "  B", :red)
       point(4, 0, "  C", :red)
       point(r1, 0, " D=r1", :red, :left, :bottom)
       point(0, r1, " E=r1", :red, :left, :bottom)
       point(x2, 0, "x2", :blue, :left, :bottom)
       point(x3, 0, "x3", :brown, :left, :bottom)
   end
end;

算額（その40）

算額（その40）

神奈川県川崎市多摩区 須賀神社祖師堂 文政6年
http://www.wasan.jp/kanagawa/suga.html

外円の中に正方形，正三角形，2 種類の小円が収まっている。各小円の径を求めよ。

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

using SymPy
@syms r1::positive, r2::positive, x1::positive;

上部にある円の半径: -1/2 + sqrt(2)/2

r1 = solve(sqrt(Sym(2))r1 + r1 - 1//2)[1]
r1 |> println

   -1/2 + sqrt(2)/2

左右にある円の半径と中心座標

eq4 = (x1 + r2/sqrt(Sym(2))) - x1 - (1-(x1 + r2/sqrt(Sym(2))))
eq5 = x1 - 1/sqrt(Sym(3)) - 2r2/sqrt(Sym(3))
res = solve([eq4, eq5], (r2, x1));

半径

res[r2] |> println

   -sqrt(3) - sqrt(6)/2 + 1 + 3*sqrt(2)/2

円の中心の x 座標

res[x1] |> println

   -2 - sqrt(2) + sqrt(3) + sqrt(6)

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)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

radian(degree) = degree / 180 * π;

function plotpolygon2(ox, oy, r, n; φ=90, col=:black, lty=:solid, lwd=1, fcol="")
   θ = range(0, 2π, length=n+1) .+ radian(φ)
   if fcol == ""
       plot!(r .* cos.(θ) .+ ox, r .* sin.(θ) .+ oy,
             linecolor=col, linestyle=lty, linewidth=lwd)
   else
       plot!(r .* cos.(θ) .+ ox, r .* sin.(θ) .+ oy,
             linecolor=col, linestyle=lty, linewidth=lwd,
             seriestype=:shape, fillcolor=fcol)
   end
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   circle(0, 0, 1, :black)
   plotpolygon2(0, 0, 1, 3, φ=150, lwd=0.25, col=:blue)
   plotpolygon2(0, 0, 1, 4, φ=90, lwd=0.25, col=:blue)
   r1 = (sqrt(2)-1)/2
   circle(0, 0.5+r1, r1)
   r2 = -sqrt(3) - sqrt(6)/2 + 1 + 3*sqrt(2)/2
   x1 = -2 - sqrt(2) + sqrt(3) + sqrt(6)
   circle(x1, 0, r2, :magenta)
   circle(-x1, 0, r2, :magenta)
   println("r1 = $r1, r2 = $r2")
   if more
       point(0, 0.5+r1, "0.5+r1 ", :red, :center)
       point(x1, 0, " x1", :red, :left)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
end;

 

 

算額（その39）

算額（その39）

香川県三豊市豊中町 本山寺 万延元年
http://www.wasan.jp/kagawa/honzanji1.html

外円の中に 4 種類の円が収まっている。青，白の円の径はそれぞれ 69 寸，15 寸である。外円の径を求めよ。



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

using SymPy
@syms r0::positive, r1::positive, r2::positive, r3::positive, r4::positive,
     x3::positive, y3::positive, x4::positive, y4::positive;

r1 = 69
r4 = 15
eq1 = x3^2 + (2r0-r1-y3)^2 - (r1+r3)^2
eq2 = x3^2 + (r0-y3)^2 - (r0-r3)^2
eq3 = x4^2 + (r0-y4)^2 - (r0-r4)^2
eq4 = x3^2 + (y3-r2)^2 - (r2+r3)^2
eq5 = x4^2 + (r2-y4)^2 - (r2+r4)^2
eq6 = (x3-x4)^2 + (y3-y4)^2 - (r3+r4)^2
eq7 = r1 + r2 - r0;

res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r0, r2, r3, x3, y3, x4, y4));

name = ["r0", "r2", "r3", "x3", "y3", "x4", "y4"]
for i in 1:7
   println("$(name[i]) = $(res[1][i])")
end

   r0 = 115
   r2 = 46
   r3 = 690/19
   x3 = 1380/19
   y3 = 1610/19
   x4 = 60
   y4 = 35

求める外円の径 r0 は 115 である。

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)
   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")
   r1 = 69
   r2 = 15
   (r0, r2, r3, x3, y3, x4, y4) = (115, 46, 690/19, 1380/19, 1610/19, 60, 35)
   println("r0 = $r0,  r1 = $r1,  r2 = $r2,  r3 = $r3,  r4 = $r4")
   println("x3 = $x3,  y3 = $y3,  x4 = $x4,  y4 = $y4")
   plot()
   circle(0, r2, r2, :red)
   circle(0, 2r0-r1, r1, :blue)
   circle(0, r0, r0, :black)
   circle(x3, y3, r3, :brown)
   circle(-x3, y3, r3, :brown)
   circle(x4, y4, r4, :magenta)
   circle(-x4, y4, r4, :magenta)
   if more
       point(0, 2r0-r1, "2r0-r1 ", :blue, :right)
       point(0, r0, "r0 ", :black, :right)
       point(0, r2, "r2 ", :red, :right)
       point(x3, y3, "(x3,y3),r3", :brown, :center)
       point(x4, y4, "(x4,y4),r4", :magenta, :center)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
end;

算額（その38）

算額（その38）

新潟県糸魚川市 天津神社 寛政12年
http://www.wasan.jp/niigata/tensin.html

図のように，外円の中に 5 種類の径（木，火，土，金，水）を持つ円が収まっている。土円，金円の径はそれぞれ 9分9釐（厘），一寸9分9釐（厘）である。水円の径を求めよ。

図のように記号を定め，方程式を解く。水円の径は r3 である。（プログラムでは，すべて直径ではなく半径を使っている。）

using SymPy
@syms　r1::positive, r2::positive, y2::positive, r0::positive, y3::positive, r3::positive, r4::positive, x4::positive, y4::positive, r5::positive, x5::positive, y5::positive;

r1 = 99
r2 = 198
y3 = r3
x5 = r5
x4 = r4
eq1 = x4^2+(2r0-r1-y4)^2-(r1+r4)^2
eq2 = x4^2+(y4-y2)^2-(r2+r4)^2
eq3 = x5^2+(y2-y5)^2-(r2+r5)^2
eq4 = x5^2+(y5-y3)^2-(r3+r5)^2
eq5 = x4^2+(y4-r0)^2-(r0-r4)^2b
eq6 = x5^2+(r0-y5)^2-(r0-r5)^2
eq7 = (x5 - x4)^2 + (y4 - y5)^2 - (r4 + r5)^2;

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (y2, r0, r3, r4, y4, r5, y5))

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)
   y2, r0, r3, r4, y4, r5, y5 = u
   return [
       r4^2 - (r4 + 99)^2 + (2*r0 - y4 - 99)^2
       r4^2 - (r4 + 198)^2 + (-y2 + y4)^2
       r5^2 - (r5 + 198)^2 + (y2 - y5)^2
       r5^2 + (-r3 + y5)^2 - (r3 + r5)^2
       r4^2 + (-r0 + y4)^2 - (r0 - r4)^2
       r5^2 - (r0 - r5)^2 + (r0 - y5)^2
       (-r4 + r5)^2 - (r4 + r5)^2 + (y4 - y5)^2
 ];
end;

iniv = [2260., 1521, 580, 326, 2670, 752, 1679];

res = nls(h, ini=iniv)

y2, r0, r3, r4, y4, r5, y5 = res[1]

   7-element Vector{Float64}:
    2260.2616848711004
    1521.243310737389
     580.0010348927411
     326.36779179808775
    2670.68298737871
     752.3774162147455
    1679.6194053038619

求める水円の径は 580.0010348927411 である。

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)
   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")
   r1 = 99
   r2 = 198
   y2, r0, r3, r4, y4, r5, y5 = [2260.2616848711004, 1521.243310737389, 580.0010348927411, 326.36779179808775, 2670.68298737871, 752.3774162147455, 1679.6194053038619]
   y3 = r3
   x4 = r4
   x5 = r5
   println("r0 = $r0,  r1 = $r1,  r2 = $r2,  r3 = $r3,  r4 = $r4,  r5 = $r5")
   println("y2 = $y2,  y3 = $y3,  x4 = $x4,  y4 = $y4,  x5 = $x5,  y5 = $y5")
   plot()
   circle(0, r0, r0, :green)
   circle(0, 2r0-r1, r1)
   circle(0, y2, r2, :blue)
   circle(0, y3, r3, :brown)
   circle(x4, y4, r4, :magenta)
   circle(-x4, y4, r4, :magenta)
   circle(x5, y5, r5, :black)
   circle(-x5, y5, r5, :black)
   if more
       point(0, r0, "r0 ", :green, :right)
       point(0, 2r0-r1, "r1 ", :red, :right)
       point(0, y2, "y2 ", :blue, :right)
       point(0, y3, "r3 ", :brown, :right)
       point(x4, y4, "(x4,y4)", :magenta, :center)
       point(x5, y5, "(x5,y5)", :black, :center)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
end;

 

算額（その37）

算額（その37）

新潟県三島郡 与板八幡宮 寛政12年
新潟県柏崎市与板 与板三柱神社 であろう
http://www.wasan.jp/niigata/yoitahatiman1.html

「長」，「平」，「高」の 3 辺を持つ直方体において，

• 「長」，「平」 の差は 1 寸
• 「斜」は 13 寸

直方体の体積が最も大きくなるのは，「長」がいくつのときか。

直方体の底面の対角線を L として，以下の式を立てる。

using SymPy
@syms 長::positive, 平::positive, 高::positive, 斜::positive, L::positive;

斜 = 13
平 = 長 - 1
eq2 = 長^2 + 平^2 - L^2;
eq3 = L^2 + 高^2 - 斜^2;

res = solve([eq2, eq3], (L, 高))

   1-element Vector{Tuple{Sym, Sym}}:
    (sqrt(2*長^2 - 2*長 + 1), sqrt(2)*sqrt(-長^2 + 長 + 84))

volume = 長 * 平 * res[1][2]
volume |> println

   sqrt(2)*長*(長 - 1)*sqrt(-長^2 + 長 + 84)

diff(volume) |> solve |> println

   Sym[1/2, 8]

1/2 は不適切解なので，正解は 8 （寸）である。

ちなみに，高さは

sqrt(2)*sqrt(-長^2 + 長 + 84)(長 => 8).evalf()

   7.48331477354788

 

算額（その36）

算額（その36）

新潟県三島郡 与板八幡宮 寛政12年
新潟県柏崎市与板 与板三柱神社 であろう
http://www.wasan.jp/niigata/yoitahatiman2.html

直線と互いに接する大円と小円があり，更に小円と大円と直線に接する連続する小円がある。大円の径が225寸，一番小さい小円の径が4寸のとき，赤で示した面積が最大になるときの一番大きい小円の径はいくつになるか。

答えを先に述べておく。

算額によれば，径 4 寸から径 100 寸の最大の小円まで全部で 5 個とあるが，本当は 7 個ではないか?

この問題は，算額（その35）の後半に示したプログラムで解くことができる。
大円の半径を R，小さい方から 2 つの小円の半径と中心の x 座標を r1, x1, r2, x2 と置き，以下の方程式を立てる。

using SymPy
@syms R::positive, r1::positive, r2::positive, x1::positive, x2::positive;

R = 225
eq1 = x1^2 + (R - r1)^2 - (R + r1)^2;
eq2 = x2^2 + (R - r2)^2 - (R + r2)^2;
eq3 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2;

solve(eq1, x1) |> println

   Sym[2*sqrt(R)*sqrt(r1)]

res = solve([eq1, eq2, eq3], (r1, x1));
res[2][1] |> simplify |> println
res[2][2] |> simplify |> println

   x2^2*(2*sqrt(R)*sqrt(r2) + R + r2)/(4*(R - r2)^2)
   x2*(sqrt(R)*sqrt(r2) + R)/(R - r2)

SymPy はここまで。JupyterLab などを使っている場合は，ここでカーネルをリスタートする。

次に大きい円を求める関数を定義する。

function nextcircle(r2, x2; R = 225)
   r = x2^2*(2*sqrt(R)*sqrt(r2) + R + r2)/(4*(R - r2)^2)
   x = x2*(sqrt(R)*sqrt(r2) + R)/(R - r2)
   return r, x
end;

この漸化式によれば，半径が 4 寸からスタートして 7 番目の小円の半径がほぼ 100 になることが確認できる。

算額によれば，径 4 寸から径 100 寸の最大の小円まで全部で 5 個とあるが，本当は 7 個ではないか?

R = 225
r = 4 # r1
x = 2*sqrt(R)*sqrt(r) # x1

for i = 2:7
   r, x = nextcircle(r, x, R = 225)
   println("$i, $r, $x")
end

   2, 5.325443786982248, 69.23076923076923
   3, 7.438016528925618, 81.81818181818181
   4, 11.111111111111112, 99.99999999999999
   5, 18.3673469387755, 128.57142857142856
   6, 35.99999999999999, 179.99999999999997
   7, 99.99999999999996, 299.99999999999994

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)
   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")
   R = 225
   plot()
   circle(0, R, R, :green)
   r = 4 # r1
   x = 2*sqrt(R)*sqrt(r) # x1
   circle(x, r, r)
   println("1, r = $r, x = $x")

   for i = 2:7
       r, x = nextcircle(r, x, R=R)
       println("$i, r = $r, x = $x")
       circle(x, r, r)
   end
   hline!([0], color=:black, linewidth=0.25)
   if more
       point(0, 0, "0 ", :black, :right)
       point(0, r3, "r3 ", :black, :right)
       point(0, r3+r1, "r3+r1 ", :blue, :right)
       point(0, 1-r2, "1-r2 ", :red, :right)
       point(x4, r4, "(x4,r4)", :magenta, :center)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
end;

draw(false)

   1, r = 4, x = 60.0
   2, r = 5.325443786982248, x = 69.23076923076923
   3, r = 7.438016528925618, x = 81.81818181818181
   4, r = 11.111111111111112, x = 99.99999999999999
   5, r = 18.3673469387755, x = 128.57142857142856
   6, r = 35.99999999999999, x = 179.99999999999997
   7, r = 99.99999999999996, x = 299.99999999999994

 

 

算額（その35）

算額（その35）

新潟県三島郡 与板八幡宮 寛政12年
新潟県柏崎市与板 与板三柱神社 であろう
http://www.wasan.jp/niigata/yoitahatiman1.html

直線と互いに接する大円と小円があり，更に小円と大円と直線に接する連続する小円がある。大円の径が10寸，一番小さい小円の径が1分のとき，小円の数はいくつか。

大円の半径を R，大きい方から 2 つの小円の半径と中心の x 座標を r1, x1, r2, x2 と置き，以下の方程式を立てる。

using SymPy
@syms R::positive, r1::positive, r2::positive, x1::positive, x2::positive;

R = 1
eq1 = x1^2 + (R - r1)^2 - (R + r1)^2;
eq2 = x2^2 + (R - r2)^2 - (R + r2)^2;
eq3 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2;

一番大きい小円の中心の x 座標は，eq1 を x1 について解けば求まる。

solve(eq1, x1)

   2*sqrt(r1)

eq1, eq2, eq3 を r2, x2 について解けば，ある小円の次に小さい小円の半径と中心の x 座標を求める漸化式が得られる。なお，2番めの解は不適切解である。

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

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

res[1][1] |> simplify |> println # r2

   x1^2*(-2*sqrt(r1) + r1 + 1)/(4*(r1 - 1)^2)

res[1][2] |> simplify |> println # x2

   x1*(sqrt(r1) - 1)/(r1 - 1)

小円 (r1, x1) から次の小円 (r2, x2) を得る関数を書く。

function nextcircle(r1, x1; R = 1)
   r2 = x1^2*(-2*sqrt(r1) + r1 + 1)/(4*(r1 - 1)^2)
   x2 = x1*(sqrt(r1) - 1)/(r1 - 1)
   return r2, x2
end;

この関数を使い，初期値を 1 にするとエラーになるが，1 に極めて近い値を設定することで 10 番目の小円の半径が大円のほぼ 1/100 = 0.01 になることがわかる。

r = 0.99999 # r1
x = 2sqrt(r) # x1

for i = 2:10
   r, x = nextcircle(r, x)
   println("$i, $r, $x")
end

   2, 0.2499986308990838, 0.999997499986309
   3, 0.11111075838251418, 0.6666656084800361
   4, 0.062499851192534595, 0.499999404769784
   5, 0.03999992381055049, 0.399999619052571
   6, 0.027777733686650655, 0.3333330687864656
   7, 0.020408135499460415, 0.28571409135329967
   8, 0.015624981399050211, 0.24999985119235743
   9, 0.012345665948302833, 0.22222210464580552
   10, 0.009999990476312013, 0.19999990476309745

図示することで，正しく描画されることを確認する。

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)
   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")
   R = 1
   plot()
   circle(0, R, R, :green)
   r = 0.99999 # r1
   x = 2sqrt(r) # x1
   circle(x, r, r)

   for i = 2:11
       r, x = nextcircle(r, x, R=R)
       println("$i, r = $r, x = $x")
       circle(x, r, r)
   end
   hline!([0], color=:black, linewidth=0.25)
   if more
       point(0, 0, "0 ", :black, :right)
       point(0, r3, "r3 ", :black, :right)
       point(0, r3+r1, "r3+r1 ", :blue, :right)
       point(0, 1-r2, "1-r2 ", :red, :right)
       point(x4, r4, "(x4,r4)", :magenta, :center)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
end;

draw(false)

   2, r = 0.2499986308990838, x = 0.999997499986309
   3, r = 0.11111075838251418, x = 0.6666656084800361
   4, r = 0.062499851192534595, x = 0.499999404769784
   5, r = 0.03999992381055049, x = 0.399999619052571
   6, r = 0.027777733686650655, x = 0.3333330687864656
   7, r = 0.020408135499460415, x = 0.28571409135329967
   8, r = 0.015624981399050211, x = 0.24999985119235743
   9, r = 0.012345665948302833, x = 0.22222210464580552
   10, r = 0.009999990476312013, x = 0.19999990476309745
   11, r = 0.008264455654629146, x = 0.18181810310999452

-----------
先の方程式を，小円 r2, x2 から次に大きい小円 r1, x1 を求める漸化式を作るように解くと以下のようになる。

using SymPy
@syms R::positive, r1::positive, r2::positive, x1::positive, x2::positive;

R = 1
eq1 = x1^2 + (R - r1)^2 - (R + r1)^2;
eq2 = x2^2 + (R - r2)^2 - (R + r2)^2;
eq3 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2;

solve(eq2, x2)

   2*sqrt(r2)

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

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

function nextcircle(r2, x2; R = 1)
   r = -x2*(-2*sqrt(r2)*x2/(r2 - 1) - x2 - 2*x2/(r2 - 1))/(4*(r2 - 1))
   x = -sqrt(r2)*x2/(r2 - 1) - x2/(r2 - 1)
   return r, x
end;

この漸化式によれば，半径が 0.01 からスタートして 10 番目の小円の半径がほぼ 1 になることが確認できる。

r = 0.01 # r1
x = 2sqrt(r) # x1

for i = 2:10
   r, x = nextcircle(r, x)
   println("$i, $r, $x")
end

   2, 0.012345679012345682, 0.22222222222222227
   3, 0.01562500000000001, 0.25000000000000006
   4, 0.02040816326530613, 0.28571428571428575
   5, 0.02777777777777779, 0.3333333333333334
   6, 0.04000000000000003, 0.40000000000000013
   7, 0.06250000000000006, 0.5000000000000002
   8, 0.1111111111111112, 0.666666666666667
   9, 0.25000000000000033, 1.0000000000000007
   10, 1.0000000000000027, 2.0000000000000027

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)
   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")
   R = 1
   plot()
   circle(0, R, R, :green)
   r = 0.01 # r1
   x = 2sqrt(r) # x1
   circle(x, r, r)
   println("1, r = $r, x = $x")

   for i = 2:10
       r, x = nextcircle(r, x, R=R)
       println("$i, r = $r, x = $x")
       circle(x, r, r)
   end
   hline!([0], color=:black, linewidth=0.25)
   if more
       point(0, 0, "0 ", :black, :right)
       point(0, r3, "r3 ", :black, :right)
       point(0, r3+r1, "r3+r1 ", :blue, :right)
       point(0, 1-r2, "1-r2 ", :red, :right)
       point(x4, r4, "(x4,r4)", :magenta, :center)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
end;

draw(false)

   1, r = 0.01, x = 0.2
   2, r = 0.012345679012345682, x = 0.22222222222222227
   3, r = 0.01562500000000001, x = 0.25000000000000006
   4, r = 0.02040816326530613, x = 0.28571428571428575
   5, r = 0.02777777777777779, x = 0.3333333333333334
   6, r = 0.04000000000000003, x = 0.40000000000000013
   7, r = 0.06250000000000006, x = 0.5000000000000002
   8, r = 0.1111111111111112, x = 0.666666666666667
   9, r = 0.25000000000000033, x = 1.0000000000000007
   10, r = 1.0000000000000027, x = 2.0000000000000027

二番目に小さい小円の半径が 0.012345679012345682 というのは，意味深だなあ。

 

算額（その34）

算額（その34）

三重県四日市市 川島神明神社 天保15年
http://www.wasan.jp/mie/sinmei1.html

直径121寸8分の大円と直線に接する円，およびこの円と大円に接する小さな円がある。小さな円の直径が1分になるのは何番目の円か。

以下のプログラムでは半径1218分の円と考える。求めるのは小円の半径が1分の円である。

答えを最初に書いておく。

16 番目の小円の半径が 0.9918566461328249 である。

まず大円と直線に接する円の直径と中心の x 座標の漸化式を求める。

using SymPy
@syms xi, xi1, ri, ri1, R;

eq1 = (xi1 - xi)^2 + (ri - ri1)^2 - (ri + ri1)^2
eq2 = (R - xi)^2 + (R - ri)^2 - (R + ri)^2
eq3 = (R - xi1)^2 + (R - ri1)^2 - (R + ri1)^2;

res = solve([eq1, eq2, eq3], (xi1, ri1))

   2-element Vector{Tuple{Sym, Sym}}:
    ((R^2 - xi^2 + (-4*R + 4*ri)*(-sqrt(R*ri)*(R - xi)^2/(2*(R^2 - 2*R*ri + ri^2)) + (R + ri)*(R - xi)^2/(4*(R - ri)^2)))/(2*(R - xi)), -sqrt(R*ri)*(R - xi)^2/(2*(R^2 - 2*R*ri + ri^2)) + (R + ri)*(R - xi)^2/(4*(R - ri)^2))
    ((R^2 - xi^2 + (-4*R + 4*ri)*(sqrt(R*ri)*(R - xi)^2/(2*(R^2 - 2*R*ri + ri^2)) + (R + ri)*(R - xi)^2/(4*(R - ri)^2)))/(2*(R - xi)), sqrt(R*ri)*(R - xi)^2/(2*(R^2 - 2*R*ri + ri^2)) + (R + ri)*(R - xi)^2/(4*(R - ri)^2))

res[2] は不適切解である。

function f(x, r; R=1218) # 次の円のx座標と半径
   return (-R*r + R*x + R*sqrt(R*r) - x*sqrt(R*r))/(R - r), (R - x)^2*(R + r - 2*sqrt(R*r))/(4*(R - r)^2)
end;

function nextcircle(n) # 後で使うために配列に収める
   xa = zeros(n)
   ra = zeros(n)
   x = 0
   r = 609/2
   xa[1] = x
   ra[1] = r
   for i = 2:n
       x, r = f(x, r)
       #println("$i: x = $x, r = $r")
       xa[i] = x
       ra[i] = r
   end
   return xa, ra
end;

次に，大円と前述の円に接する小円の半径と中心座標を求める漸化式を作る。

using SymPy
@syms R::positive, r1::positive, r2::positive, r3::positive, x2::positive, x3::positive, y3::positive;
@syms X2::positive, Y2::positive, ra2::positive;
eq1 = (R - x2)^2 + (R - r2)^2 - (R + r2)^2
eq2 = (x3 - x2)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq3 = (R - x3)^2 + (R - r3)^2 - (R + r3)^2
eq4 = (X2 - x2)^2 + (Y2 - r2)^2 - (ra2 + r2)^2
eq5 = (x3 - X2)^2 + (Y2 - r3)^2 - (ra2 + r3)^2
eq6 = (R - X2)^2 + (R - Y2)^2 - (R + ra2)^2;

res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (X2, Y2, ra2))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    ((R^2*r2 - R^2*r3 - R*x2^2 + R*x3^2 - r2*x3^2 + r3*x2^2)/(2*(R*r2 - R*r3 - R*x2 + R*x3 - r2*x3 + r3*x2)), -(-R^4*r2^2 + 2*R^4*r2*r3 - R^4*r3^2 + R^4*x2^2 - 2*R^4*x2*x3 + R^4*x3^2 + 4*R^3*r2^2*x3 - 4*R^3*r2*r3*x2 - 4*R^3*r2*r3*x3 + 2*R^3*r2*x2^2 - 4*R^3*r2*x2*x3 + 2*R^3*r2*x3^2 + 4*R^3*r3^2*x2 + 2*R^3*r3*x2^2 - 4*R^3*r3*x2*x3 + 2*R^3*r3*x3^2 - 2*R^3*x2^3 + 2*R^3*x2^2*x3 + 2*R^3*x2*x3^2 - 2*R^3*x3^3 - 6*R^2*r2^2*x3^2 + 2*R^2*r2*r3*x2^2 + 8*R^2*r2*r3*x2*x3 + 2*R^2*r2*r3*x3^2 - 4*R^2*r2*x2^2*x3 + 8*R^2*r2*x2*x3^2 - 4*R^2*r2*x3^3 - 6*R^2*r3^2*x2^2 - 4*R^2*r3*x2^3 + 8*R^2*r3*x2^2*x3 - 4*R^2*r3*x2*x3^2 + 6*R^2*x2^3*x3 - 12*R^2*x2^2*x3^2 + 6*R^2*x2*x3^3 + 4*R*r2^2*x3^3 - 4*R*r2*r3*x2^2*x3 - 4*R*r2*r3*x2*x3^2 + 2*R*r2*x2^2*x3^2 - 4*R*r2*x2*x3^3 + 2*R*r2*x3^4 + 4*R*r3^2*x2^3 + 2*R*r3*x2^4 - 4*R*r3*x2^3*x3 + 2*R*r3*x2^2*x3^2 - 2*R*x2^4*x3 + 2*R*x2^3*x3^2 + 2*R*x2^2*x3^3 - 2*R*x2*x3^4 - r2^2*x3^4 + 2*r2*r3*x2^2*x3^2 - r3^2*x2^4 + x2^4*x3^2 - 2*x2^3*x3^3 + x2^2*x3^4)/(4*(-R + x2)*(-R + x3)*(x2 - x3)*(R*r2 - R*r3 - R*x2 + R*x3 - r2*x3 + r3*x2)), -(R^4*r2^2 - 2*R^4*r2*r3 + R^4*r3^2 + R^4*x2^2 - 2*R^4*x2*x3 + R^4*x3^2 - 4*R^3*r2^2*x3 + 4*R^3*r2*r3*x2 + 4*R^3*r2*r3*x3 - 2*R^3*r2*x2^2 + 4*R^3*r2*x2*x3 - 2*R^3*r2*x3^2 - 4*R^3*r3^2*x2 - 2*R^3*r3*x2^2 + 4*R^3*r3*x2*x3 - 2*R^3*r3*x3^2 - 2*R^3*x2^3 + 2*R^3*x2^2*x3 + 2*R^3*x2*x3^2 - 2*R^3*x3^3 + 6*R^2*r2^2*x3^2 - 2*R^2*r2*r3*x2^2 - 8*R^2*r2*r3*x2*x3 - 2*R^2*r2*r3*x3^2 + 4*R^2*r2*x2^2*x3 - 8*R^2*r2*x2*x3^2 + 4*R^2*r2*x3^3 + 6*R^2*r3^2*x2^2 + 4*R^2*r3*x2^3 - 8*R^2*r3*x2^2*x3 + 4*R^2*r3*x2*x3^2 + 2*R^2*x2^4 - 2*R^2*x2^3*x3 - 2*R^2*x2*x3^3 + 2*R^2*x3^4 - 4*R*r2^2*x3^3 + 4*R*r2*r3*x2^2*x3 + 4*R*r2*r3*x2*x3^2 - 2*R*r2*x2^2*x3^2 + 4*R*r2*x2*x3^3 - 2*R*r2*x3^4 - 4*R*r3^2*x2^3 - 2*R*r3*x2^4 + 4*R*r3*x2^3*x3 - 2*R*r3*x2^2*x3^2 - 2*R*x2^4*x3 + 2*R*x2^3*x3^2 + 2*R*x2^2*x3^3 - 2*R*x2*x3^4 + r2^2*x3^4 - 2*r2*r3*x2^2*x3^2 + r3^2*x2^4 + x2^4*x3^2 - 2*x2^3*x3^3 + x2^2*x3^4)/(4*(-R + x2)*(-R + x3)*(x2 - x3)*(R*r2 - R*r3 - R*x2 + R*x3 - r2*x3 + r3*x2)))

次の小円を求める関数を定義する。

function nextcircle2(x, r, i; R=1218)
   X = (R^2*r[i] - R^2*r[i+1] - R*x[i]^2 + R*x[i+1]^2 - r[i]*x[i+1]^2 + r[i+1]*x[i]^2)/(2*(R*r[i] - R*r[i+1] - R*x[i] + R*x[i+1] - r[i]*x[i+1] + r[i+1]*x[i]))
   Y = (R^4*r[i]^2/4 - R^4*r[i]*r[i+1]/2 + R^4*r[i+1]^2/4 - R^4*x[i]^2/4 + R^4*x[i]*x[i+1]/2 - R^4*x[i+1]^2/4 - R^3*r[i]^2*x[i+1] + R^3*r[i]*r[i+1]*x[i] + R^3*r[i]*r[i+1]*x[i+1] - R^3*r[i]*x[i]^2/2 + R^3*r[i]*x[i]*x[i+1] - R^3*r[i]*x[i+1]^2/2 - R^3*r[i+1]^2*x[i] - R^3*r[i+1]*x[i]^2/2 + R^3*r[i+1]*x[i]*x[i+1] - R^3*r[i+1]*x[i+1]^2/2 + R^3*x[i]^3/2 - R^3*x[i]^2*x[i+1]/2 - R^3*x[i]*x[i+1]^2/2 + R^3*x[i+1]^3/2 + 3*R^2*r[i]^2*x[i+1]^2/2 - R^2*r[i]*r[i+1]*x[i]^2/2 - 2*R^2*r[i]*r[i+1]*x[i]*x[i+1] - R^2*r[i]*r[i+1]*x[i+1]^2/2 + R^2*r[i]*x[i]^2*x[i+1] - 2*R^2*r[i]*x[i]*x[i+1]^2 + R^2*r[i]*x[i+1]^3 + 3*R^2*r[i+1]^2*x[i]^2/2 + R^2*r[i+1]*x[i]^3 - 2*R^2*r[i+1]*x[i]^2*x[i+1] + R^2*r[i+1]*x[i]*x[i+1]^2 - 3*R^2*x[i]^3*x[i+1]/2 + 3*R^2*x[i]^2*x[i+1]^2 - 3*R^2*x[i]*x[i+1]^3/2 - R*r[i]^2*x[i+1]^3 + R*r[i]*r[i+1]*x[i]^2*x[i+1] + R*r[i]*r[i+1]*x[i]*x[i+1]^2 - R*r[i]*x[i]^2*x[i+1]^2/2 + R*r[i]*x[i]*x[i+1]^3 - R*r[i]*x[i+1]^4/2 - R*r[i+1]^2*x[i]^3 - R*r[i+1]*x[i]^4/2 + R*r[i+1]*x[i]^3*x[i+1] - R*r[i+1]*x[i]^2*x[i+1]^2/2 + R*x[i]^4*x[i+1]/2 - R*x[i]^3*x[i+1]^2/2 - R*x[i]^2*x[i+1]^3/2 + R*x[i]*x[i+1]^4/2 + r[i]^2*x[i+1]^4/4 - r[i]*r[i+1]*x[i]^2*x[i+1]^2/2 + r[i+1]^2*x[i]^4/4 - x[i]^4*x[i+1]^2/4 + x[i]^3*x[i+1]^3/2 - x[i]^2*x[i+1]^4/4)/((R - x[i])*(R - x[i+1])*(x[i] - x[i+1])*(R*r[i] - R*r[i+1] - R*x[i] + R*x[i+1] - r[i]*x[i+1] + r[i+1]*x[i]))
   rr = (-R^4*r[i]^2/4 + R^4*r[i]*r[i+1]/2 - R^4*r[i+1]^2/4 - R^4*x[i]^2/4 + R^4*x[i]*x[i+1]/2 - R^4*x[i+1]^2/4 + R^3*r[i]^2*x[i+1] - R^3*r[i]*r[i+1]*x[i] - R^3*r[i]*r[i+1]*x[i+1] + R^3*r[i]*x[i]^2/2 - R^3*r[i]*x[i]*x[i+1] + R^3*r[i]*x[i+1]^2/2 + R^3*r[i+1]^2*x[i] + R^3*r[i+1]*x[i]^2/2 - R^3*r[i+1]*x[i]*x[i+1] + R^3*r[i+1]*x[i+1]^2/2 + R^3*x[i]^3/2 - R^3*x[i]^2*x[i+1]/2 - R^3*x[i]*x[i+1]^2/2 + R^3*x[i+1]^3/2 - 3*R^2*r[i]^2*x[i+1]^2/2 + R^2*r[i]*r[i+1]*x[i]^2/2 + 2*R^2*r[i]*r[i+1]*x[i]*x[i+1] + R^2*r[i]*r[i+1]*x[i+1]^2/2 - R^2*r[i]*x[i]^2*x[i+1] + 2*R^2*r[i]*x[i]*x[i+1]^2 - R^2*r[i]*x[i+1]^3 - 3*R^2*r[i+1]^2*x[i]^2/2 - R^2*r[i+1]*x[i]^3 + 2*R^2*r[i+1]*x[i]^2*x[i+1] - R^2*r[i+1]*x[i]*x[i+1]^2 - R^2*x[i]^4/2 + R^2*x[i]^3*x[i+1]/2 + R^2*x[i]*x[i+1]^3/2 - R^2*x[i+1]^4/2 + R*r[i]^2*x[i+1]^3 - R*r[i]*r[i+1]*x[i]^2*x[i+1] - R*r[i]*r[i+1]*x[i]*x[i+1]^2 + R*r[i]*x[i]^2*x[i+1]^2/2 - R*r[i]*x[i]*x[i+1]^3 + R*r[i]*x[i+1]^4/2 + R*r[i+1]^2*x[i]^3 + R*r[i+1]*x[i]^4/2 - R*r[i+1]*x[i]^3*x[i+1] + R*r[i+1]*x[i]^2*x[i+1]^2/2 + R*x[i]^4*x[i+1]/2 - R*x[i]^3*x[i+1]^2/2 - R*x[i]^2*x[i+1]^3/2 + R*x[i]*x[i+1]^4/2 - r[i]^2*x[i+1]^4/4 + r[i]*r[i+1]*x[i]^2*x[i+1]^2/2 - r[i+1]^2*x[i]^4/4 - x[i]^4*x[i+1]^2/4 + x[i]^3*x[i+1]^3/2 - x[i]^2*x[i+1]^4/4)/((R - x[i])*(R - x[i+1])*(x[i] - x[i+1])*(R*r[i] - R*r[i+1] - R*x[i] + R*x[i+1] - r[i]*x[i+1] + r[i+1]*x[i]))
   return X, Y, rr
end;

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)
   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")
   R = 1218
   (r1, r2, r3, x2, x3, y3) = (609/2, 406/3, 87/2, 406, 348, 609/2)
   println("r1 = $r1;  r2 = $r2;  r3 = $r3")
   println("x2= $x2;  x3 = $x3;  y3 = $y3")
   plot()
   circle(R, R, R, :green, beginangle=180, endangle=270)
   circle(0, r1, r1, :magenta)
   circle(x2, r2, r2, :blue)
   circle(x3, y3, r3, :red)
   xa, ra = nextcircle(20)
   println("デカルトの円定理から，小円の半径のみを求める")
   kR = 1/R
   ka = 1 ./ ra
   for i = 1:19
       k4 = kR + ka[i] + ka[i+1] + 2sqrt(kR*ka[i] + ka[i]*ka[i+1] + ka[i+1]*kR)
       r = 1/k4
       println("$i, $r")
   end
       
   for i = 1:20
       circle(xa[i], ra[i], ra[i], :blue)
   end
   println("小円の半径，中心座標を求める")
   for i = 1:19
       X, Y, rr = nextcircle2(xa, ra, i)
       circle(X, Y, rr)
       println("i = $i, x = $X, y = $Y, r = $rr")
   end
   hline!([0], color=:black, linewidth=0.25)
end;

draw(false)

   r1 = 304.5;  r2 = 135.33333333333334;  r3 = 43.5
   x2= 406;  x3 = 348;  y3 = 304.5
   デカルトの円定理から，小円の半径のみを求める
   1, 43.5
   2, 23.423076923076916
   3, 14.499999999999996
   4, 9.822580645161292
   5, 7.0813953488372094
   6, 5.342105263157894
   7, 4.171232876712329
   8, 3.3461538461538454
   9, 2.7432432432432425
   10, 2.2894736842105234
   11, 1.9394904458598707
   12, 1.663934426229509
   13, 1.4431279620853108
   14, 1.2634854771784254
   15, 1.1153846153846176
   16, 0.9918566775244327
   17, 0.8877551020408169
   18, 0.7992125984251943
   19, 0.7232779097387141
   小円の半径，中心座標を求める
   i = 1, x = 348.0, y = 304.49999999999994, r = 43.499999999999986
   i = 2, x = 562.1538461538463, y = 163.96153846153766, r = 23.423076923075435
   i = 3, x = 695.9999999999999, y = 101.50000000001576, r = 14.499999999987837
   i = 4, x = 785.8064516129033, y = 68.75806451612152, r = 9.82258064507806
   i = 5, x = 849.7674418604653, y = 49.56976744177949, r = 7.08139534902618
   i = 6, x = 897.4736842105261, y = 37.39473684132309, r = 5.342105263101324
   i = 7, x = 934.3561643835608, y = 29.198630137693904, r = 4.171232877205791
   i = 8, x = 963.6923076923063, y = 23.423076928931703, r = 3.346153845145927
   i = 9, x = 987.5675675675733, y = 19.202702699235346, r = 2.7432432427479063
   i = 10, x = 1007.3684210526319, y = 16.02631580210383, r = 2.2894736816354584
   i = 11, x = 1024.0509554140167, y = 13.576433082760703, r = 1.9394904363737546
   i = 12, x = 1038.2950819672076, y = 11.6475410215815, r = 1.663934370306387
   i = 13, x = 1050.5971563981027, y = 10.101895600212682, r = 1.443127868224802
   i = 14, x = 1061.3278008298776, y = 8.844398562571168, r = 1.2634854854616941
   i = 15, x = 1070.7692307692128, y = 7.807692179665674, r = 1.1153844774220047
   i = 16, x = 1079.140065146608, y = 6.942996352466482, r = 0.9918566461328249
   i = 17, x = 1086.6122448979554, y = 6.214286399607383, r = 0.8877555400354336
   i = 18, x = 1093.322834645682, y = 5.594488988858351, r = 0.7992131560152675
   i = 19, x = 1099.3824228028282, y = 5.062945169849251, r = 0.7232791103597648

16 番目の小円の半径が 0.9918566461328249 である。

 

 

算額（その33）

算額（その33）

岩手県一関市　八幡神社  天保9年(1838)　
http://www.wasan.jp/iwate/itinosekiyahata2.html

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

算額の図では中央にある円と上下端にある円は違う色で塗られているが，これは同じ径を持つ円でなければ計算が合わない。

using SymPy
@syms r1::positive, r2::positive, r3::positive, r4::positive, x4::positive;

eq1 = x4^2 + r4^2 - (r3 + r4)^2;
eq2 = x4^2 + (r4 - r1 - r3)^2 - (r1 + r4)^2;
eq3 = r3 + 2r1 + 2r2 - 1;
eq4 = 1 - 2r4 - r3;
eq5 = r2 - r3;

solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, r3, r4, x4))

   1-element Vector{NTuple{5, Sym}}:
    (1/8, 1/4, 1/4, 3/8, 1/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)
   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")
   r1, r2, r3, r4, x4 = (1/8, 1/4, 1/4, 3/8, 1/2)
   println("r1 = $r1;  r2 = $r2;  r3 = $r3;  r4 = $r4")
   println("x4= $x4")
   plot()
   circle(0, 0, r3, :green)
   circle(0, r1+r3, r1, :blue)
   circle(0, -r1-r3, r1, :blue)
   circle(0, r3 + 2r1 + r2, r2, :red)
   circle(0, -r3 - 2r1 - r2, r2, :red)
   circle(x4, r4, r4, :magenta)
   circle(x4, -r4, r4, :magenta)
   circle(-x4, r4, r4, :magenta)
   circle(-x4, -r4, r4, :magenta)
   if more
       point(0, 0, "0 ", :black, :right)
       point(0, r3, "r3 ", :black, :right)
       point(0, r3+r1, "r3+r1 ", :blue, :right)
       point(0, 1-r2, "1-r2 ", :red, :right)
       point(x4, r4, "(x4,r4)", :magenta, :center)
       vline!([0], color=:black, linewidth=0.25)
       hline!([0], color=:black, linewidth=0.25)
   end
   circle(0, 0, 1, :black)
end;

算額（その32）

算額（その32）

岩手県一関市花泉町 大門神社・大門観世音菩薩 明治 33 年
http://www.wasan.jp/iwate/daimon5.html

図のように，外円の中に 3 種類の円が入っている。そのうちの 3 つは直線に接している。

外円，小円，中円，大円の半径と中心座標を {1, (0, y3)}, {r1, (0, r1)}, {r2, (x2, r2)}，{r3, (0, 2r1 + r3)} と置き，4 本の方程式を立てる。
しかし，これでは未知数 5 個に対して条件が一つ足りない。そこで，r1 を何らかの値に設定する。

r2, r3 も比較的きれいな分数式になるのは r1 = 7//64, 12//64, 15//64, 16//64 の 4 通りである。

算額的に最もきれいなのは次の解（r1 = 1/4）であろう。このとき，大円の径は小円の径の 2 倍になる。
   r1 = 0.25,  r2 = 0.375,  r3 = 0.5
   x2 = 0.6123724356957945,  y3 = 0.5

using SymPy
@syms r1::positive, r2::positive, r3::positive, x2::positive, y3::positive;

r1 = 16 // 64
eq1 = x2^2 + (r2 - r1)^2 - (r1 + r2)^2
eq2 = x2^2 + (2r1 + r3 - r2)^2 - (r2 + r3)^2
eq3 = x2^2 + (y3 - r2)^2 - (1 - r2)^2
eq4 = 2r1 + 2r3 - y3 - 1;
res = solve([eq1, eq2, eq3, eq4], (r2, r3, x2, y3))
println("parameter = [$r1, $(res[1][1]), $(res[1][2]), $(res[1][3]), $(res[1][4])]")

parameter = [1//4, 3/8, 1/2, sqrt(6)/4, 1/2]
draw(parameter, false)

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)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(parameter, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1, r2, r3, x2, y3 = parameter
   println("r1 = $r1,  r2 = $r2,  r3 = $r3")
   println("x2 = $x2,  y3 = $y3")
   plot()
   circle(0, r1, r1, :green)
   circle(x2, r2, r2, :blue)
   circle(0,  2r1 + r3, r3, :magenta)
   if more
       point(0, 2r1+r3, "2r1+r3 ", :magenta, :right)
       point(0, r1, "r1 ", :green, :right)
       point(x2, r2, "(x2,r2)", :blue, :center)
       point(0, y3, "y3 ", :black, :right)
       vline!([0], color=:black, linewidth=0.25)
   end
   circle(0, y3, 1, :black)
   hline!([0], color=:black, linewidth=0.25)
end;

Julia v1.8.3 リリース

Current stable release: v1.8.3 (November 14, 2022)

算額（その31）

算額（その31）

大阪府茨木市 井於神社 弘化3年(1846)
http://www.wasan.jp/osaka/iyo.html

(3) 奈良県大和郡山市小泉町 耳成山口神社 嘉永7年(1854)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.

キーワード：円5個，直線上

図のように，大円，中円とそれらに接する 3 個の円（甲，乙，丙）がある。大円，中円の径を三寸五分，二寸 としたとき，甲，乙，丙の円の径を求めよ。

それぞれの円の中心座標と半径を (0, r1, r1), (x2, r2, r2), (x3, r3, r3), (x4, r4, r4), (x5, r5, r5) とする。

using SymPy

@syms r1, r2, r3::positive, r4::positive, r5::positive
@syms x2::positive, x3::positive, x4::positive, x5::positive;

r1 = 35//10;
r2 = 2;
eq1 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2;
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2;
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2;
eq4 = x4^2 + (r1 - r4)^2 - (r1 + r4)^2;
eq5 = (x3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2;
eq6 = (x2 - x5)^2 + (r2 - r5)^2 - (r2 + r5)^2;
eq7 = (x5 - x3)^2 + (r3 - r5)^2 - (r3 + r5)^2;

res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r3, r4, r5, x2, x3, x4, x5))

   3-element Vector{NTuple{7, Sym}}:
    (154/9 - 56*sqrt(7)/9, -505064*sqrt(7)*sqrt(296 - 110*sqrt(7))/177147 - 1326104*sqrt(296 - 110*sqrt(7))/177147 - 56*sqrt(7)/81 + 242/81 + 513812*sqrt(2072 - 770*sqrt(7))/177147 + 195692*sqrt(7)*sqrt(2072 - 770*sqrt(7))/177147, 7/9 - 7*sqrt(7)/36, 2*sqrt(7), -28/3 + 14*sqrt(7)/3, 2*sqrt(7)/9 + 4*sqrt(14)*sqrt(148 - 55*sqrt(7))/(9*(-1 + sqrt(7))) + 28/9, -7/3 + 7*sqrt(7)/3)
    (154/9 - 56*sqrt(7)/9, -195692*sqrt(7)*sqrt(2072 - 770*sqrt(7))/177147 - 513812*sqrt(2072 - 770*sqrt(7))/177147 - 56*sqrt(7)/81 + 242/81 + 1326104*sqrt(296 - 110*sqrt(7))/177147 + 505064*sqrt(7)*sqrt(296 - 110*sqrt(7))/177147, 7/9 - 7*sqrt(7)/36, 2*sqrt(7), -28/3 + 14*sqrt(7)/3, -4*sqrt(14)*sqrt(148 - 55*sqrt(7))/(9*(-1 + sqrt(7))) + 2*sqrt(7)/9 + 28/9, -7/3 + 7*sqrt(7)/3)
    (56*sqrt(7)/9 + 154/9, -195692*sqrt(7)*sqrt(770*sqrt(7) + 2072)/177147 - 505064*sqrt(7)*sqrt(110*sqrt(7) + 296)/177147 + 56*sqrt(7)/81 + 242/81 + 1326104*sqrt(110*sqrt(7) + 296)/177147 + 513812*sqrt(770*sqrt(7) + 2072)/177147, 7*sqrt(7)/36 + 7/9, 2*sqrt(7), 28/3 + 14*sqrt(7)/3, -28/9 + 2*sqrt(7)/9 + 4*sqrt(14)*sqrt(55*sqrt(7) + 148)/(9*(1 + sqrt(7))), 7/3 + 7*sqrt(7)/3)

題意を満たすのは res[2] である。それにしても，r4 の式はものすごい。

name = ["r3", "r4", "r5", "x2", "x3", "x4", "x5"]
j = 2
for i in 1:7
   println("$(name[i]) = $(res[j][i]) = $(res[j][i].evalf())")
end

   r3 = 154/9 - 56*sqrt(7)/9 = 0.648658508931436
   r4 = -195692*sqrt(7)*sqrt(2072 - 770*sqrt(7))/177147 - 513812*sqrt(2072 - 770*sqrt(7))/177147 - 56*sqrt(7)/81 + 242/81 + 1326104*sqrt(296 - 110*sqrt(7))/177147 + 505064*sqrt(7)*sqrt(296 - 110*sqrt(7))/177147 = 0.316985841490936
   r5 = 7/9 - 7*sqrt(7)/36 = 0.263326133959663
   x2 = 2*sqrt(7) = 5.29150262212918
   x3 = -28/3 + 14*sqrt(7)/3 = 3.01350611830142
   x4 = -4*sqrt(14)*sqrt(148 - 55*sqrt(7))/(9*(-1 + sqrt(7))) + 2*sqrt(7)/9 + 28/9 = 2.10660907167730
   x5 = -7/3 + 7*sqrt(7)/3 = 3.84008639248404

算額の解答では，
甲: 0.625，
乙: 0.36525，
丙: 0.19315 
となっているが，かなりの誤差がある。

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")
   r1, r2 = 35/10, 2
   r3, r4, r5, x2, x3, x4, x5 = (154/9 - 56*sqrt(7)/9, -195692*sqrt(7)*sqrt(2072 - 770*sqrt(7))/177147 - 513812*sqrt(2072 - 770*sqrt(7))/177147 - 56*sqrt(7)/81 + 242/81 + 1326104*sqrt(296 - 110*sqrt(7))/177147 + 505064*sqrt(7)*sqrt(296 - 110*sqrt(7))/177147, 7/9 - 7*sqrt(7)/36, 2*sqrt(7), -28/3 + 14*sqrt(7)/3, -4*sqrt(14)*sqrt(148 - 55*sqrt(7))/(9*(-1 + sqrt(7))) + 2*sqrt(7)/9 + 28/9, -7/3 + 7*sqrt(7)/3)
   println("r1 = $r1;  r2 = $r2;  r3 = $r3;  r4 = $r4;  r5 = $r5")
   plot()
   circle(0, r1, r1, :red)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :brown)
   circle(x4, r4, r4, :green)
   circle(x5, r5, r5, :red)
   hline!([0], color=:black, linewidth=0.25)
   if more
       point(0, r1, "大円 r1")
       point(x2, r2,"中円 r2")
       point(x3, r3,"甲")
       point(x4, r4,"乙")
       point(x5, r5,"丙")
   end
end;

なお，フォードの円はデカルトの円定理を簡約化したものであるが，それを順次用いることで円の径を求めることができる。しかし，この方法では円の中心座標についての情報は得られないので，図を描くためには別途中心座標を求める必要がある。

using SymPy

@syms r1, r2, r3::positive, r4::positive, r5::positive
r1 = 35//10
r2 = 2
eq1 = 1/sqrt(r1) + 1/sqrt(r2) - 1/sqrt(r3);
r3 = solve(eq1)[1]
    0.648658508931432

eq2 = 1/sqrt(r1) + 1/sqrt(r3) - 1/sqrt(r4);
r4 = solve(eq2)[1]
    0.316985841490935

eq3 = 1/sqrt(r2) + 1/sqrt(r3) - 1/sqrt(r5);
r5 = solve(eq3)[1]
    0.263326133959661

デカルトの円定理では「曲率=円の径の逆数」を用いる。

using SymPy
@syms r1, r2, r3::positive, r4::positive, r5::positive
@syms k1, k2, k3, k4, k5
r1 = 35//10
r2 = 2
k1, k2, k3, k4, k5 = 1 ./ [r1, r2, r3, r4, r5]

まず r1, r2 から r3 を求める。

k3 =  k1 + k2 + 2√(k1*k2)
1/k3.evalf()
    0.648658508931436

ついで，r1, r3 から r4 を求める。

k4 = k1 + k3 + 2√(k1*k3)
1/k4.evalf()
    0.316985841490937

更に r2,r3 から r5 を求める。

k5 = k2 + k3 + 2√(k2*k3)
1/k5.evalf()
    0.263326133959663

 

 

算額（その30）

算額（その30）

奈良県奈良市 円満寺 天保 15 年
http://www.wasan.jp/nara/enman.html

奈良県奈良市下山町 円満寺 天保15年(1844)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.

図のように，1 個の大円と，6 個の小円がある。それぞれの径を求めよ。

見ただけでわかる。小円の径は大円の径の 1/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 draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1, r2 = 1, 1/3
   println("r1 = $r1;  r2 = $r2")
   plot()
   circle(0, 0, r1)
   for degree = 30:60:330
       circle(-4r1/3*cosd(degree), 4r1/3*sind(degree), r2, :blue)
   end
   plot!([0, 2r1*cosd(210), 2r1*cosd(330), 0], [2r1, 2r1*sind(210), 2r1*sind(330), 2r1], color=:black, linewidth=0.25)
   plot!([0, 2r1*cosd(210), 2r1*cosd(330), 0], [-2r1, -2r1*sind(210), -2r1*sind(330), -2r1], color=:black, linewidth=0.25)
end;