算額（その355）

算額（その355）

福島県楢葉町北田 北田神社 明治27年(1894)
ならはの絵馬―村人の祈り― 018/036page
http://is2.sss.fukushima-u.ac.jp/fks-db/txt/10088.002/html/00018.html
福島県楢葉町北田には，北田天満宮と大山祇神社があるが北田神社はないようだ。

正三角形内に松円，松円の半円，竹円，梅円を入れる。梅円直径が 3 寸のとき，松円の直径を求めよ。

正三角形の一辺の長さを 1 とする。
松円の半径，中心座標を r1, (0, y1)
竹円の半径，中心座標を r2, (0, r2)
梅円の半径，中心座標を r3, (r3, y3)
とおき，以下の連立方程の解を求める。

include("julia-source.txt");

using SymPy

@syms r1::positive, y1::positive, r2::positive, r3::positive, x3::positive;

eq1 = (1 - r1/sqrt(Sym(3)))^2 + (y1 - r1)^2 - 4r1^2
eq2 = (1 - r1/sqrt(Sym(3)))^2 + (r1 - r2)^2 - (r1 + r2)^2
eq3 = (1 - r1/sqrt(Sym(3)) - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq4 = x3^2 + (r2 - r3)^2 - (r2 + r3)^2
eq5 = (sqrt(Sym(3)) - y1)/2 - r1;
res = solve([eq1, eq2, eq3, eq4, eq5], (r1, y1, r2, r3, x3))

   1-element Vector{NTuple{5, Sym}}:
    (sqrt(3)/4, sqrt(3)/2, 3*sqrt(3)/16, -9 + 21*sqrt(3)/4, -9/4 + 3*sqrt(3)/2)

res[1][1] |> println  # r1: 松円の半径
res[1][4] |> println  # r3: 梅円の半径

   sqrt(3)/4
   -9 + 21*sqrt(3)/4

res[1][1]/res[1][4] |> simplify |> println

   4*sqrt(3)/3 + 7/3

松円の半径（直径）は梅円の半径（直径）の (4√3 + 7)/3 倍である。
梅円の直径が 3 寸ならば，松円の直径は 4√3 + 7 寸である。

4√3 + 7

   13.928203230275509

   r1 = 0.433013; y1 = 0.866025; r2 = 0.32476; r3 = 0.0932667; x3 = 0.348076
   梅円が 3 寸のとき，松円は r1/r3 * 3 = 13.9282 寸

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, y1, r2, r3, x3) = (sqrt(3)/4, sqrt(3)/2, 3*sqrt(3)/16, -9 + 21*sqrt(3)/4, -9/4 + 3*sqrt(3)/2)
   @printf("r1 = %g; y1 = %g; r2 = %g; r3 = %g; x3 = %g\n", r1, y1, r2, r3, x3)
   @printf("梅円が 3 寸のとき，松円は r1/r3 * 3 = %g 寸\n", 3r1/r3)
   plot([1, 0, -1, 1], [0, sqrt(3), 0, 0], color=:black, lw=0.5)
   circle(0, y1, r1, :green)
   circle(1 - r1/sqrt(3), r1, r1, :green, beginangle=120, endangle=300)
   circle(r1/sqrt(3) - 1, r1, r1, :green, beginangle=240, endangle=420)
   circle(0, r2, r2, :blue)
   circle(x3, r3, r3, :red)
   circle(-x3, r3, r3, :red)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, y1, " y1:\n 松円:r1", :green, :left, :vcenter)
       point(1 - r1/sqrt(3), r1, "(1-r1/√3,r1)", :green, :right, :vcenter)
       point(0, r2, " r2:\n 竹円:r2", :blue, :left, :vcenter)
       point(x3, r3, " 梅円:r3\n (x3,r3)", :red, :left, :vcenter)
       point(0, √3, " √3", :black)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その354）

算額（その354）

福島県楢葉町北田 北田神社 明治27年(1894)
ならはの絵馬―村人の祈り― 018/036page
http://is2.sss.fukushima-u.ac.jp/fks-db/txt/10088.002/html/00018.html
福島県楢葉町北田には，北田天満宮と大山祇神社があるが北田神社はないようだ。

外円内に春夏秋冬の6円を入れる。春円，夏円の直径をそれぞれ 3 寸，4 寸としたとき，秋円の直径を求めよ。

外円の半径，中心座標を r0, (0, 0)
春円の半径，中心座標を r1, (0, r0 - r1)
夏円の半径，中心座標を r2, (x2, y2)
秋円の半径，中心座標を r3, (r3, y3)
冬円の半径，中心座標を r4, (0, r0 - 2r1 - r4)
とおき，以下の連立方程式の数値解を nlsolve() で求める。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive, x2::positive,
     y2::positive, r3::positive, y3::negative, r4::positive;

(r1, r2) = (3, 4) .// 2
eq1 = x2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + (r0 - 2r1 - r4 - y2)^2 - (r2 + r4)^2
eq3 = x2^2 + y2^2 - (r0 - r2)^2
eq4 = (x2 - r3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq5 = r3^2 + y3^2 - (r0 - r3)^2
eq6 = r3^2 + (r0 - 2r1 - r4 - y3)^2 - (r3 + r4)^2;
# solve([eq1, eq2, eq3, eq4, eq5, 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)
   (r0, x2, y2, r3, y3, r4) = u
   return [
       x2^2 + (r0 - y2 - 3/2)^2 - 49/4,  # eq1
       x2^2 - (r4 + 2)^2 + (r0 - r4 - y2 - 3)^2,  # eq2
       x2^2 + y2^2 - (r0 - 2)^2,  # eq3
       (-r3 + x2)^2 - (r3 + 2)^2 + (y2 - y3)^2,  # eq4
       r3^2 + y3^2 - (r0 - r3)^2,  # eq5
       r3^2 - (r3 + r4)^2 + (r0 - r4 - y3 - 3)^2,  # eq6
   ]
end;

iniv = [28.0, 7.2, 8.0, 6.0, -5.2, 2.8]
res = nls(H, ini=iniv);

   r0 = 8.69972; x2 = 3.23607; y2 = 5.86635; r3 = 4.34164; y3 = -0.378172; r4 = 1.77267

秋円の直径は 8.68328 寸である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (3, 4) .// 2
   (r0, x2, y2, r3, y3, r4) = res[1]
   @printf("r0 = %g; x2 = %g; y2 = %g; r3 = %g; y3 = %g; r4 = %g\n", r0, x2, y2, r3, y3, r4)
   @printf("秋円直径 = %g 寸\n", 2r3)
   plot()
   circle(0, 0, r0, :black)
   circle(0, r0 - r1, r1, :green)
   circle(x2, y2, r2)
   circle(-x2, y2, r2)
   circle(r3, y3, r3, :gray)
   circle(-r3, y3, r3, :gray)
   circle(0, r0 - 2r1 - r4, r4, :orange)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, r0 - r1, "春:r1\nr0-r1", :green, :left, :vcenter)
       point(x2, y2, "夏:r2,(x2,y2)", :red, :center, delta=-delta)
       point(0, r0 - 2r1 - r4, "冬:r4\nr0-2r1-r4", :orange, :left, :vcenter)
       point(r3, y3, "秋:r3,(r3,y3)", :gray, :center, delta=-delta)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その353）

算額（その353）

群馬県高崎市山名町 山名八幡宮 文化11年(1814)
http://takasakiwasan.web.fc2.com/gunnsann/g029-4.html

鉤股弦（直角三角形）内に，大円 1 個と小円 3 個を入れる。鉤の長さが 3 寸のとき，股の長さはいかほどか。

include("julia-source.txt");

using SymPy

@syms 鉤::positive, 股::positive, x1::positive, r1::positive,
     r2::positive, x2::positive, y2::positive;

鉤 = 3
eq1 = (x1 - r2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x1 - r2)^2 + (3r2 - r1)^2 - (r1 + r2)^2
eq3 = (x1 - x2)^2 + (y2 - r1)^2 - (r1 + r2)^2
eq4 = (x2 - r2)^2 + (y2 - 3r2)^2 - 4r2^2
eq5 = distance(0, 鉤, 股, 0, x1, r1) - r1^2
eq6 = distance(0, 鉤, 股, 0, x2, y2) - r2^2;
# solve([eq1, eq2, eq3, eq4, eq5, 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, x1, r2, x2, y2) = u
   return [
       (r1 - r2)^2 - (r1 + r2)^2 + (-r2 + x1)^2,  # eq1
       (-r1 + 3*r2)^2 - (r1 + r2)^2 + (-r2 + x1)^2,  # eq2
       (-r1 + y2)^2 - (r1 + r2)^2 + (x1 - x2)^2,  # eq3
       -4*r2^2 + (-3*r2 + y2)^2 + (-r2 + x2)^2,  # eq4
       -r1^2 + (r1 - 3*(3*r1 - x1*股 + 股^2)/(股^2 + 9))^2 + (x1 - 股*(-3*r1 + x1*股 + 9)/(股^2 + 9))^2,  # eq5
       -r2^2 + (x2 - 股*(x2*股 - 3*y2 + 9)/(股^2 + 9))^2 + (y2 - 3*(-x2*股 + 3*y2 + 股^2)/(股^2 + 9))^2,  # eq6
   ]
end;

iniv = [big"4.0", 1, 2, 0.5, 1.2, 2]
res = nls(H, ini=iniv);

   股 = 3.71231;  r1 = 0.782752;  x1 = 1.49835;  r2 = 0.391376;  x2 = 0.883366;  y2 = 1.78293

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   鉤 = 3
   (股, r1, x1, r2, x2, y2) = res[1]
   @printf("股 = %g;  r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", 股, r1, x1, r2, x2, y2)
   plot([0, 股, 0, 0], [0, 0, 鉤, 0], color=:black, lw=0.5)
   circle(x1, r1, r1, :blue)
   circle(r2, r2, r2)
   circle(r2, 3r2, r2)
   circle(x2, y2, r2)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(x1, r1, "大円:r1,(x1, r1)", :blue, :center, delta=-delta)
       point(r2, r2, "小円:r2\n(r2,r2)", :green, :center, delta=-delta)
       point(r2, 3r2, "(r2,3r2)", :red, :center, delta=-delta)
       point(x2, y2, "(x2,x2)", :red, :center, delta=-delta)
       point(股, 0, " 股", :black, :left, :bottom)
       point(0, 鉤, " 鉤", :black, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その352）

算額（その352）

群馬県桐生市天神町 天満宮 文化6年(1830)

大山誠: 桐生市天満宮の算額題についてhttps://www.jstage.jst.go.jp/article/jjsme/80/7/80_20/_pdf

図のように大円と直角三角形が交わっており，等円を 3 個配置する。大円の直径が 10 寸のとき，等円の直径を求めよ。

大円の直径を r0 = 10 ，AC の長さを a，B の x 座標を b（b < 0)，等円の直径を r，その中心座標を (r0 - r, 0), (x1, r - a), (x2, y2) とおく。

以下の連立方程式を nlsolve() で解く。
eq6 は，∠EBF = θ とすると，∠OAD = π/4 - θ および tan(π/4 - θ) = (1 - tan(θ)/(1 + tan(θ) であることによる。

include("julia-source.txt");

using SymPy

@syms r0::positive, r::positive, a::positive, b::negative, x1::negative, x2::negative, y2::positive;

r0 = 10
eq1 = (r0 - 2r)^2 + a^2 - r0^2
eq2 = x2^2 + y2^2 - (r0 - r)^2
eq3 = x1^2 + (r - a)^2 - (r0 + r)^2
eq4 = distance(b, -a, r0 - 2r, a, x1, r - a) - r^2
eq5 = distance(b, -a, r0 - 2r, a, x2, y2) - r^2
tanθ = r/(x1 - b)
eq6 = (r0 - 2r)/a - (1 - tanθ)/(1 + tanθ);  # (r0 - 2r)/a - tan(PI/4 - θ)
# θ =atan(tanθ)
# eq6 = (r0 - 2r)/a - tan(PI/4 - θ);
# res = solve([eq1, eq2, eq3, eq4, eq5], (r, a, b, x1, x2, y2))

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)
   (r, a, b, x1, x2, y2) = u
   return [
       a^2 + (10 - 2*r)^2 - 100,  # eq1
       x2^2 + y2^2 - (10 - r)^2,  # eq2
       x1^2 + (-a + r)^2 - (r + 10)^2,  # eq3
       -r^2 + (x1 - (4*a^2*b - 2*a*b*r - 4*a*r^2 + 20*a*r + b^2*x1 + 4*b*r*x1 - 20*b*x1 + 4*r^2*x1 - 40*r*x1 + 100*x1)/(4*a^2 + b^2 + 4*b*r - 20*b + 4*r^2 - 40*r + 100))^2 + (-a - a*(-4*a^2 + 4*a*r + b^2 - 2*b*x1 - 4*r^2 - 4*r*x1 + 40*r + 20*x1 - 100)/(4*a^2 + b^2 + 4*b*r - 20*b + 4*r^2 - 40*r + 100) + r)^2,  # eq4
       -r^2 + (x2 - (b*(4*a^2 + b^2 + 4*b*r - 20*b + 4*r^2 - 40*r + 100) - (b + 2*r - 10)*(2*a^2 + 2*a*y2 + b^2 + 2*b*r - b*x2 - 10*b - 2*r*x2 + 10*x2))/(4*a^2 + b^2 + 4*b*r - 20*b + 4*r^2 - 40*r + 100))^2 + (-a*(4*a*y2 + b^2 - 2*b*x2 - 4*r^2 - 4*r*x2 + 40*r + 20*x2 - 100)/(4*a^2 + b^2 + 4*b*r - 20*b + 4*r^2 - 40*r + 100) + y2)^2,  # eq5
       -(-r/(-b + x1) + 1)/(r/(-b + x1) + 1) + (10 - 2*r)/a,  # eq6
   ]
end;

iniv = [big"2.8", 8.9, -19.5, -11.2, -4.3, 5.8]
res = nls(H, ini=iniv);
println(res);
     (BigFloat[2.758765525590156335644381495793494653227916752992539435376878299670476979753896, 8.939097947940123521918150817467879027312574890368149878678356808089001454295962, -19.47025838954737912975466155729263249608233106765657839605148838383817497557, -11.16197065424549685348050573396238073769521442013558168263557635706947113807764, -4.33134044315998946461745333303927732461882711841070523331857400641089934088305, 5.803013585959302297582280566300900177092654477073446038766027387779471835660339], true)

   r = 2.75877;  a = 8.9391;  b = -19.4703;  x1 = -11.162;  x2 = -4.33134;  y2 = 5.80301

大円の直径が 10 寸のとき，等円の直径は 2.75877 寸である。

include("julia-source.txt");

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 10
   (r, a, b, x1, x2, y2) = res[1]
   @printf("r = %g;  a = %g;  b = %g;  x1 = %g;  x2 = %g;  y2 = %g\n", r, a, b, x1, x2, y2)
   plot()
   circle(0, 0, r0, :blue)
   circle(x1, r - a, r)
   circle(x2, y2, r)
   circle(r0 - r, 0, r)

   if more
       plot!([r0 - 2r, r0 - 2r, b], [-a, a, -a], color=:black, lw=0.5)
       segment(0, 0, r0 - 2r, a, :orange)
       segment(b, -a, x1, r - a, :orange)
       segment(x1, r - a, x1, -a, :orange)
       point(r0, 0, " r0", :blue, :left, :bottom)
       point(r0 - r, 0, "r0-r", :red, :center, :bottom)
       point(r0 - 2r, 0, "r0-2r ", :red, :right, :bottom)
       point(r0 - 2r, a, " A:(r0-2r,a)", :black, :left, :bottom)
       point(b, -a, " B:(b,-a)", :black)
       point(r0 - 2r, -a, " C:(r0-2r,-a)", :black, :left)
       point(r0 - 2r, 0, "D ", :black, :right)
       point(0, 0, " O", :black)
       point(x1, r - a, "E:(x1,r-a)", :red, :center, :bottom, mark=false)
       point(x1, -a, " F", :black)
       point(x2, y2, "(x2,y2)", :red, :center)
       hline!([0, -a], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!([b, r0 - 2r, r0 - 2r, b], [-a, -a, a, -a], color=:black, lw=0.5)
       plot!(showaxis=false)
   end
end;

注: r の解析解は以下の 3次式=0 の解の 10 倍である。
SymPy では虚数解になるが虚部は非常に小さく（実質0），3 個の実数解が得られる。この内適解は t = 0.275876552559016 なので，r = 2.75876552559016 である。

using SymPy
@syms t::positive
eq = 8t^3 - 8sqrt(Sym(2))t^2 +5t - (12 - 8sqrt(Sym(2)))
res2 = solve(eq)
res2[1].evalf() |> println
res2[2].evalf() |> println
res2[3].evalf() |> println

   0.45518044897151 + 0.e-21*I
   0.275876552559016 + 0.e-24*I
   0.683156560842569 - 0.e-22*I

f(t) = 8t^3 - 8sqrt(Sym(2))t^2 +5t - (12 - 8sqrt(Sym(2))) を図にすると以下のようになる。

using Plots
pyplot(size=(500, 500), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(eq, xlims=(0.25, 0.7))
hline!([0])
vline!([0.275876552559016, 0.45518044897151, 0.683156560842569])

算額（その351）

算額（その351）

岐阜県大垣市赤坂町 金生山明星輪寺 元治2年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

正三角形と外円との間に，緑円，赤円，白円を計 16 個をいれる。正三角形の内接円の半径が与えられたとき白円の直径を求めよ。

記述を簡潔にできるように元の図を上下反転させて考える。
緑円の半径と中心座標 r1, (x, y)
赤円の半径と中心座標 r2, (0, 0), (0, 2r2 + 4r3), (0, 4r2 + 4r3)
白円の半径と中心座標 r3, (0, r2 + r3), (0, r2 + 3r3)
外円の半径と中心座標 r0, (0, 0);  r0 = 5r2 + 4r3
正三角形に内接する円の半径は 3r2 + 4r3

外円の半径を 1 として他の変数の値を求め，r2, r3 を求め，r2 と 3r2 + 4r3 の比を求める。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive, r3::positive, x::positive, y::positive;

r0 = 1
y = r2 + 2r3
r1 = 2r2 + 2r3
eq1 = x^2 + y^2 - (r1 + r2)^2
eq2 = x^2 + (y - r2 - r3)^2 - (r1 + r3)^2
eq3 = 5r2 + 4r3 - r0

res = solve([eq1, eq2, eq3], (r2, r3, x))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (1/7, 1/14, 2*sqrt(3)/7)

   r0 = 1; r1 = 0.428571; r2 = 0.142857; r3 = 0.0714286; x = 0.494872, y = 0.285714
   正三角形に内接する円の半径 = 0.714286
   白円の半径 = 0.0714286

赤円，白円の半径はそれぞれ 1/7, 1/14。
正三角形の内接円の半径は 3r2 + 4r3 = 3/7 + 4/14 = 10/14。
よって，白円の半径は正三角形の内接円の半径の 1/10 である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3, x) = res[1]
   r0 = 1
   y = r2 + 2r3
   r1 = 2r2 + 2r3
   @printf("r0 = %g; r1 = %g; r2 = %g; r3 = %g; x = %g, y = %g\n", r0, r1, r2, r3, x, y)
   @printf("正三角形に内接する円の半径 = %g\n", 3r2 + 4r3)
   @printf("白円の半径 = %g\n", r3)
   plot()
   circle(0, 0, r0, :black)
   circle(0, 0, r2)
   rotate(x, y, r1, :green)
   rotate(0, r2 + r3, r3, :blue)
   rotate(0, r2 + 3r3, r3, :blue)
   rotate(0, 2r2 + 4r3, r2)
   rotate(0, 4r2 + 4r3, r2)
   ic = 3r2 + 4r3
   plot!([√3ic, -√3ic, 0, √3ic], [ic, ic, -2ic, ic], color=:black, lw=0.5)
   if more
       point(x, y, "(x,y)")
       point(0, r2 + r3, "  r2+r3", :blue, :left, :vcenter)
       point(0, r2 + 3r3, "  r2+3r3", :blue, :left, :vcenter)
       point(0, 2r2 + 4r3, " 2r2+4r3", :red, :left, :vcenter)
       point(0, 3r2 + 4r3, "   3r2+4r3", :red, :left, :bottom)
       point(0, 4r2 + 4r3, " 4r2+4r3", :red, :left, :vcenter)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その350）

算額（その350）

岐阜県大垣市赤坂町 金生山明星輪寺 元治2年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

外円内に青，緑，赤，白の 12 個の円をそれぞれ接するように入れる。外円の直径が与えられたときに白円の直径を求めよ。

2 つの青円の中心を y 軸，2 つの緑円の中心を x 軸に置く。
x 軸，y 軸で点対称なので，第1象限の図形を対象とする。

外円の半径と中心座標 r0, (0, 0)
青円の半径と中心座標 r1, (0, r1)
緑円の半径と中心座標 r2, (x2, 0)
赤円の半径と中心座標 r3, (x3, y3)
白円の半径と中心座標 r4, (x4, y4)

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

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive, x2::positive, r3::positive,
     x3::positive, y3::positive, r4::positive, x4::positive, y4::positive;

r0 = 1
r1 = r0//2
eq1 = x2^2 + (r0 - r1)^2 - (r1 + r2)^2 |> expand
eq2 = x3^2 + (r0 - r1 - y3)^2 - (r1 + r3)^2 |> expand
eq3 = (x3 - x2)^2 + y3^2 - (r2 + r3)^2 |> expand
eq4 = (x4 - x2)^2 + r4^2 - (r2 + r4)^2 |> expand
eq5 = (x4 - x3)^2 + (y3 - r4)^2 - (r3 + r4)^2 |> expand
eq6 = x3^2 + y3^2 - (r0 - r3)^2 |> expand
eq7 = x4^2 + r4^2 - (r0 - r4)^2 |> expand;

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7]) #, (r2, x2, r3, x3, y3, r4, x4))

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

iniv = [big"0.5", 0.5, 0.5, 0.5, 0.5, 0.5, 0.5]
res = nls(H, ini=iniv);
println(res);

 (BigFloat[0.2307692307692307945960496972420115643146885404357311695540829279753122088871492, 0.5329387100211930608365475760066592887331159643976213583701227822622771799660959, 0.1999999999999999665985768137693666773172486596065619280957172091107786570697997, 0.6928203230275509687226737223425773030846798934032123402136212243156691574043034, 0.400000000000000100204269558691899968048254021180314215712850083601413507819824, 0.1249999999999999591818515058195745225114059751553107390846466073349432317897547, 0.8660254037844387219680669190319846379474904384664942387376085192902656644536619], true)

外円の半径が 1 のとき，緑円の半径 = 0.230769;  赤円の半径 = 0.2;  白円の半径 = 0.125

0.230769 = 3/13

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 1
   r1 = 1//2
   (r2, x2, r3, x3, y3, r4, x4) = res[1]
   @printf("外円の半径が %g のとき，緑円の半径 = %g;  赤円の半径 = %g;  白円の半径 = %g\n", r0, r2, r3, r4)
   plot()
   circle(0, 0, r0, :black)
   circle(0, r0 - r1, r1, :blue)
   circle(0, r1 - r0, r1, :blue)
   circle(x2, 0, r2, :green)
   circle(-x2, 0, r2, :green)
   circle4(x3, y3, r3)
   circle4(x4, r4, r4, :gray)
   if more
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その349）

算額（その349）

岐阜県大垣市赤坂町 金生山明星輪寺 元治2年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

台形内に青の半円と，赤，白，黄，浅青，緑の合計 7 個の円を入れる。赤円の直径が与えられたときに，緑円の直径を求めよ。

台形の下底（図の右側の辺），上底の長さをそれぞれ a, b とする。青円の中心を座標原点とする。。
青円の半径と中心座標 r1, (0, 0)
赤円の半径と中心座標 r2, (r2, r2)
黃円の半径と中心座標 r4, (0, r1 - r4)
緑円の半径と中心座標 r6, (x6, y6)
白円の半径と中心座標 r3, (r1 - r3, y3)
浅青円の半径と中心座標 r5, (r5 - r1, y5)

r2 を既知とすれば，残りの未知数は 11 変数なので，11 元連立方程式になるが，条件は 10 個しかない。a, b は確定できないので，b も与えないと解けない。いずれにせよ solve() では一度には解けないので，白円と浅青円を除いてまず解を求め，しかる後に得られた解を既知として描画に必要なパラメータを求める。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive, y3::positive, r4::positive,
     r5::positive, y5::positive, r6::positive, x6::positive, y6::positive;

r2 = 1
eq1 = 2r2^2 - (r1 - r2)^2
eq3 = r2^2 + (r1 - r4 - r2)^2 - (r2 + r4)^2
eq5 = (x6 - r2)^2 + (y6 - r2)^2 - (r2 + r6)^2
eq6 = x6^2 + (y6 - r1 + r4)^2 - (r4 + r6)^2
eq7 = x6^2 + y6^2 - (r1 - r6)^2

res = solve([eq1, eq3, eq5, eq6, eq7], (r1, r4, r6, x6, y6))

   1-element Vector{NTuple{5, Sym}}:
    (1 + sqrt(2), -1 + sqrt(2), 1/5, 3/5, 1 + 4*sqrt(2)/5)

赤円の半径を 1 としたとき，緑円の半径は 1/5 である（赤円の直径の 1/5 である）。

残る 6 変数 r3, y3, r5, y5, a, b に対して，条件式は 5 個しかない。
図を見ればわかるが，a と b は互いに従属しており，a が小さくなれば b は大きくなる。
ここでは，b = 2 として解を求める。

using SymPy

@syms r1::positive, r2::positive, r3::positive, y3::positive, r4::positive,
     r5::positive, y5::positive, r6::positive, x6::positive, y6::positive,
     a::positive, b::positive;

r2 = 1
b = 2
(r1, r4, r6, x6, y6) = (1 + sqrt(Sym(2)), -1 + sqrt(Sym(2)), 1//5, 3//5, 1 + 4*sqrt(Sym(2))/5)
eq2 = (r1 - r3)^2 + y3^2 - (r1 + r3)^2
eq4 = (r1 - r5)^2 + y5^2 - (r1 + r5)^2
eq8 = r3/(a - y3) - r1/a
eq9 = r5/(b - y5) - r1/b
eq10 = 4r1^2 + (a - b)^2 - (a + b)^2
res = solve([eq2, eq4, eq8, eq9, eq10], (r3, y3, r5, y5, a))

   1-element Vector{NTuple{5, Sym}}:
    (-66932*sqrt(236*sqrt(2) + 334) - 94628*sqrt(118*sqrt(2) + 167) - 7 + 9*sqrt(2) + 66912*sqrt(2)*sqrt(118*sqrt(2) + 167) + 47328*sqrt(2)*sqrt(236*sqrt(2) + 334), -4 + 2*sqrt(118*sqrt(2) + 167)/(2*sqrt(2) + 3), (-sqrt(2)/2 - 1/2)*(-3 - 2*sqrt(2) + sqrt(20*sqrt(2) + 29)) + 1 + sqrt(2), -3 - 2*sqrt(2) + sqrt(20*sqrt(2) + 29), sqrt(2) + 3/2)

二段階に分けて求めたパラメータで図を描き，解が正しいことを確認する。

   r1 = 2.41421;  r2 = 1;  r4 = 0.414214
   r6 = 0.2;  x6 = 0.6;  y6 = 2.13137
   r3 = 0.533631;  y3 = 2.27006
   r5 = 0.313594;  y5 = 1.74021
   a = 2.91421;  b = 2

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1
   (r1, r4, r6, x6, y6) = (1 + sqrt(2), -1 + sqrt(2), 1/5, 3/5, 1 + 4*sqrt(2)/5)
   @printf("r1 = %g;  r2 = %g;  r4 = %g\n", r1, r2, r4)
   @printf("r6 = %g;  x6 = %g;  y6 = %g\n", r6, x6, y6)
   b = 2
   (r3, y3, r5, y5, a) = (-66932*sqrt(236*sqrt(2) + 334) - 94628*sqrt(118*sqrt(2) + 167) - 7 + 9*sqrt(2) + 66912*sqrt(2)*sqrt(118*sqrt(2) + 167) + 47328*sqrt(2)*sqrt(236*sqrt(2) + 334), -4 + 2*sqrt(118*sqrt(2) + 167)/(2*sqrt(2) + 3), (-sqrt(2)/2 - 1/2)*(-3 - 2*sqrt(2) + sqrt(20*sqrt(2) + 29)) + 1 + sqrt(2), -3 - 2*sqrt(2) + sqrt(20*sqrt(2) + 29), sqrt(2) + 3/2)
   @printf("r3 = %g;  y3 = %g\n", r3, y3)
   @printf("r5 = %g;  y5 = %g\n", r5, y5)
   @printf("a = %g;  b = %g\n", a, b)
   plot([r1, r1, -r1, -r1, r1], [0, a, b, 0, 0], color=:black, lw=0.5)
   circlef(0, 0, r1, :dodgerblue, beginangle=0, endangle=180)
   circlef(r2, r2, r2, :firebrick1)
   circlef(-r2, r2, r2, :firebrick1)
   circlef(0, r1 - r4, r4, :yellow)
   circlef(x6, y6, r6, :limegreen)
   circlef(-x6, y6, r6, :limegreen)
   circle(r1 - r3, y3, r3, :gray)
   circlef(r5 -r1, y5, r5, :skyblue1)
   if more
       point(r1, 0, " r1", :black, :left, :bottom)
       point(r1, a, "(r1,a)", :black, :right, :bottom)
       point(-r1, b, "(-r1,b)\n", :black, :left, :bottom)
       point(r2, r2, "(r2,r2)", :black)
       point(0, r1 - r4, " r1-r4", :black)
       point(x6, y6, "(x6,y6)", :black)
       point(r1 - r3, y3, "(r1-r3,y3)", :black)
       point(r5 - r1, y5, "(r5-r1,y5)", :black)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       plot!(xlims=(-2.5, 2.7))
   else
       plot!(showaxis=false)
   end
end;

健康保険証を廃止しマイナンバーカードに一本化する政府方針

ChatGPT に答えてもらいました

 

問:　来年(2023)秋に健康保険証を廃止しマイナンバーカードに一本化する政府方針に関してなんらかのコメントはありますか

 

 

ChatGPT

 

私は2021年9月までの情報しか持っていないため、来年秋に健康保険証を廃止しマイナンバーカードに一本化するという政府の方針についての最新の情報は知りません。ただし、一般的な観点からコメントを述べることはできます。

政府が健康保険証をマイナンバーカードに一本化するという方針を採用する場合、いくつかのメリットと懸念点が考えられます。

メリット：

1. 利便性向上：マイナンバーカードは個人識別情報を一元化し、各種サービスへのアクセスが容易になります。健康保険証と他のカードを持ち歩く必要がなくなり、手続きの簡素化が期待できます。
2. 行政効率化：個人情報が一元化されることで、行政手続きの効率が向上する可能性があります。情報の重複を防ぎ、行政サービスのスピードと正確性が向上するでしょう。
3. セキュリティ強化：個人情報を統合することで、情報漏洩や不正利用を防ぐセキュリティ対策が強化されることが期待されます。

懸念点：

1. プライバシー保護：マイナンバーカードは個人の多くの情報を一元化するため、情報漏洩やプライバシーの侵害へのリスクが高まります。適切なセキュリティ対策と情報管理が必要です。
2. トラブルへの対応：システムの障害や誤った情報登録などにより、個人が適切なサービスを受けられないリスクが考えられます。その場合、迅速な対応が求められます。
3. 高齢者やデジタル格差：高齢者やデジタル技術に慣れていない人々は、新しいシステムに適応するのが難しい場合があります。社会的包摂を考慮した対策が必要です。

この方針が実施される場合、政府は十分な情報提供やサポートを行い、市民の利便性向上と個人情報の保護を両立させるよう努める必要があります。また、市民の意見を十分に反映させることも重要です。最新の情報や具体的な方針については、政府の公式発表や報道を注視する必要があります。

算額（その348）

算額（その348）

岐阜県大垣市赤坂町 金生山明星輪寺 元治2年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

長方形内に青の半円 2 個，赤円 1 個，緑円 2 個，浅青円 1 個，白円 2 個，および 薄紅円 2 個，無名円 2 個を入れる。緑円は青，赤，浅青のいずれにも外接し，緑，浅青円は長方形の辺にも接している。

注：算額では，緑円と赤円，浅青円と赤円が接していないような，誤解を与えかねない図になっている。

長方形の左下隅を原点とし，以下のように記号を定める。
青円の半径と中心座標 r1, (r1, 0)
赤円の半径と中心座標 r2, (2r1 - r2, y2)
緑円の半径と中心座標 r3, (r3, y3)
浅青円の半径と中心座標 r5, (r5, y2)
白円の半径と中心座標 r4, (2r1 - r4, y4)
薄紅円の半径と中心座標 r6, (x6, y6)
無名円の半径と中心座標 r7, (x7, y7)
長方形の幅は 2r1，高さは 2y2 である。

r1 を既知とすれば，残りの未知数は 13 変数なので，13 元連立方程式になるが，solve() では一度には解けないので，薄紅円と無名円を除いてまず解を求め，しかる後に得られた解を既知として，薄紅円と無名円のパラメータを求める。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, y2::positive, r3::positive, y3::positive,
     r4::positive, y4::positive, r5::positive, r6::positive, x6::positive,
     y6::positive, r7::positive, x7::positive, y7::positive;

eq1 = (r1 - r2)^2 + y2^2 - (r1 + r2)^2
eq2 = (r1 - r3)^2 + y3^2 - (r1 + r3)^2
eq3 = (r1 - r4)^2 + y4^2 - (r1 + r4)^2
eq5 = (2r1 - r2 - r3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq6 = (r2 - r4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq9 = (r3 - r5)^2 + (y2 - y3)^2 - (r3 + r5)^2
eq13 = 2r5 + 2r2 - 2r1
res1 = solve([eq1, eq2, eq3, eq5, eq6, eq9, eq13], (r2, y2, r3, y3, r4, y4, r5))
solve([eq1, eq2, eq3, eq5, eq6, eq9, eq13, eq4, eq7, eq8, eq10, eq11, eq12], (r2, y2, r3, y3, r4, y4, r5, r6, x6, y6, r7, x7, y7))

   2-element Vector{NTuple{7, Sym}}:
    (3*r1/4, sqrt(3)*r1, r1/3, 2*sqrt(3)*r1/3, 3*r1*(7 - 4*sqrt(3)), 2*r1*(-3 + 2*sqrt(3)), r1/4)
    (3*r1/4, sqrt(3)*r1, r1/3, 2*sqrt(3)*r1/3, 3*r1*(4*sqrt(3) + 7), 2*r1*(3 + 2*sqrt(3)), r1/4)

最初の組のものが適解である。
浅青円の半径は r1/4 すなわち，青円の半径の 1/4 である。

得られた r2, y2, r3, y3, r4, y4, r5 を既知として，次に，r6, x6, y6, r7, x7, y7 に関する連立方程式を解く。

using SymPy

@syms r1::positive, r2::positive, y2::positive, r3::positive, y3::positive,
     r4::positive, y4::positive, r5::positive, r6::positive, x6::positive,
     y6::positive, r7::positive, x7::positive, y7::positive;

(r2, y2, r3, y3, r4, y4, r5) = (3*r1/4, sqrt(Sym(3))*r1, r1/3, 2*sqrt(Sym(3))*r1/3, 3*r1*(7 - 4*sqrt(Sym(3))), 2*r1*(-3 + 2*sqrt(Sym(3))), r1/4)
eq4 = (x6 - r1)^2 + y6^2 - (r1 + r6)^2
eq7 = (2r1 - r2 - x6)^2 + (y2 - y6)^2 - (r2 + r6)^2
eq8 = (2r1 - r2 - x7)^2 + (y2 - y7)^2 - (r2 + r7)^2
eq10 = (x6 - r3)^2 + (y6 - y3)^2 - (r3 + r6)^2
eq11 = (x7 - r3)^2 + (y7 - y3)^2 - (r3 + r7)^2
eq12 = (x7 - r5)^2 + (y2 - y7)^2 - (r5 + r7)^2
res2 = solve([eq4, eq7, eq8, eq10, eq11, eq12], (r6, x6, y6, r7, x7, y7))

   1-element Vector{NTuple{6, Sym}}:
    (3*r1*(-8 + 5*sqrt(3))/22, r1*(32 - 9*sqrt(3))/22, 3*r1/11 + 5*sqrt(3)*r1/11, 3*r1*(-25 + 16*sqrt(3))/143, r1*(109 - 24*sqrt(3))/143, 6*r1/143 + 122*sqrt(3)*r1/143)

二段階に分けて求めたパラメータで図を描き，解が正しいことを確認する。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1
   (r2, y2, r3, y3, r4, y4, r5) = (3*r1/4, sqrt(3)*r1, r1/3, 2*sqrt(3)*r1/3, 3*r1*(7 - 4*sqrt(3)), 2*r1*(-3 + 2*sqrt(3)), r1/4)
   (r6, x6, y6, r7, x7, y7)= (3*r1*(-8 + 5*sqrt(3))/22, r1*(32 - 9*sqrt(3))/22, 3*r1/11 + 5*sqrt(3)*r1/11, 3*r1*(-25 + 16*sqrt(3))/143, r1*(109 - 24*sqrt(3))/143, 6*r1/143 + 122*sqrt(3)*r1/143)
   plot([0, 2r1, 2r1, 0, 0], [0, 0, 2y2, 2y2, 0], color=:black, lw=0.5)
   circlef(r1, 0, r1, :dodgerblue, beginangle=0, endangle=180)
   circlef(r1, 2y2, r1, :dodgerblue, beginangle=180, endangle=360)
   circlef(2r1 - r2, y2, r2, :red)
   circlef(r3, y3, r3, :green)
   circlef(r3, 2y2 - y3, r3, :green)
   circle(2r1 - r4, y4, r4, :orange)
   circle(2r1 - r4, 2y2 - y4, r4, :orange)
   circlef(r5, y2, r5, :skyblue1)
   circlef(x6, y6, r6, :pink)
   circlef(x6, 2y2 - y6, r6, :pink)
   circlef(x7, y7, r7, :orange)
   circlef(x7, 2y2 - y7, r7, :orange)
   if more
       point(r1, 0, "r1", :black, :center, :bottom)
       point(2r1 - r2, y2, "(2r1-r2,y2)", :black, :center, :bottom)
       point(r3, y3, "(r3,y3)", :black, :center, :bottom)
       point(2r1 - r4, y4, "(2r1-r4,y4)   ", :black, :center, :bottom)
       point(r5, y2, "(r5,y2)", :black, :center, :bottom)
       point(x6, y6, "  (x6,y6)", :black, :left, :bottom)
       point(x7, y7, "  (x7,y7)", :black, :left, :vcenter)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その347）

算額（その347）

岐阜県大垣市赤坂町 金生山明星輪寺 元治2年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

菱形の中に交わる赤円が 2 個，白円が 2 個，青円が 5 個入っている。白円の直径が与えられているとき，青円の直径を求めよ。

菱形の中心を原点とし，縦，横の対角線の長さをそれぞれ 2y, 2x とする。
菱形の一辺の長さを z とする。z = sqrt(x^2 + y^2) である。
青，白，赤の円の直径を r1, r2, r3 とする。
上の青円の中心座標を (0, y1) とする。
右の青円の中心座標は (2r1 + 2r2, 0)
右の白円の中心座標は (r1 + r2, 0)
右の赤円の中心座標は (r1 + 2r2 - r3, 0)
である。

白円の半径を 1 として，連立方程式をたてる。簡単な4元連立方程式であるが，solve() で一度に求めるのは無理である。nlsolve() で数値解を求める過程について後半に記す。

しかし，図の見方を変えることにより，3 つの円の関係は非常に簡単に求めることができる。

1. 右の赤円の左端が接する y 軸に平行な補助線を引けば，三角ADEは正三角形になる。
2. ⊿ABC は直角三角形になる
   a. AB は青円と赤円の中心を結ぶ線分で，長さは r1 + r3
   b. AC は菱形の辺に平行な線分で，長さは 2sqrt(r1\*r3)
   c. BC は x 軸に 60°で交わる線分で，長さは AB/2
3. 赤円の直径は青円の直径と白園の直径の和に等しい。

これらの関係について連立方程式をたてる。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive
eq1 = (r1 + r3)^2/4 + 4r1*r3 - (r1 + r3)^2
eq2 = r1 + r2 - r3
res = solve([eq1, eq2], (r1, r3))

   1-element Vector{Tuple{Sym, Sym}}:
    (r2/2, 3*r2/2)

r1 は r2 の 1/2，r3 はr2 の 3/2 である。

白円の直径が与えられるので，r1, r3 は確定され，図を描くために必要なのは，菱形の縦横の対角線の長さと上にある青円の中心の y 座標だけである。

白円の半径 r2 = 1 とすれば，青円，赤円の半径は r1 = 1/2, r3 = 3/2 である。

using SymPy

@syms r1::positive, r2::positive, r3::positive,
     x::positive, y::positive, y1::positive
(r1, r2, r3) = (1//2, 1, 3//2)
eq1 = (r1 + 2r2 - r3)^2 + y1^2 - (r1 + r3)^2
eq2 = (r1*y/x + 4sqrt(Sym(r1*r3)) + r1*x/y)^2 - x^2 - y^2
eq3 = r3*x - y*(2sqrt(Sym(r1*r3)) + r1*x/y)
solve([eq1, eq2, eq3], (x, y, y1))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (4, 4*sqrt(3)/3, sqrt(3))

using Plots

function draw()
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = (1//2, 1, 3//2)
   (x, y, y1) = (4, 4*sqrt(3)/3, sqrt(3))
   plot([x, 0, -x, 0, x], [0, y, 0, -y, 0], color=:black, lw=0.5)
   circle(0, 0, r1, :blue)
   circle(2r1 + 2r2, 0, r1, :blue)
   circle(-2r1 - 2r2, 0, r1, :blue)
   circle(0, y1, r1, :blue)
   circle(0, -y1, r1, :blue)
   circle(r1 + r2, 0, r2, :green)
   circle(-r1 - r2, 0, r2, :green)
   circle(r1 + 2r2 - r3, 0, r3)
   circle(-r1 - 2r2 + r3, 0, r3)
   plot!(showaxis=false)
end;
draw()

---

円の大きさの関係と図を描くために必要なパラメータを一括して求める

solve() では無理なので，nlsolve() で数値街を求める。

using SymPy

@syms r1, r2, r3, x, y, y1
@syms r1::positive, r2::positive, r3::positive, x::positive, y::positive, y1::positive

r2 = 1
r3 = r1 + r2

eq1 = (r1 + 2r2 - r3)^2 + y1^2 - (r1 + r3)^2
eq2 = (r1*y/x + 4sqrt(r1*r3) + r1*x/y)^2 - x^2 - y^2
eq3 = r3*x - y*(2sqrt(r1*r3) + r1*x/y)
eq4 = r1*sqrt(x^2 + y^2) - y*(x - 2r1 - 2r2);

# solve([eq1, eq2, eq3, eq4], (r1, x, y, y1))

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

   y1^2 - (2*r1 + 1)^2 + 1,  # eq1
   -x^2 - y^2 + (4*sqrt(r1)*sqrt(r1 + 1) + r1*x/y + r1*y/x)^2,  # eq2
   x*(r1 + 1) - y*(2*sqrt(r1)*sqrt(r1 + 1) + r1*x/y),  # eq3
   r1*sqrt(x^2 + y^2) - y*(-2*r1 + x - 2),  # eq4

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, x, y, y1) = u
   return [
       y1^2 - (2*r1 + 1)^2 + 1,  # eq1
       -x^2 - y^2 + (4*sqrt(r1)*sqrt(r1 + 1) + r1*x/y + r1*y/x)^2,  # eq2
       x*(r1 + 1) - y*(2*sqrt(r1)*sqrt(r1 + 1) + r1*x/y),  # eq3
       r1*sqrt(x^2 + y^2) - y*(-2*r1 + x - 2),  # eq4
   ]
end;

iniv = [big"0.45", 4.1, 2.2, 1.6]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[0.5000000000000000000002647732231042221299746017703653458755774117300395986631534, 4.00000000000000000000150018442196727793019547881937562510553208700008372428873, 2.309401076758503058036288681801685476240636088014479770497112246296876844518601, 1.732050807568877293528878493906762400441924270519252071343713219148935154969712], true)

r1 = 0.5;  x = 4;  y = 2.3094;  y1 = 1.73205
x/y = 1.73205;  x/r2 = 4
白円直径 = 2;  赤円直径 = 3;  青円直径 = 1
青円の直径は白円の直径の 1/2 である。
なお，菱形は横の対角線の長さが縦の対角線の長さの √3 倍なので，正三角形を二個くっつけたものである。

using Plots

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1
   (r1, x, y, y1) = [0.5, 4, 2.30940107675850305803, 1.73205080756887729352]  # res[1]
   r3 = r1 + r2
   @printf("r1 = %g;  x = %g;  y = %g;  y1 = %g\n", r1, x, y, y1)
   @printf("x/y = %g;  x/r2 = %g\n", x/y, x/r2)
   @printf("白円直径 = %g;  赤円直径 = %g;  青円直径 = %g\n", 2r2, 2r3, 2r1)
   plot([x, 0, -x, 0, x], [0, y, 0, -y, 0], color=:black, lw=0.5)
   circle(0, 0, r1, :blue)
   circle(2r1 + 2r2, 0, r1, :blue)
   circle(-2r1 - 2r2, 0, r1, :blue)
   circle(0, y1, r1, :blue)
   circle(0, -y1, r1, :blue)
   circle(r1 + r2, 0, r2, :green)
   circle(-r1 - r2, 0, r2, :green)
   circle(r1 + 2r2 - r3, 0, r3)
   circle(-r1 - 2r2 + r3, 0, r3)
   if more
       vline!([-r1], color=:orange, lw=0.5)
       segment(0, y1, r1 + 2r2 - r3, 0, :orange)
       (ox, oy) = (r1 + 2r2 - r3 + (r3 - r1)*cos(pi/3), (r3 - r1)*sin(pi/3))
       segment(r1 + 2r2 - r3, 0, r1 + 2r2 - r3 + r3*cos(pi/3), r3*sin(pi/3), :orange)
       segment(2r1 + 2r2, 0, 0, y1, :orange)
       segment(2r1 + 2r2, 0, 0, -y1, :orange)
       segment(0, y1, r1*cos(pi/3), y1 + r1*sin(pi/3), :orange)
       segment(2r1 + 2r2, 0, 2r1 + 2r2 + r1*cos(pi/3), r1*sin(pi/3), :orange)
       point(ox, oy, "(ox,oy)", :orange)
       point(ox, oy, "C\n", :orange, :right, :bottom, mark=false)
       point(r1, 0, " r1", :blue, :left, :bottom)
       point(r1 + 2r2 - r3, 0, "  r1+2r2-r3", :red, :center, :top)
       point(r1 + 2r2 - r3, 0, "B\n", :orange, :center, :bottom, mark=false)
       point(r1 + r2, 0, "  r1+r2", :green, :center, :bottom)
       point(2r1 + 2r2, 0, " 2r1+2r2", :blue, :center, :bottom)
       point(2r1 + 2r2, 0, "  D", :orange, :center, :top, mark=false)
       point(x, 0, " x", :black, :left, :bottom)
       point(0, y1, " y1", :blue)
       point(0, y1, "A ", :orange, :right, mark=false)
       point(0, -y1, "E ", :orange, :right)
       point(0, y, " y", :black, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その346）

算額（その346）

岐阜県大垣市赤坂町 金生山明星輪寺 元治2年(1865)http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

扇の中に赤円，緑円，白円を入れる。扇の半径を与えて，白円が最大となるとき，赤円と白円との直径の和を求めよ。

扇の要を原点とし，扇の半径を R
赤円の直径を r1, (0, R - r1)
白円の直径を r2, (x2, y2)あ
緑円の直径を r3, (x3, y3)
とする。

include("julia-source.txt");

using SymPy

R = 1 とおき，赤円，白円，緑円の半径を決める。

緑円の半径は R - 2r1 + 2r2 + 2r3 = R なので，r3 = r1 - r2 である。

@syms R::positive, x::positive, r1::positive, r2::positive, x2::positive, y2::positive,
     r3::positive, x3::positive, y3::positive

solve(R - 2r1 + 2r2 + 2r3 - R, r3);
r3 = r1 - r2

白円と緑円の中心は原点を通る直線上にある。この直線と扇の右端の直線のなす角を θ とする。
原点から白円と緑円中心までの距離はそれぞれ 1 - 2r1 + r2, 1 - 2r1 + 2r2 + r3 である。
よって
r2 = (1 - 2r1 + r2)sin(θ)
r3 = (1 - 2r1 + 2r2 + r3)sin(θ) である（つまり，r2, r3 を一辺とする直角三角形は相似である）。
これより，r2/r3 = (1 - 2r1 + r2)/(1 - r3) である。

eq0 = r2/r3 - (1 - 2r1 + r2)/(1 - r3) |> expand
eq0 |> println

   2*r1/(-r1 + r2 + 1) - r2/(-r1 + r2 + 1) + r2/(r1 - r2) - 1/(-r1 + r2 + 1)

この式を r2 について解くと，2つの解が得られるが 2 番目のもの r1 + sqrt(1 -2r1)/2 -1/2 が適解である。

solve(eq0, r2)
r2 = r1 + sqrt(1 -2r1)/2 - 1//2

これを r1 の関数とみて，極大値・最大値を求めるには，まず r1 について微分し，微分係数 = 0 となる r1 を求める。

plot(r2, xlims=(0, 0.5), ylims=(0, 0.15), xlabel="r1", ylabel="r2", aspectratio=:none)

diff(r2)
d = 1 - 1/2sqrt(1 - 2r1)
plot(d, xlims=(0, 0.5), ylims=(-0.1, 0.5), xlabel="r1", ylabel="r2'", aspectratio=:none)
hline!([0])
plot!()

d = 1 - 1/2sqrt(1 - 2r1) = 0 を解いて，r1 = 3/8 のときに r2 が最大値 1/8 をとる。

r1 = 3/8 のとき極大かつ最大値をとる。

solve(eq0(r1 => 3//8))[1] |> println  # r2

   1/8

赤円の直径と白円の直径の和は 2\*3/8 + 2\*1/8 = 1。つまり，白円が最大となるとき，赤円と白円の直径の和は扇の半径に等しい。

また，r3 = r1 - r2 なので，r3 = 2/8 である。

次に，各円の中心座標と扇の右端の点の座標を求める。

using SymPy

@syms R::positive, x::positive, x2::positive, y2::positive,
     r3::positive, x3::positive, y3::positive

R = 1
(r1, r2) = (3//8, 1//8)
r3 = r1 - r2
eq1 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
eq2 = x3^2 + (R - r1 - y3)^2 - (r1 + r3)^2
eq4 = x2^2 + y2^2 - (r2 + R - 2r1)^2
eq5 = x3^2 + y3^2 - (R - r3)^2
eq6 = distance(0, 0, x, sqrt(R^2 - x^2), x2, y2) - r2^2;
eq7 = distance(0, 0, x, sqrt(R^2 - x^2), x3, y3) - r3^2;

solve([eq1, eq2, eq4, eq5, eq6], (x, x2, y2, x3, y3))

   2-element Vector{NTuple{5, Sym}}:
    (-1/5 + 8*sqrt(2)/15, 3/10, 9/40, 3/5, 9/20)
    (1/5 + 8*sqrt(2)/15, 3/10, 9/40, 3/5, 9/20)

2 番目のものが適解である。

扇の中心角 = 145.203°
赤円半径 = 0.375;  白円半径 = 0.125;  緑円半径 = 0.25
x = 0.954247;  x2 = 0.3;  x3 = 0.225;  x3 = 0.6;  y4 = 0.45
赤円直径 + 白円直径 = 1 すなわち，赤円の直径と白円の直径の和は扇の半径に等しい

以上の結果に基づき図を描く。

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 1
   (r1, r2) = (3//8, 1//8)
   r3 = r1 - r2
   (x, x2, y2, x3, y3) = (1/5 + 8*sqrt(2)/15, 3/10, 9/40, 3/5, 9/20)
   y = sqrt(R^2 - x^2)
   r3 = r1 - r2
   θ = atand(y, x)
   @printf("扇の中心角 = %g°\n", 180 - 2θ)
   @printf("赤円半径 = %g;  白円半径 = %g;  緑円半径 = %g\n", r1, r2, r3)
   @printf("x = %g;  x2 = %g;  x3 = %g;  x3 = %g;  y4 = %g\n", x, x2, y2, x3, y3)
   @printf("赤円直径 + 白円直径 = %g\n", 2r1 + 2r2)
   plot()
   circle(0, 0, R, :black, beginangle=θ, endangle=180 - θ)
   circle(0, 0, R - 2r1, :black, beginangle=θ, endangle=180 - θ)
   segment(0, 0, x, y, :black)
   segment(0, 0, -x, y, :black)
   circle(0, R - r1, r1, :red)
   circle(x2, y2, r2, :blue)
   circle(-x2, y2, r2, :blue)
   circle(x3, y3, r3, :green)
   circle(-x3, y3, r3, :green)
   if more
       point(0, R, " R", :black)
       point(x, y, "(x,y)", :black)
       point(0, R - r1, " 赤円:r1\n R-r1", :red, :left, :vcenter)
       point(0, R - 2r1, " R-2r1", :red, :center, :bottom)
       point(x3, y3, " 緑円:r3,(x3,y3)", :green, :center)
       point(x2, y2, " 白円:r2,(x2, y2)", :blue)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その345）

算額（その345）

美濃・飛騨の国の和算の歴史 -算額の問題に挑戦しよう-
http://hamaguri.sakura.ne.jp/minohidawasan.htm

岐阜県郡上市八幡町 郡上八幡神社（小野天満宮） 嘉永3年(1850)

鉤股弦に内接する円がある。鉤 + 股 = 17，(鉤^3 + 股^3 + 弦^3)/(弦 - 円の直径) = 213 のとき，弦の長さはいかほどか。

一般解を求めるために，与えられる 2 つの定数を a, b とおく。
鉤 + 股 = a
(鉤^3 + 股^3 + 円の直径^3)/(弦 - 円の直径) = b
eq3 は鉤股弦に内接する円の直径，eq4 は鉤股弦に関するピタゴラスの定理

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, 鉤::positive, 股::positive, 弦::positive, 円の直径::positive

(a, b) = (17, 213)
eq1 = 鉤 + 股 - a
eq2 = (鉤^3 + 股^3 + 円の直径^3)/(弦 - 円の直径) - b
eq3 = 鉤 + 股 - 弦 - 円の直径
eq4 = 鉤^2 + 股^2 - 弦^2
solve([eq1, eq2, eq3, eq4], (鉤, 股, 弦, 円の直径))

   4-element Vector{NTuple{4, Sym}}:
    (a/2 - sqrt(13*a^2 - 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, a/2 + sqrt(13*a^2 - 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, 2*a - sqrt(3*a^2 - 2*b), -a + sqrt(3*a^2 - 2*b))
    (a/2 + sqrt(13*a^2 - 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, a/2 - sqrt(13*a^2 - 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, 2*a - sqrt(3*a^2 - 2*b), -a + sqrt(3*a^2 - 2*b))
    (a/2 - sqrt(13*a^2 + 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, a/2 + sqrt(13*a^2 + 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, 2*a + sqrt(3*a^2 - 2*b), -a - sqrt(3*a^2 - 2*b))
    (a/2 + sqrt(13*a^2 + 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, a/2 - sqrt(13*a^2 + 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2, 2*a + sqrt(3*a^2 - 2*b), -a - sqrt(3*a^2 - 2*b))

鉤 < 股 なので，最初の組のものが適解。
弦は 2*a - sqrt(3*a^2 - 2*b)
a = 17, b = 213 のとき，2*a - sqrt(3*a^2 - 2*b) = 13

(a, b) = (17, 213)
a/2 - sqrt(13*a^2 - 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2,
a/2 + sqrt(13*a^2 - 8*a*sqrt(3*a^2 - 2*b) - 4*b)/2,
2*a - sqrt(3*a^2 - 2*b),
-a + sqrt(3*a^2 - 2*b)

   (5.0, 12.0, 13.0, 4.0)

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鉤, 股, 弦, 円の直径) = (5.0, 12.0, 13.0, 4.0)
   @printf("鉤 = %g, 股 = %g, 弦 = %g, 円の直径 = %g\n", 鉤, 股, 弦, 円の直径)
   plot([0, 股, 0, 0], [0, 0, 鉤, 0], color=:black, lw=0.5)
   circle(円の直径/2, 円の直径/2, 円の直径/2, :blue)
   if more
       point(円の直径/2, 円の直径/2, "", :blue)
       point(股, 0, " 股", :black, :left, :bottom)
       point(0, 鉤, " 鉤", :black, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その344）

算額（その344）

美濃・飛騨の国の和算の歴史 -算額の問題に挑戦しよう-
http://hamaguri.sakura.ne.jp/minohidawasan.htm
米山忠興:等円術(II) : 三斜等円術(一般公式)
https://toyo.repo.nii.ac.jp/?action=repository_uri&item_id=2534&file_id=18&file_no=1

岐阜県郡上市八幡町 郡上八幡神社（小野天満宮） 嘉永3年(1850)

三角形内に界斜を隔てて二等円を入れる。中斜は 257寸、小斜は 68寸、界斜は 40寸である。
大斜はいかほどか。

等円の半径を r, 大斜を x とする他，図のように記号を定め，以下の連立方程式を解く。

solve() で代数解を求めるのは困難なので，nlsolve() で数値解を求める。

include("julia-source.txt");

using SymPy

@syms r, x0, y0, x, x1, x2, x3, dummy

(中斜, 小斜, 界斜) = (257, 68, 40)
eq1 = (小斜 + 界斜 + x3)r +(中斜 + 界斜 +(x - x3))r - x*y0
eq2 = x0^2 + y0^2 - 小斜^2
eq3 = (x - x0)^2 + y0^2 - 中斜^2
eq4 = distance(0, 0, x0, y0, x1, r) - r^2
eq5 = distance(x0, y0, x3, 0, x1, r) - r^2
eq6 = distance(x0, y0, x3, 0, x2, r) - r^2
eq7 = distance(x0, y0, x, 0, x2, r) - r^2;
eq4 = apart(eq4, dummy) |> numerator
eq5 = apart(eq5, dummy) |> numerator
eq6 = apart(eq6, dummy) |> numerator
eq7 = apart(eq7, dummy) |> numerator;

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

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

   r*(x3 + 108) + r*(x - x3 + 297) - x*y0,  # eq1
   x0^2 + y0^2 - 4624,  # eq2
   y0^2 + (x - x0)^2 - 66049,  # eq3
   -r^2*y0^2 - 2*r*x0*x1*y0 + x1^2*y0^2,  # eq4
   -r^2*y0^2 - 2*r*x0*x1*y0 + 2*r*x0*x3*y0 + 2*r*x1*x3*y0 - 2*r*x3^2*y0 + x1^2*y0^2 - 2*x1*x3*y0^2 + x3^2*y0^2,  # eq5
   -r^2*y0^2 - 2*r*x0*x2*y0 + 2*r*x0*x3*y0 + 2*r*x2*x3*y0 - 2*r*x3^2*y0 + x2^2*y0^2 - 2*x2*x3*y0^2 + x3^2*y0^2,  # eq6
   -r^2*y0^2 - 2*r*x^2*y0 + 2*r*x*x0*y0 + 2*r*x*x2*y0 - 2*r*x0*x2*y0 + x^2*y0^2 - 2*x*x2*y0^2 + x2^2*y0^2,  # eq7

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)
   (r, x0, y0, x, x1, x2, x3) = u
   return [
       r*(x3 + 108) + r*(x - x3 + 297) - x*y0,  # eq1
       x0^2 + y0^2 - 4624,  # eq2
       y0^2 + (x - x0)^2 - 66049,  # eq3
       -r^2*y0^2 - 2*r*x0*x1*y0 + x1^2*y0^2,  # eq4
       -r^2*y0^2 - 2*r*x0*x1*y0 + 2*r*x0*x3*y0 + 2*r*x1*x3*y0 - 2*r*x3^2*y0 + x1^2*y0^2 - 2*x1*x3*y0^2 + x3^2*y0^2,  # eq5
       -r^2*y0^2 - 2*r*x0*x2*y0 + 2*r*x0*x3*y0 + 2*r*x2*x3*y0 - 2*r*x3^2*y0 + x2^2*y0^2 - 2*x2*x3*y0^2 + x3^2*y0^2,  # eq6
       -r^2*y0^2 - 2*r*x^2*y0 + 2*r*x*x0*y0 + 2*r*x*x2*y0 - 2*r*x0*x2*y0 + x^2*y0^2 - 2*x*x2*y0^2 + x2^2*y0^2,  # eq7
   ]
end;

iniv = [big"34.98", 67.27, 117.06, 291.98, 60.54, 119.75, 141.28]
res = nls(H, ini=iniv);
println(res);

 (BigFloat[14.00000000000000000002832495142943662706356368549023146540001887379607263571971, 59.99999999999999999997271770164078264394833470869088845531230198310065464574821, 32.0000000000000000000659299852261495342990267083663162723582624564073423302806, 314.9999999999999999999660433365560846765482392433240243631619025246373680740372, 55.99999999999999999997472857143202733080502439021613717408206542155355708244299, 90.99999999999999999998943873473162951157245262928774364469574387503292726000683, 83.99999999999999999998599488059260011385774611340144279061843511421554191942925], true)

大斜は 315 寸である。
図を見るとわかるが，この三角形は随分と扁平で，算額の図とはイメージが異なる（しかし，これが実際の姿である）。

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x0, y0, x, x1, x2, x3) = res[1]
   @printf("r = %g, x0 = %g, y0 = %g, x = %g, x1 = %g, x2 = %g, x3 = %g\n",
       r, x0, y0, x, x1, x2, x3)
   (中斜, 小斜, 界斜) = (257, 68, 40)
   @printf("大斜 = %g;  術による答え = %g\n", x, sqrt((中斜 + 小斜)^2 - 4界斜^2))
   plot([0, x, x0, 0], [0, 0, y0, 0], color=:black, lw=0.5)
   segment(x0, y0, x3, 0, :orange)
   circle(x1, r, r, :blue)
   circle(x2, r, r, :blue)
   if more
       point(x0, y0, "(x0,y0)", :black, :center, :bottom)
       point(x, 0, " x", :black, :left, :bottom)
       point(x3, 0, "x3", :black, :center, :top)
       point(x1, r, "(x1,r) ", :blue, :right, :vcenter)
       point(x2, r, " (x2,r)", :blue, :left, :vcenter)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

 

算額（その343）

算額（その343）

和算で遊ぼう！！ 「三春まちなか寺子屋」2017レポート
https://miharukoma.com/wp-content/uploads/2018/01/%E4%B8%89%E6%98%A5%E3%81%BE%E3%81%A1%E3%81%AA%E3%81%8B%E5%AF%BA%E5%AD%90%E5%B1%8B2017%E3%83%AC%E3%83%9D%E3%83%BC%E3%83%88.pdf

10月 七草木天神社 

長方形内に斜線を描き，甲円 1 個，乙円 2 個，丙円 1 個を入れる。乙円の直径が 9 寸のとき，丙円の直径を求めよ。

長方形の幅を a，斜線の端点の座標を (0, 0), (b, 2r1) とする。
甲円の半径，中心座標を r1, (a - r1, r1)
乙円の半径，中心座標を r2, (2r1 - r2), (x2, r2)
丙円の半径，中心座標を r3, (x3, 2r1 - r3)
とおき，以下の連立方程式を nlsole() で解き，数値解を求める。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r1::positive, r2::positive, x2::positive, r3::positive, x3::positive
@syms a, b, r1, r2, x2, r3, x3
@syms d

eq1 = (a - r1 - x3)^2 + (2r1 - r3 - r1)^2 - (r1 + r3)^2 |> expand |> simplify
eq2 = (a - r1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2 |> expand |> simplify
eq3 = distance(0, 0, b, 2r1, r2, 2r1 - r2) - r2^2
eq4 = distance(0, 0, b, 2r1, x2, r2) - r2^2
eq5 = distance(0, 0, b, 2r1, x3, 2r1 - r3) - r3^2
eq6 = distance(0, 0, b, 2r1, a - r1, r1) - r1^2;
eq3 = apart(eq3, d) |> numerator |> expand |> simplify
eq4 = apart(eq4, d) |> numerator |> expand |> simplify
eq5 = apart(eq5, d) |> numerator |> expand |> simplify
eq6 = apart(eq6, d) |> numerator |> expand |> simplify;

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

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

r2 = 9//2
iniv = [big"40.0", 10, 20, 5, 2, 10]
res = nls(H, ini=iniv);
println(res);
   (BigFloat[42.00000000000000000000000000000000002256826092573487185561562244382886438918157, 10.49999999999999999999999999999999999460987895327460597548276631046112533791355, 18.00000000000000000000000000000000001793984185362760822286679887253225521134492, 5.999999999999999999999999999999999994624492660452202698353623697006093101679257, 1.999999999999999999999998432434804260885043360258691473226485827072352597107708, 12.00000000000000000000000313847108063424084221460499189708023376629061330011589], true)

a = 42; b = 10.5; r1 = 18; x2 = 6;  r3 = 2;  x3 = 12
丙円の直径(2r3) = 4

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 9//2
   (a, b, r1, x2, r3, x3) = res[1]
   @printf("r2 = %g;  a = %g;  b = %g;  r1 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", r2, a, b, r1, x2, r3, x3)
   @printf("丙円の直径(2r3) = %g\n", 2r3)
   plot([0, a, a, 0, 0], [0, 0, 2r1, 2r1, 0], color=:black, lw=0.5)
   segment(0, 0, b, 2r1, :orange)
   circle(a - r1, r1, r1, :green)
   circle(r2, 2r1 - r2, r2, :blue)
   circle(x2, r2, r2, :blue)
   circle(x3, 2r1 - r3, r3)
   if more
       point(a, 0, "a ", :black, :right, :bottom)
       point(b, 2r1, "b ", :black, :right, :bottom)
       point(a, 2r1, "(a,2r1)", :black, :right, :bottom)
       point(a - r1, r1, " 甲円:r1,(a-r1,r1)", :green)
       point(r2, 2r1 - r2, "乙円;r2\n(r2,2r1-r2)", :blue, :center)
       point(x2, r2, "乙円;r2\n(x2,r2)", :blue, :center)
       point(x3, 2r1 - r3, " 丙円:r3,(x3,2r1-r3)", :red)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その342）

算額（その342）

和算で遊ぼう！！ 「三春まちなか寺子屋」2017レポート
https://miharukoma.com/wp-content/uploads/2018/01/%E4%B8%89%E6%98%A5%E3%81%BE%E3%81%A1%E3%81%AA%E3%81%8B%E5%AF%BA%E5%AD%90%E5%B1%8B2017%E3%83%AC%E3%83%9D%E3%83%BC%E3%83%88.pdf
8月 直毘神社、沫なぎ神社

扇面（円弧）内に，図のように相交わる大円の一部（弧）と小円があり，その隙間に内接する全円を入れる。大円の直径が 12寸，小円の直径が 4寸，扇長（扇の半径）が 4寸5分のとき，全円の直径を求めよ。

扇面の円弧と大円の円弧の交点を (x, y) とする。
大円の半径，中心座標を r0, (0, y0)
扇の半径，中心座標を r2, (0, 0)
小円の円の半径，中心座標を r3, (0, r2 - r3)
全円の円の半径，中心座標を r4, (0, r2 - r4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms y0::positive, x::positive, y::positive, r0::positive, r2::positive, r3::positive, r4::positive

(r0, r2, r3) = (6, 9//2, 2)
eq1 = x^2 + (y - y0)^2 - r0^2
eq2 = x^2 + y^2 - r2^2
eq3 = r2 - 2r4 - y0 + r0
eq4 = r0 - 2r2 + 2r4;

res = solve([eq1, eq2, eq3, eq4], (r4, y0, x, y))

   2-element Vector{NTuple{4, Sym}}:
    (-(r0 - 2*r2)/2, 2*r0 - r2, -sqrt(3)*r0*sqrt(-(r0 - 2*r2)*(3*r0 - 2*r2))/(2*(2*r0 - r2)), (3*r0^2 - 4*r0*r2 + 2*r2^2)/(2*(2*r0 - r2)))
    (-(r0 - 2*r2)/2, 2*r0 - r2, sqrt(3)*r0*sqrt(-(r0 - 2*r2)*(3*r0 - 2*r2))/(2*(2*r0 - r2)), (3*r0^2 - 4*r0*r2 + 2*r2^2)/(2*(2*r0 - r2)))

2 番目の対が適解である。また，若干簡約化できる。
これによれば，全円の半径は「扇長 - 大円の半径/2 = 4.5 - 6/2 = 1.5」なので，直径は 3 である。
全円の直径は「2扇長 - 大円の半径 = 2*4.5 - 6 = 3」でもよい。

術では，全円の直径は「(扇長＋大円の直径) - (扇長×大円の直径)/小円の直径」としている。「小円の直径 = 2*扇長*(4/9)」なので
全円の直径は簡約化され，「扇長 - 大円の直径/8 = 4.5 - 12/8 = 3」となる。

   r4 = 1.5;  y0 = 7.5;  x = 3.6;  y = 2.7

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r2, r3) = (6, 9//2, 2)
   r4 = -r0/2 + r2
   y0 = 2*r0 - r2
   x = r0*sqrt(-9*r0^2 + 24*r0*r2 - 12*r2^2)/(2*(2*r0 - r2))
   y = (3*r0^2/2 - 2*r0*r2 + r2^2)/(2*r0 - r2)
   @printf("r4 = %g;  y0 = %g;  x = %g;  y = %g\n", r4, y0, x, y)
   θ = round(Int64, atand(y, x))
   θ2 = round(Int64, atand(y0 - y, x))
   plot()
   circle(0, y0, r0, beginangle=180+θ2, endangle=360-θ2)
   circle(0, 0, r2, :black, beginangle=θ, endangle=180-θ)
   segment(0, 0, x, y)
   segment(0, 0, -x, y)
   circle(0, r2 - r4, r4, :green)
   circle(0, r2 - r3, r3, :blue)
   circle(0, 0, r2 - 2r4, :black, beginangle=θ, endangle=180-θ)
   if more
       point(x, y, "(x,y)")
       point(0, y0, " y0", :red)
       point(0, r2, " r2", :red, :left, :bottom)
       point(0, r0-r2, " r0-r2", :red, :left, :bottom)
       point(0, r2 - r4, " r2-r4", :green)
       point(0, r2 - 2r4, " r2-2r4", :green)
       point(0, r2 - r3, " r2-r3", :blue)
       point(0, r2-2r3, " r2-2r3", :blue)
       #point(0, y0 - r0, " y0-r0", :blue)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;