算額（その240）

算額（その240）

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

外円の中に天円 4 個，地円 6 個，人円 4 個が入っている。天円の径が与えられたとき，人円の径を求めよ。

外円，天円，地円，人円の半径と中心座標を以下のようにおく。

外円: R, (0, 0)
天円: r1, (x1, r1)
地円: r2, (0, 0) および (x2, r2)
人円: r3, (x3, y3)

天円の半径が与えれれるので，r1 = 1 とおく。
また，x2 = 2x1 である。
注: 天円は，y 軸とは接していない。

include("julia-source.txt");

以下の 3 本の方程式は，x1, r2, R についてのもので，solve() で解を求めることができる。

using SymPy

@syms r1::positive, x1::positive, r2::positive, x2::positive, r3::positive, x3::positive, y3::positive, R::positive;

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

solve([eq1, eq2, eq3], (x1, r2, R))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (sqrt(28 - 12*sqrt(5)), 7 - 3*sqrt(5), -2 + 2*sqrt(5))

しかし，x1, r2, R が決まったあと，r3, x3, y3 を求めるための残りの 3 本の方程式は solve() で解くことができない。

eq4 = x3^2 + y3^2 - (R - r3)^2
eq5 = (x3 - x1)^2 + (r1 - y3)^2 - (r1 + r3)^2
eq6 = (x3 - x2)^2 + (y3 - r2)^2 - (r3 + r2)^2;

それで，結局のところ 6 本の連立方程式を nlsolve() で解くことにする。

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

   x1^2 - (R - 1)^2 + 1,  # eq1
   r2^2 + 4*x1^2 - (R - r2)^2,  # eq2
   x1^2 + (1 - r2)^2 - (r2 + 1)^2,  # eq3
   x3^2 + y3^2 - (R - r3)^2,  # eq4
   (1 - y3)^2 - (r3 + 1)^2 + (-x1 + x3)^2,  # eq5
   (-r2 + y3)^2 - (r2 + r3)^2 + (-2*x1 + x3)^2,  # eq6

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)
   r1 = 1
   (x1, r2, r3, x3, y3, R) = u
   # println(u)
   return [
x1^2 - (R - 1)^2 + 1,  # eq1
r2^2 + 4*x1^2 - (R - r2)^2,  # eq2
x1^2 + (1 - r2)^2 - (r2 + 1)^2,  # eq3
x3^2 + y3^2 - (R - r3)^2,  # eq4
(1 - y3)^2 - (r3 + 1)^2 + (-x1 + x3)^2,  # eq5
(-r2 + y3)^2 - (r2 + r3)^2 + (-2*x1 + x3)^2,  # eq6
      ]
end;

iniv = [big"1.1", 0.3, 0.14, 2.25, 0.75, 2.5]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[1.080363026950906962773783583238333577245845592041721754016512840267774995583403, 0.2917960675006310117247564105634979118556556424332539734005336141870321167000515, 0.1519553880558686263724085395974730657224106667858737670157781675475329832318468, 2.201116552320873166994524966387851524880342069547970579400755429309693477168625, 0.7337055174402905594576728270074182849955867099200287269565993335856897130533752, 2.472135954999580224200317035334459839068561852275454446131625344134375443426156], true)

   x1 = 1.08036;  r2 = 0.29180;  r3 = 0.15196;  x3 = 2.20112;  y3 = 0.73371;  R = 2.47214

人円の半径は「天円の径 × 0.15196」である。

using Printf
function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x1, r2, r3, x3, y3, R) = res[1]
   @printf("x1 = %.5f; r2 = %.5f; r3 = %.5f; x3 = %.5f; y3 = %.5f; R = %.5f\n", x1, r2, r3, x3, y3, R)
   r1 = 1
   x2 = 2x1
   plot()
   circle(0, 0, R)
   circle4(x1, r1, r1, :brown)
   circle(0, r2, r2, :green)
   circle(0, -r2, r2, :green)
   circle4(x2, r2, r2, :green)
   circle4(x3, y3, r3, :blue)
   if more == true
       point(x1, r1, "(x1,r1)", :brown)
       point(0, r2, " r2")
       point(x2, r2, "(x2,r2)", :green, :right)
       point(x3, y3, "(x3,y3)  ", :blue, :right)
       point(0, R, " R", :red)
       point(1.2x1, 1.2r1, "天円", :brown, mark=false)
       point(0.3r2, 1.7r2, "地円", :green, mark=false)
       point(x2+0.3r2, 1.7r2, "地円", :green, mark=false)
       point(x3+0.3r3, y3+1.7r3, "人円", :blue, mark=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その239）

算額（その239）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（225-224)
長野県東筑摩郡筑北村青柳 碩水寺豊川稲荷社 慶応4年(1868)

鉤股の中に全円 1 個，大円 3 個，小円 2 個がある。大円はそれぞれ部分的に重なっており，その重なった部分に小円がある。小円は大円に内接する。大円は鉤股弦に内接する。

股が 50 寸，大円の径が 15 寸，全円の径が 25 寸のとき，小円の径はいかほどか。

鉤を y，股を x，弦を z とする。全円，大円，小円の半径をそれぞれ r0, r1, r2 とする。

include("julia-source.txt");

using SymPy

@syms r2::positive, y::positive, z::positive;

x = 50
r0 = 25 // 2
r1 = 15 // 2

eq1 = x^2 + y^2 - z^2
eq2 = x + y - sqrt(x^2 + y^2) - 2r0;

鉤股弦の股と弦を求め，鉤股弦に内接する円の直径が「鉤+股-弦」になることから弦を求める。

(y, z) = solve([eq1, eq2], (y, z))[1]
(y, z) |> println

   (75/2, 125/2)

一番左の大円の中心座標 (x - (5r1 - 4r2), r1) で鉤股弦を3分割したときの鉤股弦の面積が，x*y/2 に等しいことから r2 を求める。

eq3 = x * r1 + z * r1 + (5*r1 - 4r2)y - x*y
r2 = solve(eq3, r2)[1]
r2 |> println

   5/2

小円の直径は 5 寸である。

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   x = 50
   r0 = 25 // 2
   r1 = 15 // 2
   (y, z) = (75/2, 125/2)
   r2 = 5 // 2
   plot([0, x, x, 0], [0, 0, y, 0], color=:black, lw=0.5)
   circle(x - r0, r0, r0)
   circle(x - r1, r1, r1, :green)
   circle(x - (3r1 - 2r2), r1, r1, :green)
   circle(x - (5r1 - 4r2), r1, r1, :green)
   circle(x - (2r1 - r2), r1, r2, :blue)
   circle(x - (4r1 - 3r2), r1, r2, :blue)
   if more == true
       point(x, y, " (x,y)", :black)
       point(x - r0, r0, "(x-r0,r0)", :red, :center)
       point(x - (5r1 - 4r2), 0.9r1, "(x-(5r1-4r2),r1)", :green, :center, mark=false)
       point(x - (4r1 - 3r2), 1.1r1, "(x-(4r1-3r2),r1)", :blue, :center, :bottom, mark=false)
       point(x - (5r1 - 4r2), r1, "")
       point(x - (4r1 - 3r2), r1, "", :blue)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;