算額（その146）

算額（その146）

岐阜県大垣市西外側町 大垣八幡宮 天保年間
http://www.wasan.jp/gifu/ogakihatiman.html

第1問: 等脚台形に甲乙の 2 円を入れる。その隙間に順次，丙円，丁円，戊円...を入れていく。甲円の径と最後に入れた円（黒円と呼ぼう）の径が分かっているときに，乙円から黒円まで何個の円が入ったか。

この問題は見た目は違うが，90度時計回りに回転すると本質的には算額（その36）と同じものである。

算額の「術」では，「円の個数は ⌊√(甲径 / 黒径)⌋」とされている。甲径 / 黒径 が平方数の場合には，「等脚台形」は「長方形」になる。

隙間に入れる円の半径と中心の x 座標に関して漸化式を作り，順次円を埋めていく。

漸化式の作成

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

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

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)

次の円の半径 と x 座標を求める関数

function nextcircle(r2, x2, R) # 
  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;

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(R, r, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   println("個数 = $(floor(Int, sqrt(R/r)))")
   circle(0, R, R, :blue)
   x = 2*sqrt(R)*sqrt(r) # x1
   i = 0
   while true
       i += 1
       println("i = $i, x = $x, r = $r")
       circle(x, r, r, :red)
       (r, x) = nextcircle(r, x, R)
       round(r, digits=10) > R && break
   end
   if more == true
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
   hline!([0], color=:black, lw=0.5)
end;

draw(225, 1, true)
savefig("fig1.png");

   個数 = 15
   i = 1, x = 30.0, r = 1
   i = 2, x = 32.142857142857146, r = 1.1479591836734695
   i = 3, x = 34.61538461538462, r = 1.3313609467455623
   i = 4, x = 37.500000000000014, r = 1.562500000000001
   i = 5, x = 40.90909090909093, r = 1.8595041322314063
   i = 6, x = 45.00000000000002, r = 2.2500000000000018
   i = 7, x = 50.00000000000002, r = 2.7777777777777803
   i = 8, x = 56.25000000000002, r = 3.515625000000003
   i = 9, x = 64.2857142857143, r = 4.5918367346938815
   i = 10, x = 75.00000000000003, r = 6.250000000000005
   i = 11, x = 90.00000000000003, r = 9.000000000000007
   i = 12, x = 112.50000000000003, r = 14.062500000000012
   i = 13, x = 150.00000000000003, r = 25.00000000000002
   i = 14, x = 225.00000000000006, r = 56.25000000000004
   i = 15, x = 450.00000000000034, r = 225.00000000000034

draw(225, 2, false)

   個数 = 10
   i = 1, x = 42.42640687119285, r = 2
   i = 2, x = 46.8427872018762, r = 2.438051903155841
   i = 3, x = 52.285444912527154, r = 3.0375197218899017
   i = 4, x = 59.159137903470494, r = 3.8886706638687114
   i = 5, x = 68.11368676693468, r = 5.154971472204546
   i = 6, x = 80.26252312010512, r = 7.157858464006006
   i = 7, x = 97.68589245763658, r = 10.602815094716597
   i = 8, x = 124.77119327573541, r = 17.297611857167695
   i = 9, x = 172.63857263936492, r = 33.11564084773032
   i = 10, x = 280.0943102542604, r = 87.16980292978876

draw(10000, 5, false)

   個数 = 44
   i = 1, x = 447.21359549995793, r = 5
   i = 2, x = 457.44231665828704, r = 5.231336826742515
   i = 3, x = 468.1498952905389, r = 5.479108111513565
   i = 4, x = 479.37076393767734, r = 5.744908232954808
   i = 5, x = 491.142737399151, r = 6.030529712498285
   i = 6, x = 503.5074384809699, r = 6.33799351514169
   i = 7, x = 516.5107897148246, r = 6.669584897295796
   i = 8, x = 530.203583290577, r = 7.027895993354192
   i = 9, x = 544.6421441115267, r = 7.415876628560024
   i = 10, x = 559.8891042209045, r = 7.836895225632172
   i = 11, x = 576.0143110525872, r = 8.294812163434667
   i = 12, x = 593.0958972857454, r = 8.794068584429588
   i = 13, x = 611.2215468749545, r = 9.339794484105308
   i = 14, x = 630.4900005461035, r = 9.93794101971564
   i = 15, x = 651.0128553214608, r = 10.595443444845037
   i = 16, x = 672.9167273238959, r = 11.320423047807564
   i = 17, x = 696.345866399034, r = 12.122439141275532
   i = 18, x = 721.4653366451872, r = 13.01280579951383
   i = 19, x = 748.4649110978014, r = 14.004993078615996
   i = 20, x = 777.5638749236822, r = 15.115139489658288
   i = 21, x = 809.016994374947, r = 16.362712429686827
   i = 22, x = 843.1219955098877, r = 17.771367482814387
   i = 23, x = 880.2290178099703, r = 19.37007809486762
   i = 24, x = 920.7526791297857, r = 21.194637403116957
   i = 25, x = 965.1876341291538, r = 23.289679226895828
   i = 26, x = 1014.1288661817565, r = 25.71143393057737
   i = 27, x = 1068.2984826283346, r = 28.531541199650054
   i = 28, x = 1128.5815822186587, r = 31.842409693079276
   i = 29, x = 1196.0749925986145, r = 35.76488469799439
   i = 30, x = 1272.154608282929, r = 40.459433684387335
   i = 31, x = 1358.5701736362862, r = 46.14282291735322
   i = 32, x = 1457.5815120307168, r = 53.11359660534376
   i = 33, x = 1572.1590071720439, r = 61.79209859580466
   i = 34, x = 1706.2867080130786, r = 72.78535824855273
   i = 35, x = 1865.435060426439, r = 86.9961991167048
   i = 36, x = 2057.325407741825, r = 105.81469583350167
   i = 37, x = 2293.220441761242, r = 131.4714998627906
   i = 38, x = 2590.217418383382, r = 167.73065686241685
   i = 39, x = 2975.5884730214275, r = 221.3531690194499
   i = 40, x = 3495.672632567669, r = 305.4931788520645
   i = 41, x = 4236.0679774997825, r = 448.6067977499774
   i = 42, x = 5374.379909090633, r = 722.0989851809261
   i = 43, x = 7349.267758474197, r = 1350.2934146437085
   i = 44, x = 11618.723119204686, r = 3374.8681730185353

draw(225, 5, false)

   個数 = 6
   i = 1, x = 67.0820393249937, r = 5
   i = 2, x = 78.8339038551072, r = 6.905315996706985
   i = 3, x = 95.57784803962724, r = 10.150138928762306
   i = 4, x = 121.35254915624216, r = 16.362712429686848
   i = 5, x = 166.16178515947303, r = 30.677487608203194
   i = 6, x = 263.434588481236, r = 77.10864712030906

 

算額（その145）

算額（その145）

岐阜県大垣市西外側町 大垣八幡宮 天保年間
http://www.wasan.jp/gifu/ogakihatiman.html

第6問: 正三角形の頂点から対辺へ2本の斜線を描き，甲乙2種類の円を入れる。甲径を知って乙径を求めよ。

正三角形の一辺の長さを 2 とし，甲円，乙円の半径と中心の x 座標を r1, x1, r2, x2 とする。
B の x 座標を x とする。BO = OC。
各円の中心点から斜線及び斜辺への二乗距離が半径に等しいことを表す方程式を立て，解く。

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 x::positive, x1::positive, r1::positive, x2::positive, r2::positive;
eq1 = distance(0, sqrt(Sym(3)), 1, 0, x1, r1) - r1^2   # 右の甲円の中心から AD への二乗距離
eq2 = distance(0, sqrt(Sym(3)), x, 0, x2, r2) - r2^2   # 右の乙円の中心から AC への二乗距離
eq3 = distance(0, sqrt(Sym(3)), x, 0, 0, r2) - r2^2    # 真ん中の乙円から AC への二乗距離
eq4 = distance(0, sqrt(Sym(3)), -x, 0, x1, r1) - r1^2  # 右の甲円の中心から AB への二乗距離
eq5 = distance(0, sqrt(Sym(3)), 1, 0, x2, r2) - r2^2;  # 右の乙円の中心から AD への二乗距離

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

   1-element Vector{NTuple{5, Sym}}:
    (3^(1/3)*(-1 + 3^(1/3))/2, -3^(2/3)/4 - 1/2 + 3^(1/3)/4 + sqrt(3^(2/3) + 3*3^(1/3) + 6)/4, sqrt(3)*(-sqrt(3^(2/3)*(1 - 3^(1/3))^2 + 12) - 3^(1/3) + 3^(2/3) + 6)/12, 3^(1/3)*(1 - 3^(1/3))*sqrt(3^(2/3)*(1 - 3^(1/3))^2 + 12)/4 + 3^(2/3)*(1 - 3^(1/3))^2/4 + 1, 3^(5/6)*(1 - 3^(1/3))*(-sqrt(3^(2/3)*(1 - 3^(1/3))^2 + 12) + 3^(1/3)*(-1 + 3^(1/3)))/12)

算額は「甲径を知って乙径を求めよ」なので，「乙径 / 甲径」を求める。

f = res[1][5] / res[1][3]
f |>  println

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

容易には簡単な式にならない。数値としては以下のようになる。

f.evalf() |> println

   0.590541436813876

算額の「術文」では 1 / (1/∛3 + 1) である。もう少し簡単にすれば ∛3/(1+∛3) である。これが前の長々しい複雑な式と同じ値になる。

1 / (1/3^(1/3) + 1) |> println

   0.590541436813876

∛3/(1+∛3) |> println

   0.590541436813876

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, 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")
   (x, x1, r1, x2, r2) = res[1]
   s3 = sqrt(3)
   plot([1, 0, -1, 1], [0, s3, 0, 0], color=:black, lw=0.5)
   circle(0, r2, r2, :red)
   circle(x2, r2, r2, :red)
   circle(-x2, r2, r2, :red)
   circle(x1, r1, r1, :blue)
   circle(-x1, r1, r1, :blue)
   segment(0, s3, x, 0, :magenta, linewidth=0.25)
   segment(0, s3, -x, 0, :magenta, linewidth=0.25)
   if more == true
       point(0, s3, " A")
       point(-x, 0, " B")
       point(x, 0, " C")
       point(1, 0, " D")
       point(0, 0, " O")
       point(x1, r1, "甲:(x1,r1)", :green, :left, :bottom)
       point(x2, r2, "乙:(x2,r2)", :green, :center, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5, ylims=(-0.1, 1.8))
   else
       plot!(showaxis=false)
   end
end;

 

算額（その144）

算額（その144）

岐阜県大垣市西外側町 大垣八幡宮 天保年間
http://www.wasan.jp/gifu/ogakihatiman.html

第5問: 長方形の中に2本の斜線を入れ，甲乙丙丁の4種類の円を入れる。長方形の短辺の長さが与えられるとき，長辺の長さを求めよ。

長方形の短辺の長さを y = 1 とする。甲乙丙丁の円の半径を r1, r2, r3, r4 とする。また，右下の乙丙の中心の x 座標をx2, x3，右上の丁円の中心座標を x4 とする。
以下の方程式を立て，解く。

using SymPy

@syms x::positive, y::positive, x2::positive, x3::positive, x4::positive, r1::positive, r2::positive, r3::positive, r4::positive;
y = 1
eq1 = r1/x - r2/(x-x2)  # 右下の甲，乙，丙の半径について tanθ = r1/x との関係
eq2 = r1/x - r3/(x - x3)  # 右下の甲，乙，丙の半径について 同上
eq3 = r1/x - (r3 - r4)/x4  # 上の丙，丁の半径について 同上
eq4 = x2^2 + (r1 - r2)^2- (r1 + r2)^2  # 甲乙が外接
eq5 = (x3 - x2)^2 + (r2 - r3)^2 - (r2 + r3)^2  # 乙丙が外接
eq6 = x4^2 + (r3 - r4)^2 - (r3 + r4)^2  # 丙丁が外接
eq7 = r4 - r1^3/x^2  # 甲丁の半径について sin(π-2θ) = (r1-r4)/(r1+r4)=cos(2θ) を r4 について解く
eq8 = 5　//　2 * r1　/　(1　-　r1^2　/　x^2) - y ## x*tanθ + x/4 * tanθ = y
eq9 = 20　//　7　*　r1 - y;  # 精査すれば r1 は単純に求められる

eq9 の代わりに eq8 を使うと solve() では解けないので，nlsolve() で数値解を求める。

# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (x, x2, x3, x4, r1, r2, r3, r4))

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

   r1/x - r2/(x - x2),
   r1/x - r3/(x - x3),
   r1/x - (r3 - r4)/x4,
   x2^2 + (r1 - r2)^2 - (r1 + r2)^2,
   (r2 - r3)^2 - (r2 + r3)^2 + (-x2 + x3)^2,
   x4^2 + (r3 - r4)^2 - (r3 + r4)^2,
   -r1^3/x^2 + r4,
   5*r1/(2*(-r1^2/x^2 + 1)) - 1,

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
       println(params...)
       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)
   (x, x2, x3, x4, r1, r2, r3, r4) = u
   return [
       r1/x - r2/(x - x2),
       r1/x - r3/(x - x3),
       r1/x - (r3 - r4)/x4,
       x2^2 + (r1 - r2)^2 - (r1 + r2)^2,
       (r2 - r3)^2 - (r2 + r3)^2 + (-x2 + x3)^2,
       x4^2 + (r3 - r4)^2 - (r3 + r4)^2,
       -r1^3/x^2 + r4,
       5*r1/(2*(-r1^2/x^2 + 1)) - 1,
   ]
end;

iniv = [1, 0.5, 0.7, 0.1, 0.4, 0.2, 0.1, 0.4]
res = nls(H, ini=iniv)
println(res)

   ([0.9899494936611649, 0.494974746830583, 0.742462120245875, 0.12374368670764607, 0.3499999999999999, 0.17499999999999988, 0.08750000000000022, 0.04375000000000009], true)

eq8 の代わりに eq9 を使うと，solve() で解析解が求まる。

res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq9], (x, x2, x3, x4, r1, r2, r3, r4))

   1-element Vector{NTuple{8, Sym}}:
    (7*sqrt(2)/10, 7*sqrt(2)/20, 21*sqrt(2)/40, 7*sqrt(2)/80, 7/20, 7/40, 7/80, 7/160)

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, 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")
   y = 1
   (x, x2, x3, x4, r1, r2, r3, r4) = res[1]
   plot([x, x, -x, -x, x], [0, y, y, 0, 0], color=:black, lw=0.5)
   circle(0, r1, r1, :red)
   circle(x2, r2, r2, :blue)
   circle(-x2, r2, r2, :blue)
   circle(x3, r3, r3, :orange)
   circle(-x3, r3, r3, :orange)
   circle(0, 2r1 + r4, r4, :green)
   circle(0, y - r3, r3, :orange)
   circle(x4, y - r4, r4, :green)
   circle(-x4, y - r4, r4, :green)
   if more == true
       println("長方形の長辺 = $(2x)")
       println("甲: 半径 = $r1")
       println("乙: 半径 = $r2")
       println("丙: 半径 = $r3")
       println("丁: 半径 = $r4")
       point(0, r1, "甲", :red)
       point(x2, r2, "乙", :blue)
       point(x3, r3, "丙", :orange)
       point(0, 2r1 + r4, "丁", :green)
       point(0, y - r3, "丙", :orange)
       point(x4, y - r4, "丁", :green)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
   segment(x, 0, -x/4, y, :black)
   segment(-x, 0, x/4, y, :black)
end;

   長方形の長辺 = 7*sqrt(2)/5
   甲: 半径 = 7/20
   乙: 半径 = 7/40
   丙: 半径 = 7/80
   丁: 半径 = 7/160