算額（その255）

算額（その255）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県外(285)
東京都あきる野市 二宮神社 寛政6年(1794)

勾股弦（直角三角形）に 5 個の等円（直径が同じ円）が弦（直角に対する変）及び鉤（直角を挟む短い方の辺），股（同じく長い方の辺）に接している。内接する全円の径が 44 寸のとき，等円の径を求めよ。

全円の半径を r1 = 22 寸 とする。
等円の半径を r2 とする。

以下の方程式をとき，まず勾，股，弦の長さを求める。

include("julia-source.txt");

using SymPy

# 全円半径 r1 = 22 寸、
# 等円半径 r2
@syms 鉤, 股, 弦, r1::positive, r2::positive,  x::positive, y::positive
弦 = 110
eq1 = 鉤 + 股 - 弦 - 44  # 鉤股弦と内接円の径の関係
eq2 = 鉤^2 + 股^2 - 弦^2  # ピタゴラスの定理
(鉤, 股) = solve([eq1, eq2], (鉤, 股))[1]

   (66, 88)

eq3 = distance(0, 鉤, 股, 0, r2, y) - r2^2
eq4 = distance(0, 鉤, 股, 0, x, r2) - r2^2
eq5= (x - r2)^2 + (y - r2)^2 -(8r2)^2;

res = solve([eq3, eq4, eq5], (r2, x, y))

   4-element Vector{Tuple{Sym, Sym, Sym}}:
    (110/13, 814/13, 638/13)
    (660/53, 4884/53, 3828/53)
    (-1540/191 + 880*sqrt(15)/191, 21428/191 - 2640*sqrt(15)/191, 440*sqrt(15)/191 + 11836/191)
    (-2310/491 + 1980*sqrt(15)/491, 660*sqrt(15)/491 + 42438/491, 37026/491 - 3960*sqrt(15)/491)

names = ["等円半径", "内側の股の x 座標", "内側の鉤の y 座標"]
for i = 1:4
   println(i)
   for j = 1:3
       println("$(names[j]) = $(res[i][j]) = $(res[i][j].evalf())")
   end
end

   1
   等円半径 = 110/13 = 8.46153846153846
   内側の股の x 座標 = 814/13 = 62.6153846153846
   内側の鉤の y 座標 = 638/13 = 49.0769230769231
   2
   等円半径 = 660/53 = 12.4528301886792
   内側の股の x 座標 = 4884/53 = 92.1509433962264
   内側の鉤の y 座標 = 3828/53 = 72.2264150943396
   3
   等円半径 = -1540/191 + 880*sqrt(15)/191 = 9.78128452702894
   内側の股の x 座標 = 21428/191 - 2640*sqrt(15)/191 = 58.6561464189132
   内側の鉤の y 座標 = 440*sqrt(15)/191 + 11836/191 = 70.8906422635145
   4
   等円半径 = -2310/491 + 1980*sqrt(15)/491 = 10.9134562637285
   内側の股の x 座標 = 660*sqrt(15)/491 + 42438/491 = 91.6378187545762
   内側の鉤の y 座標 = 37026/491 - 3960*sqrt(15)/491 = 44.1730874725430

一番目の解が適切である。

   (r2, x, y) = (110/13, 814/13, 638/13)

   (8.461538461538462, 62.61538461538461, 49.07692307692308)

   等円の半径 = 8.46154;  等円の直径 = 16.92308;  x = 62.61538;  y = 49.07692

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鉤, 股, r2, x, y) = (66, 88, 110/13, 814/13, 638/13)
   @printf("等円の半径 = %.5f;  等円の直径 = %.5f;  x = %.5f;  y = %.5f\n", r2, 2r2, x, y)
   plot([0, 88, 0, 0], [0, 0, 66, 0], color=:black, lw=0.5)
   circle(22, 22, 22, :red)
   for i = 0:4
       circle(r2*(1 + 2i*4/5) , y - 2i*r2*3/5, r2, :blue)
   end
   if more == true
       plot!([r2, x, r2, r2], [r2, r2, y, r2], color=:black, lw=0.5)
       circle(r2, r2, r2, :blue)
       point(股, 0, "88", :black, :left, :bottom)
       point(0, 鉤, "  66", :black)
       point(22, 22, " 全円\n(r1,r1)", :red) 
       point(r2, y, " 等円 (r2,y)\n", :blue, :center, :bottom)
       point(x, r2, "\n (x,r2)", :blue, :center, :top)
       point(r2, r2, " (r2,r2)", :blue)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その249）

算額（その249）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（230)
長野県諏訪市下諏訪 諏訪大社下社 慶応4年(1868)

長方形内に甲，乙，丙，丁の円が 1 個ずつ入っている。丙円の径が 5寸0分6厘2毛5糸，丁円の径が 1寸4分4厘のとき，甲円の径はいかほどか。

注：以下の説明図は算額の図に近くなるように，r3 = 2, r4 = 1 としたものである。

長方形の長辺，短辺の長さをそれぞれ x，y とする。左下隅角を原点とする。
甲円の半径，中心座標を r1, (x − r1, r1) とする。y = 2r1 である。
乙円の半径，中心座標を r2, (r2, r2) とする。
丙円の半径，中心座標を r3, (x3, y - r3) とする。r3 = 50625 / 2
丁円の半径，中心座標を r4, (x4, r4) とする。r4 = 14400 / 2

以下の連立方程式を nlsolve() で解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive, x3::positive,
   r4::positive, x4::positive, x::positive;

この段階では r3, r4 も未知数とする

y = 2r1
eq1 = (x3 - r2)^2 + (y - r3 - r2)^2 - (r2 + r3)^2      # 乙円と丙円が外接
eq2 = (x - r1 - r2)^2 + (r1 - r2)^2 - (r1 + r2)^2      # 甲円と乙円が外接
eq3 = (x4 - r2)^2 + (r2 - r4)^2 - (r2 + r4)^2          # 乙円と丁円が外接
eq4 = (x - r1 - x3)^2 + (y - r3 - r1)^2 - (r1 + r3)^2  # 甲円と丙円が外接
eq5 = (x - r1 - x4)^2 + (r1 - r4)^2 - (r1 + r4)^2;     # 甲円と丁円が外接

# res = solve([eq1, eq2, eq3, eq4, eq5])

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
println(eq5, ",  # eq5")

   (-r2 + x3)^2 - (r2 + r3)^2 + (2*r1 - r2 - r3)^2,  # eq1
   (r1 - r2)^2 - (r1 + r2)^2 + (-r1 - r2 + x)^2,  # eq2
   (-r2 + x4)^2 + (r2 - r4)^2 - (r2 + r4)^2,  # eq3
   (r1 - r3)^2 - (r1 + r3)^2 + (-r1 + x - x3)^2,  # eq4
   (r1 - r4)^2 - (r1 + r4)^2 + (-r1 + x - x4)^2,  # eq5

算額の通り，r3 = 50625 / 2, r4 = 14400 / 2 とすると，算額の図とは似ても似つかない図が得られる。しかし，求まる r1, r2 は算額の通りである。

算額の多くは，実際に図を描くためのすべてのパラメータを計算しないので，実際とは異なる図を掲示してしまう。
算額では r1 = r4*(sqrt(2sqrt(r3/r4) + 1/4) + 1/2)^2
算額の図は r3/r4≒2 の場合のもの（前述の説明図）であるが，問では 5.0625/1.4400≒3.5 なので違って当然なわけである。

r3 = 50625 / 2
r4 = 14400 / 2
f(r3, r4) = r4*(sqrt(2sqrt(r3/r4) + 1/4) + 1/2)^2
f(50625 / 2, 14400 / 2), f(50625 / 2, 50625 / 4)

   (45000.0, 64331.362230694234)

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)
   (x, r1, r2, x3, x4) = u
   return [
       (-r2 + x3)^2 - (r2 + r3)^2 + (2*r1 - r2 - r3)^2,  # eq1
       (r1 - r2)^2 - (r1 + r2)^2 + (-r1 - r2 + x)^2,  # eq2
       (-r2 + x4)^2 + (r2 - r4)^2 - (r2 + r4)^2,  # eq3
       (r1 - r3)^2 - (r1 + r3)^2 + (-r1 + x - x3)^2,  # eq4
       (r1 - r4)^2 - (r1 + r4)^2 + (-r1 + x - x4)^2,  # eq5
   ]
end;

(r3, r4) = (50625 // 2, 14400 // 2)
(r3, r4) = (2, 1)
iniv = [big"200000.0", 89000, 50000, 50000, 60000]
res = nls(H, ini=iniv);
println(res);
   (BigFloat[125000.0000000000000000000000000000000078311250578101097708795336203908490939767, 45000.00000000000000000000000000000000240863585814089291883850249994466664153376, 20000.00000000000000000000000000000000145859703738035599114200719600180134501038, 12500.00000000000000000000000000000000319845103367027695329127846269509477692267, 44000.00000000000000000000000000000000423640766899408968775137416801857989402613], true)

names = ["x", "r1", "r2", "x3", "x4"]
for i in 1:length(names)
   println("$(names[i]) = $(Float64(res[1][i]))")
end

   x = 125000.0
   r1 = 45000.0
   r2 = 20000.0
   x3 = 12500.0
   x4 = 44000.0

甲円の半径は 45000 なので，元の単位での直径は 9 寸。

using Plots
using Printf

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   # res, r3, r4 は　グローバル変数
   (x, r1, r2, x3, x4) = res[1]
   @printf("r3 = %.5f;  r4 = %.5f\n", r3, r4)
   @printf("x = %.5f;  r1 = %.5f;  r2 = %.5f;  x3 = %.5f;  x4 = %.5f\n", x, r1, r2, x3, x4)
   y = 2r1
   plot([0, x, x, 0, 0], [0, 0, y, y, 0], color=:black, lw=0.5)
   circle(x - r1, r1, r1)
   circle(r2, r2, r2, :blue)
   circle(x3, y - r3, r3, :green)
   circle(x4, r4, r4, :magenta)
   if more == true
       point(x − r1, r1, " 甲円:r1\n (x−r1,r1)\n", :red, :center, :bottom)
       point(r2, r2, "\n乙円:r2\n (r2,r2)\n", :blue, :center, :bottom)
       point(x3, y - r3, "丙円:r3\n(x3,y-r3)\n", :green, :center, :bottom)
       point(x4, r4, "   丁円:r4\n      (x4,r4)", :magenta, :left, :bottom)
       point(x, 0, " x", :black, :left, :bottom)
       point(0, y, " y", :black, :left, :top)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;