算額（その312）

算額（その312）

「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県菰野町 広幡神社 嘉永5年(1852)

問題
半円の中に円弧があり，半円と円弧に挟まれた領域にいくつかの円が挟まれている。半円の直径，丙円，丁円の直径が与えられたとき，戊円の直径はいくつか。



円弧の半径と中心座標を r0, (0, y0)
半円の半径と中心座標を r,  (0, 0)
丙円の半径と中心座標を r1, (x1, y1)
丁円の半径と中心座標を r2, (x2, y2)
戊円の半径と中心座標を r3, (x3, y3)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, y0::positive, r::positive,
   r1::positive, x1::positive, y1::negative,
   r2::positive, x2::positive, y2::negative,
   r3::positive, x3::positive, y3::negative;

(r, r1, r2) = (50, 10, 5)
eq1 = x1 ^2 + (y0 - y1)^2 - (r0 + r1)^2
eq2 = x2 ^2 + (y0 - y2)^2 - (r0 + r2)^2
eq3 = x3 ^2 + (y0 - y3)^2 - (r0 + r3)^2
eq4 = x1^2 + y1^2 - (r - r1)^2
eq5 = x2^2 + y2^2 - (r - r2)^2
eq6 = x3^2 + y3^2 - (r - r3)^2
eq7 = (x2 - x1)^2 + (y2 - y1)^2 - (r2 + r1)^2
eq8 = (x3 - x2)^2 + (y3 - y2)^2 - (r3 + r2)^2
eq9 = r^2 + y0^2 - r0^2;

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

連立方程式が複雑すぎて solve() では解けないので，nlsolve() で数値解を求める。

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

   x1^2 - (r0 + 10)^2 + (y0 - y1)^2,  # eq1
   x2^2 - (r0 + 5)^2 + (y0 - y2)^2,  # eq2
   x3^2 - (r0 + r3)^2 + (y0 - y3)^2,  # eq3
   x1^2 + y1^2 - 1600,  # eq4
   x2^2 + y2^2 - 2025,  # eq5
   x3^2 + y3^2 - (50 - r3)^2,  # eq6
   (-x1 + x2)^2 + (-y1 + y2)^2 - 225,  # eq7
   -(r3 + 5)^2 + (-x2 + x3)^2 + (-y2 + y3)^2,  # eq8
   -r0^2 + y0^2 + 2500,  # eq9

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)
   (r0, y0, x1, y1, x2, y2, r3, x3, y3) = u
   return [
       x1^2 - (r0 + 10)^2 + (y0 - y1)^2,  # eq1
       x2^2 - (r0 + 5)^2 + (y0 - y2)^2,  # eq2
       x3^2 - (r0 + r3)^2 + (y0 - y3)^2,  # eq3
       x1^2 + y1^2 - 1600,  # eq4
       x2^2 + y2^2 - 2025,  # eq5
       x3^2 + y3^2 - (50 - r3)^2,  # eq6
       (-x1 + x2)^2 + (-y1 + y2)^2 - 225,  # eq7
       -(r3 + 5)^2 + (-x2 + x3)^2 + (-y2 + y3)^2,  # eq8
       -r0^2 + y0^2 + 2500,  # eq9
   ]
end;

iniv = [big"66.0", 45.0, 31.2, -24.5, 42, -13, 3, 47, -7]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[76.12612612612612612612612611757250855206487617275733389169965733209449470577765, 57.40372007954591399580886043401747772776331606928505070760935723414311831179716, 33.42516087186933536638372029599284981641171160639776150802648249557181481914727, -21.97176872010205673632683968491622029305014289094789997653090282550024009759404, 43.6384044716071878394454126067690155978138981377281139789361700897187230488509, -10.9858843600510283681634198397671466786061193950162591863412958011058978766673, 2.226463104325699745547073791346257958834702049836892015992527235142629533435189, 47.52241383500506650373257942562149702505122375947591050606909041806671507003448, -4.8919332392084731919302760073024864943304938091114520230459338837681434939186], true)

   円弧を構成する円の直径 = 152.252252, 中心座標 = (0.000000, 57.403720)
   半円の直径 = 100.000000, 中心座標 = (0.000000, 0.000000)
   円弧を構成する円の直径 = 152.252252, 中心座標 = (0.000000, 57.403720)
   丙円の直径 = 20.000000, 中心座標 = (33.425161, -21.971769)
   丁円の直径 = 10.000000, 中心座標 = (43.638404, -10.985884)
   戊円の直径 = 4.452926, 中心座標 = (47.522414, -4.891933)

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, r1, r2) = (50, 10, 5)
   (r0, y0, x1, y1, x2, y2, r3, x3, y3) = res[1]
   @printf("円弧を構成する円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r0, 0, y0)
   @printf("半円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r, 0, 0)
   @printf("円弧を構成する円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r0, 0, y0)
   @printf("丙円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r1, x1, y1)
   @printf("丁円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r2, x2, y2)
   @printf("戊円の直径 = %.6f, 中心座標 = (%.6f, %.6f)\n", 2r3, x3, y3)
   plot()
   circle(0, y0, r0, :black, beginangle=180, endangle=360)
   circle(0, 0, r, beginangle=180, endangle=360)
   circle(x1, y1, r1, :blue)
   circle(x2, y2, r2, :orange)
   circle(x3, y3, r3, :magenta)
   if more
       point(0, y0, " y0", :black)
       point(r, 0, "r ", :red, :right, :bottom)
       point(x1, y1, "  丙円:r1,(x1,y1)", :blue, :left, :bottom) 
       point(x2, y2, "   丁円:r2,(x2,y2)", :orange, :left, :bottom) 
       point(x3, y3, "  戊円:r3,(x3,y3)", :magenta, :left, :bottom) 
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;