算額（その379）

算額（その379）

二八 都川村西平 坂東九番観音堂（慈光寺境内観音堂） 文政13年(1830)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

埼玉県都幾川町西平 慈光寺 文政13年(1830)
やまぶき　和算と歴史随想
https://yamabukiwasan.sakura.ne.jp/page3.html

キーワード：円5個，円弧2個

二本の等弧の間に，天円 2 個，人円 1 個，地円 2 個が互いに接している。天円の径は 68 寸(注），地円の径が 17 寸のとき，人円の径はいかほどか。

注：「六寸八寸」と描かれているが「六十八寸」の誤記である。

その一部が「弧」になる円の半径と中心座標を r1, (r1, 0)
天円の半径と中心座標を r2, (r2, y2)
人円の半径と中心座標を r3, (0, y3)
地円の半径と中心座標を r4, (r4, y4)
とし，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, y2::positive,  r3::positive,
     y3::positive, r4::positive, y4::positive;

(r2, r4) = (68//2, 17//2)
eq1 = (r1 - r2)^2 + y2^2 - (r1 + r2)^2
eq2 = r1^2 + y3^2 - (r1 + r3)^2
eq3 = (r1 - r4)^2 + y4^2 - (r1 + r4)^2
eq4 = r2^2 + (y2 - y3)^2 - (r2 + r3)^2
eq5 = r4^2 + (y3 - y4)^2 - (r3 + r4)^2

res = solve([eq1, eq2, eq3, eq4, eq5], (r1, y2, r3, y3, y4))

   1-element Vector{NTuple{5, Sym}}:
    (272, 136*sqrt(2), 32, 96*sqrt(2), 68*sqrt(2))

   r1 = 272;  y2 = 192.333;  r3 = 32;  y3 = 135.765;  y4 = 96.1665
   人円の直径 = 2r3 = 64 寸

人円の半径は 32 寸である（直径は 64 寸）

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r4) = (68//2, 17//2)
   (r1, y2, r3, y3, y4) = res[1]
   @printf("r1 = %g;  y2 = %g;  r3 = %g;  y3 = %g;  y4 = %g\n", r1, y2, r3, y3, y4)
   @printf("人円の直径 = 2r3 = %g 寸\n", 2r3)
   plot()
   circle(r1, 0, r1)
   circle(-r1, 0, r1)
   circle(r2, y2, r2, :blue)
   circle(-r2, y2, r2, :blue)
   circle(0, y3, r3, :green)
   circle(r4, y4, r4, :magenta)
   circle(-r4, y4, r4, :magenta)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(r1, 0, " r1", :red)
       point(r2, y2, "  天円:r2,(r2,y2)", :blue, :left, :vcenter)
       point(0, y3, "  人円:r3,(0,y2)", :green, :left, :vcenter)
       point(r4, y4, "  地円:r4,(r4,y4)", :magenta, :left, :vcenter)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       plot!(xlims=(-100, 350), ylims=(-100, 275))
   else
       plot!(showaxis=false)
   end
end;

算額（その378）

算額（その378）

三六 川越市古谷本郷 古尾谷八幡神社 天保12年(1841)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

埼玉県川越市 古尾谷八幡神社 天保12年
http://tamtom.blog44.fc2.com/blog-entry-406.html

算額[古尾谷八幡神社]
http://blog.livedoor.jp/koike631/archives/50527866.html

等脚台形内に，乾円，坤円が入っている。等脚台形の上底の長さが 9 寸，乾円と坤円の直径がそれぞれ 7寸，9 寸のとき下底の長さを求めよ。

上底と下底の長さをそれぞれ 2w1，2w2，高さを h
乾円の半径と中心座標を r1, (x1, h - r1)
坤円の半径と中心座標を r2, (x2, r2)
とする。

r1 = 7/2, r2 = 4/2, w1 = 9/2 として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms h::positive, w1::positive, w2::positive,  r1::positive,
     r2::positive, x1::positive, x2::positive;

eq1 = distance(-w2, 0, -w1, h, x1, r1) - r1^2
eq2 = distance(w2, 0, -w1, h, x2, h - r2) - r2^2
eq3 = distance(w2, 0, -w1, h, x1, r1) - r1^2
eq4 = distance(w2, 0, w1, h, x2, h - r2) - r2^2;

# solve([eq1, eq2, eq3, eq4], (h, w2, x1, x2))

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)
   (h, w2, x1, x2) = u
   return [
       -r1^2 + (x1 - (-h*(h*w2 + h*x1 + r1*w1 - r1*w2) + x1*(h^2 + w1^2 - 2*w1*w2 + w2^2))/(h^2 + w1^2 - 2*w1*w2 + w2^2))^2 + (-h*(h*r1 - w1*w2 - w1*x1 + w2^2 + w2*x1)/(h^2 + w1^2 - 2*w1*w2 + w2^2) + r1)^2,  # eq1
       -r2^2 + (x2 - (-h^2*w1 + h*r2*w1 + h*r2*w2 + w1^2*x2 + 2*w1*w2*x2 + w2^2*x2)/(h^2 + w1^2 + 2*w1*w2 + w2^2))^2 + (h - h*(h^2 - h*r2 + w1*w2 - w1*x2 + w2^2 - w2*x2)/(h^2 + w1^2 + 2*w1*w2 + w2^2) - r2)^2,  # eq2
       -r1^2 + (x1 - (h^2*w2 - h*r1*w1 - h*r1*w2 + w1^2*x1 + 2*w1*w2*x1 + w2^2*x1)/(h^2 + w1^2 + 2*w1*w2 + w2^2))^2 + (-h*(h*r1 + w1*w2 - w1*x1 + w2^2 - w2*x1)/(h^2 + w1^2 + 2*w1*w2 + w2^2) + r1)^2,  # eq3
       -r2^2 + (x2 - (h*(h*w1 - h*x2 - r2*w1 + r2*w2) + x2*(h^2 + w1^2 - 2*w1*w2 + w2^2))/(h^2 + w1^2 - 2*w1*w2 + w2^2))^2 + (h - h*(h^2 - h*r2 - w1*w2 + w1*x2 + w2^2 - w2*x2)/(h^2 + w1^2 - 2*w1*w2 + w2^2) - r2)^2,  # eq4
   ]
   end;

(r1, r2, w1) = [7, 4, 9] .// 2
iniv = [big"8.1", 10, -3, 3]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[8.000000000000000289915981508997928978680584913763671176144672868231933897568114, 10.49999999999999484138720074297806292880734825561173528954255155770257311298397, -3.499999999999999683399841807629998190876460019780784537186673263801523533784222, 3.499999999999996870048545983071833654565390837799776925526657173046905812471916], true)

h = 8;  w2 = 10.5;  x1 = -3.5;  x2 = 3.5
下底の長さ = 2w2 = 20.99999999999999

算額の図は実際の数値に基づいていない。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (h, w2, x1, x2) = [9.0, 7.5, 2, -2.8]
   (h, w2, x1, x2) = res[1]
   @printf("h = %g;  w2 = %g;  x1 = %g;  x2 = %g\n", h, w2, x1, x2)
   println("下底の長さ = $(Float64(2w2))")
   plot([-w2, w2, w1, -w1, -w2], [0, 0, h, h, 0], color=:black, lw=0.5)
   circle(x1, r1, r1, :red)
   circle(x2, h - r2, r2, :green)
   segment(-w1, h, w2, 0)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(x1, r1, "乾円:r1(x1,r1)", :red, :center, :top, delta=-delta)
       point(x2, h - r2, "坤円:r2\n(x2,h-r2)", :green, :center, :bottom, delta=delta)
       point(w1, h, "(w1,h)", :black, :right, :bottom, delta=delta)
       point(w2, 0, "w2", :black, :left, :bottom, delta=delta)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その377）

算額（その377）

三六 川越市古谷本郷 古尾谷八幡神社 天保12年(1841)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

埼玉県川越市 古尾谷八幡神社 天保12年
http://tamtom.blog44.fc2.com/blog-entry-406.html

長方形内に，木，火，土，金，水の 5 円が入っている。木円の直径が 17 寸のとき，水円の直径を求めよ。

術に曰く

8 の開平方に 3 を加えたものを「極」と名付ける。
5 から 24 の開平方を引いたものに木円の径(17)を掛け，更に「極」を掛けることで水円の径を得る。

極 = √8 + 3
(5 - √24)*17*極 |> println
(5 - √24)*17*(√8 + 3) |> println

   10.009442010174723
   10.009442010174723

木円の半径と中心座標を r1, (x1, r1)
火円の半径と中心座標を r2, (x2, 2r1 - r2)
土円の半径と中心座標を r3, (x3, y3)
金円の半径と中心座標を r4, (r4, 2r1 - r4)
水円の半径と中心座標を r5, (r5, r5)
とする。

r1 = 17 として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms x1::positive, r1::positive, r2::positive,  x2::positive,
     r3::positive, x3::positive, y3::positive, r4::positive, r5::positive;

eq1 = (x1 - x2)^2 + (2r1 - r2 - r1)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq3 = (x1 - r5)^2 + (r1 - r5)^2 - (r1 + r5)^2
eq4 = (x2 - x3)^2 + (2r1 - r2 - y3)^2 - (r2 + r3)^2
eq5 = (x2 - r4)^2 + (r4 - r2)^2 - (r2 + r4)^2
eq6 = (x3 - r4)^2 + (2r1 - r4 - y3)^2 - (r3 + r4)^2
eq7 = (x3 - r5)^2 + (y3 - r5)^2 - (r3 + r5)^2
eq8 = (r5 - r4)^2 + (2r1 -r4 - r5)^2 - (r4 + r5)^2;

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

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)
   (x1, r2, x2, r3, x3, y3, r4, r5) = u
   return [
       (r1 - r2)^2 - (r1 + r2)^2 + (x1 - x2)^2,  # eq1
       (-r1 + y3)^2 - (r1 + r3)^2 + (x1 - x3)^2,  # eq2
       (r1 - r5)^2 - (r1 + r5)^2 + (-r5 + x1)^2,  # eq3
       -(r2 + r3)^2 + (x2 - x3)^2 + (2*r1 - r2 - y3)^2,  # eq4
       (-r2 + r4)^2 - (r2 + r4)^2 + (-r4 + x2)^2,  # eq5
       -(r3 + r4)^2 + (-r4 + x3)^2 + (2*r1 - r4 - y3)^2,  # eq6
       -(r3 + r5)^2 + (-r5 + x3)^2 + (-r5 + y3)^2,  # eq7
       (-r4 + r5)^2 - (r4 + r5)^2 + (2*r1 - r4 - r5)^2,  # eq8
   ]
end;

r1 = 17/2
iniv = [big"2.1", 0.25, 1, 0.2, 0.96, 1.3, 0.4, 0.6] .* r1
res = nls(H, ini=iniv);
println(res);

   (BigFloat[18.04927980072966600990743424177031286841046841344571292514916108734929792090129, 2.277568135664543005016640483057481026967779743580825713557355507836459659630661, 9.249432267243960090044883726521819586682889562263977714737446982268376356193549, 1.688090125635917522526267428949100510249937740184158307062589102275226843795516, 8.150820151636668050007516015640830020184999772933719726083627782869396869278436, 10.91198610762122951754068288011827079712431459719790662797327282314006745314288, 3.556929041116716190186087773344970300658574809456580227043502433088470316206073, 5.004721005087340180215705712235570136210476406875180438641404697785532644015275], true)

x1 = 18.0493;  r2 = 2.27757;  x2 = 9.24943;
r3 = 1.68809;  x3 = 8.15082;  y3 = 10.912;
r4 = 3.55693;  r5 = 5.00472
木円の直径 = 17.0;  水円の直径 = 10.00944201017468

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 17/2
   (x1, r2, x2, r3, x3, y3, r4, r5) = res[1]
   b = 2r1
   @printf("x1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  r5 = %g\n", x1, r2, x2, r3, x3, y3, r4, r5)
   println("木円の直径 = $(Float64(2r1));  水円の直径 = $(Float64(2r5))\n")
   plot([0, x1 + r1, x1 + r1, 0, 0], [0, 0, b, b, 0], color=:black, lw=0.5)
   circle(x1, r1, r1, :green)
   circle(x2, b - r2, r2, :red)
   circle(x3, y3, r3, :brown)
   circle(r4, b - r4, r4, :orange)
   circle(r5, r5, r5, :blue)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(x1, r1, "木円:r1,(x1,r1)", :green, :center, delta=-delta/2)
       point(x2, b - r2, "火円:r2\n(x2,b-r2)", :red, :center, :bottom, delta=delta/2)
       point(x3, y3, "土円:r3\n(x3,y3)", :brown, :center, :vcenter)
       point(r4, b - r4, "金円:r4\n(r4,b-r4)", :orange, :center, delta=-delta/2)
       point(r5, r5, "水円:r5,(r5,r5)", :blue, :center, delta=-delta/2)
       point(x1 + r1, 2r1, "(x1+r1,2r1)", :black, :right, delta=-delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;