算額（その785）

算額（その785）

福島県白河市南湖 南湖神社 昭和58年(1983)
http://www.wasan.jp/fukusima/nanko.html

半円内に 2 本の斜線を描き，天円，地円，人円を入れる。人円の直径が与えられたとき，天円の直径を求めよ。

半円の半径と中心座標を r1, (0, 0)
天円の半径と中心座標を r2, (0, r2)
地円の半径と中心座標を r3, (x3, r3)
人円の半径と中心座標を r4, (x4, r4)
斜線と半円の交点座標を (x01, y01), (x02, y02)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive,r3::positive, x3::positive,
     r4::positive, x4::positive, y4::positive,
     x01::positive, y01::positive, x02::positive, y02::positive
eq1 = x3^2 + (r2 - r3)^2 - (r2 + r3)^2
eq2 = x3^2 + r3^2 - (r1 - r3)^2
eq3 = x4^2 + y4^2 - (r1 - r4)^2
eq4 = 2sqrt(r1^2 - (r1 - 2r4)^2) - sqrt((x01 - x02)^2 + (y01 - y02)^2)
eq5 = x01^2 + y01^2 - r1^2
eq6 = x02^2 + y02^2 - r1^2
eq7 = (x01 - x02)/(y02 - y01) - y4/x4
eq8 = r1 - 2r2
eq9 = dist(x01, y01, x02, y02, 0, r2) - r2^2
eq10 = dist(x01, y01, x02, y02, x3, r3) - r3^2;

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=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (r1, r2, r3, x3, x4, y4, x01, y01, x02, y02) = u
   return [
       x3^2 + (r2 - r3)^2 - (r2 + r3)^2,  # eq1
       r3^2 + x3^2 - (r1 - r3)^2,  # eq2
       x4^2 + y4^2 - (r1 - r4)^2,  # eq3
       2*sqrt(r1^2 - (r1 - 2*r4)^2) - sqrt((x01 - x02)^2 + (y01 - y02)^2),  # eq4
       -r1^2 + x01^2 + y01^2,  # eq5
       -r1^2 + x02^2 + y02^2,  # eq6
       (x01 - x02)/(-y01 + y02) - y4/x4,  # eq7
       r1 - 2*r2,  # eq8
       -r2^2 + (-x01 - (-x01 + x02)*(-x01*(-x01 + x02) + (r2 - y01)*(-y01 + y02))/((-x01 + x02)^2 + (-y01 + y02)^2))^2 + (r2 - y01 - (-y01 + y02)*(-x01*(-x01 + x02) + (r2 - y01)*(-y01 + y02))/((-x01 + x02)^2 + (-y01 + y02)^2))^2,  # eq9
       -r3^2 + (r3 - y01 - (-y01 + y02)*((r3 - y01)*(-y01 + y02) + (-x01 + x02)*(-x01 + x3))/((-x01 + x02)^2 + (-y01 + y02)^2))^2 + (-x01 + x3 - (-x01 + x02)*((r3 - y01)*(-y01 + y02) + (-x01 + x02)*(-x01 + x3))/((-x01 + x02)^2 + (-y01 + y02)^2))^2,  # eq10
   ]
end;
r4 = 1/2
iniv = BigFloat[9.1, 4.6, 2.2, 6.3, 5.3, 6.6, 8.2, 3.6, 1.8, 8.8]
res = nls(H, ini=iniv)

   ([9.0, 4.5, 2.25, 6.363961030678928, 5.342584568965026, 6.611111111111111, 8.23517481947363, 3.630688046735422, 1.8214549574017131, 8.813756397709023], true)

人円の直径が 1 のとき，天円の直径は 9 である。

その他のパラメータは以下のとおりである。

r1 = 9;  r2 = 4.5;  r3 = 2.25;  x3 = 6.36396;  x4 = 5.34258;  y4 = 6.61111;  x01 = 8.23517;  y01 = 3.63069;  x02 = 1.82145;  y02 = 8.81376

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r4 = 1/2
   (r1, r2, r3, x3, x4, y4, x01, y01, x02, y02) = res[1]
   @printf("人円の直径が %g のとき，天円の直径は %g\n", 2r4, 2r2)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  x4 = %g;  y4 = %g;  x01 = %g;  y01 = %g;  x02 = %g;  y02 = %g\n",
       r1, r2, r3, x3, x4, y4, x01, y01, x02, y02)
   plot()
   circle(0, 0, r1, beginangle=0, endangle=180)
   circle(0, r2, r2, :blue)
   circle2(x3, r3, r3, :green)
   circle2(x4, y4, r4, :magenta)
   segment(x01, y01, x02, y02, :orange)
   segment(-x01, y01, -x02, y02, :orange)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(r1, 0, "r1  ", :red, :right, :bottom, delta=delta/2)
       point(0, 0, "半円:r1,(0,0)", :red, :center, :bottom, delta=delta/2)
       point(r1, 0, "r1  ", :red, :right, :bottom, delta=delta/2)
       point(0, r2, "天円:r2,(0,r2)", :blue, :center, :bottom, delta=delta)
       point(x3, r3, "地円:r3,(x3,r3)", :green, :center, :bottom, delta=delta)
       point(x4, y4, "人円:r4,(x4,r4) ", :black, :right, :vcenter)
       point(x01, y01, "(x01,y01)", :black, :right, :vcenter)
       point(x02, y02, "(x02,y02)", :black, :left, :bottom, delta=delta)
   end
end;

算額（その784）

算額（その784）

福島県白河市南湖 南湖神社 昭和58年(1983)
http://www.wasan.jp/fukusima/nanko.html

半円内に甲円を入れる。弧の両端から甲円に接する2直線を引き交点を頂点とする二等辺三角形を作る。斜辺に 2 点で内接し，同時に甲円に外接する乙円を入れる。乙円の直径が与えられたとき，半円の直径を求めよ。

半円の半径と中心座標を r1, (0, 0)
甲円の半径と中心座標を r2, (0, r2)
乙円の半径と中心座標を r3, (0, r3 + 2r2)
二等辺三角形の高さを h とする。
等辺のなす角度を θ とすると，sin(θ) = r1/sqrt(h^2 + r1^2) である。
以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms h::positive, r1::positive,
     r2::positive, r3::positive
sinθ = r1/sqrt(h^2 + r1^2)
eq1 = (h - r2)*sinθ - r2
eq2 = (h - r3 - 2r2)*sinθ - r3
eq3 = r1 - 2r2
res = solve([eq1, eq2, eq3], (h, r1, r2));
res[1] |> println

   (32*r3/3, 8*r3, 4*r3)

乙円の直径が 1 のとき，半円の直径は 8 である。

その他のパラメータは以下のとおりである。

r3 = 0.5;  h = 5.33333;  r1 = 4;  r2 = 2

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 1/2
   (h, r1, r2) = (32*r3/3, 8*r3, 4*r3)
   @printf("乙円の直径が %g のとき，半円の直径は %g\n", 2r3, 2r1)
   @printf("r3 = %g;  h = %g;  r1 = %g;  r2 = %g\n",  r3, h, r1, r2)
   plot([r1, 0, -r1, r1], [0, h, 0, 0], color=:black, lw=0.5)
   circle(0, 0, r1, beginangle=0, endangle=180)
   circle(0, r2, r2, :blue)
   circle(0, r3 + 2r2, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, h, " h", :red, :left, :bottom, delta=delta/2)
       point(0, 0, "半円:r1,(0,0)", :red, :center, :bottom, delta=delta/2)
       point(r1, 0, "r1  ", :red, :right, :bottom, delta=delta/2)
       point(0, r2, "甲円:r2 \n(0,r2)", :blue, :right, :bottom, delta=delta/2)
       point(0, r3 + 2r2, " 乙円:r3,(0,r3+2r2)", :green, :left, :bottom, delta=delta/2)
   end
end;

算額（その783）

算額（その783）

百四　岩手県大船渡市田茂山 根城八幡宮 天保11年(1840)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html

福島県白河市南湖 南湖神社 昭和58年(1983)
http://www.wasan.jp/fukusima/nanko.html

和算に挑戦_一関市博物館https://www.city.ichinoseki.iwate.jp/museum/wasan/h25/hard.html

外円内に正三角形と甲円，乙円を入れる。乙円の直径が与えられたとき，甲円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1,(x1, r1 - R/2)
乙円の半径と中心座標を r2, (x2, r2 - R/2), (x3, r2 - R/2), (0, r2 - R); x3 = x2 - 2r2
とおいて以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, x1::positive,
     r2::positive, x2::negative, x3::negative
R = 4r2
x3 = x2 - 2r2
eq1 = x1^2 + (r1 - R/2)^2- (R - r1)^2
eq2 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq3 = x3^2 + (r2 - R/2)^2 - (R - r2)^2
solve([eq1, eq2, eq3], (r1, x1, x2))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (2*r2, 2*r2, 2*r2*(1 - sqrt(2)))

乙円の直径が 1 のとき，甲円の直径は 2 である。
なお，甲円の中心は x 軸上にある。また，横並びの 2 個の乙円は共に x 軸に接している。

その他のパラメータは以下のとおりである。

   r2 = 0.5;  r1 = 1;  x1 = 1;  x2 = -0.414214;  x3 = -1.41421;  R = 2

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (r1, x1, x2) = (2*r2, 2*r2, 2*r2*(1 - sqrt(2)))
   R = 4r2
   x3 = x2 - 2r2
   @printf("r2 = %g;  r1 = %g;  x1 = %g;  x2 = %g;  x3 = %g;  R = %g\n",  r2, r1, x1, x2, x3, R)
   @printf("乙円の直径が %g のとき，甲円の直径は %g\n", 2r2, 2r1)
   plot([√3R/2, 0, -√3R/2, √3R/2], [-R/2, R, -R/2, -R/2], color=:black, lw=0.5)
   circle(0, 0, R)
   circle(x1, r1 - R/2, r1, :blue)
   circle(x2, r2 - R/2, r2, :green)
   circle(x3, r2 - R/2, r2, :green)
   circle(0, r2 - R, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, R, " R", :red, :left, :bottom, delta=delta/2)
       point(x1, r1 - R/2, "甲円:r1 \n(x1,r1-R/2)", :blue, :right, :bottom, delta=delta/2)
       point(x2, r2 - R/2, "乙円:r2\n(x2,r2-R/2)", :green, :center, :bottom, delta=delta/2)
       point(x3, r2 - R/2, "乙円:r2\n(x3,r2-R/2)", :green, :left, delta=-delta/2)
       point(0, r2 - R, "乙円:r2\n(0,r2-R)", :green, :center, delta=-delta/2)
       point(√3R/2, -R/2, "(√3R/2,-R/2) ", :black, :right, delta=-delta/2)
   end
end;

算額（その782）

算額（その782）

福島県白河市南湖 南湖神社 昭和58年(1983)
http://www.wasan.jp/fukusima/nanko.html

直角三角形の各辺を直径とする 3 個の半円と，それらの隙間に乾円 1 個，坤円 2 個を入れる。坤円の直径が与えられたとき，乾円の直径を求めよ。

直角三角形の直角を挟む二辺のうち短い方の長さを 2r1，長い方の長さを 2r2，対辺の長さを 2R とする。
乾円の半径と中心座標を r3, (r3, y3)
坤円の半径と中心座標を r4, (x4, y4), (r2, -r4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, r2::positive, r3::positive, y3::positive,
     r4::positive, x4::positive, y4::positive, x0::positive, y0::positive, d
R = sqrt(4r1^2 + 4r2^2)/2
eq1 = r1 + 2r4 - R
eq2 = x0^2 + (y0 - r1)^2 - r1^2
eq3 = y0/(2r2 - x0) - 2r1/2r2
eq4 = x4^2 + (r1 - y4)^2 - (r1 - r4)^2
eq5 = (r2 - r3)^2 + y3^2 - (r2 + r3)^2
eq6 = (r2 - x4)^2 + y4^2 - (r2 - r4)^2
eq7 = dist(0, 2r1, 2r2, 0, r3, y3) - r3^2
eq7 = div(numerator(apart(eq7, d)), r2);
eq8 = 2R*sqrt(x0^2 + y0^2) - 4r1*r2;

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=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (r1, r2, r3, y3, x4, y4, x0, y0) = u
   return [
       r1 + 2*r4 - sqrt(4*r1^2 + 4*r2^2)/2,  # eq1
       -r1^2 + x0^2 + (-r1 + y0)^2,  # eq2
       -r1/r2 + y0/(2*r2 - x0),  # eq3
       x4^2 - (r1 - r4)^2 + (r1 - y4)^2,  # eq4
       y3^2 + (r2 - r3)^2 - (r2 + r3)^2,  # eq5
       y4^2 - (r2 - r4)^2 + (r2 - x4)^2,  # eq6
       4*r1^2*r2 - 4*r1^2*r3 - 4*r1*r2*y3 + 2*r1*r3*y3 - r2*r3^2 + r2*y3^2,  # eq7
       -4*r1*r2 + sqrt(4*r1^2 + 4*r2^2)*sqrt(x0^2 + y0^2),  # eq8
   ]
end;

r4 = 1/2
iniv = BigFloat[1.5, 2, 0.5, 2, 0.8, 0.9, 36/25, 48/25]
res = nls(H, ini=iniv);
println(res);

   ([1.5, 2.0, 0.5, 2.0, 0.8000000000009287, 0.9000000000012384, 1.44, 1.92], true)

坤円の半径が 0.5 のとき，乾円の半径は 0.5 である（坤円の半径に等しい）。

その他のパラメータは以下のとおりである。

r1 = 1.5;  r2 = 2;  r3 = 0.5;  y3 = 2;  x4 = 0.8;  y4 = 0.9;  x0 = 1.44;  y0 = 1.92;  R = 2.5

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r4 = 0.5
   (r1, r2, r3, y3, x4, y4, x0, y0) = (1.5, 2, 0.5, 2, 0.8, 0.9, 36/25, 48/25)
   (r1, r2, r3, y3, x4, y4, x0, y0) = res[1]
   R = sqrt(4r1^2 + 4r2^2)/2
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  y3 = %g;  x4 = %g;  y4 = %g;  x0 = %g;  y0 = %g;  R = %g\n",  r1, r2, r3, y3, x4, y4, x0, y0, R)
   @printf("乾円の直径 = %g;  坤円の直径 = %g\n", r4, r3)
   plot([0, 2r2, 0, 0], [0, 0, 2r1, 0], color=:black, lw=0.5)
   circle(0, r1, r1, beginangle=270, endangle=450)
   circle(r2, 0, r2, :blue, beginangle=0, endangle=180)
   circle(r3, y3, r3, :green)
   circle(x4, y4, r4, :magenta)
   circle(r2, -r4, r4, :magenta)
   circle(r2, r1, R, :orange)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(x0, y0, " (x0,y0)", :black, :left, :vcenter)
       point(r2, r1, " (r2,r1)", :black, :left, :vcenter)
       point(r3, y3, "乾円:r3\n(r3,y3)", :green, :center, delta=-delta/2)
       point(x4, y4, "坤円:r4\n(x4,y4)", :magenta, :center, delta=-delta/2)
       point(r2, -r4, "坤円:r4\n(r2,-r4)", :magenta, :center, delta=-delta/2)
       point(0, 2r1, "2r1 ", :red, :right, :vcenter)
       point(0, r1, "r1 ", :red, :right, :vcenter)
       point(r2, 0, " r2", :blue, :left, :bottom, delta=delta/2)
       point(2r2, 0, " 2r2", :blue, :left, :bottom, delta=delta/2)
   else
       plot!(showaxis=false)
   end
end;