算額（その2097）

#和算

算額（その2097）

七十二 群馬県富岡市一ノ宮 貫前神社 嘉永2年(1849)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円2個，直角三角形，正五角形

直角三角形の中に正五角形と甲円，乙円を容れる。甲円の直径が 5 寸のとき，乙円の直径はいかほどか。

正五角形を内包する円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (r1 - R*cosd(18), R)
乙円の半径と中心座標を r2(x2, r2 -R*cosd(36))
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive,
     x2::positive;
eq1 = dist2(0, R, R*cosd(Sym(18)), R*sind(Sym(18)), r1 - R*cosd(Sym(18)), R, r1)
eq2 = dist2(0, R, R*cosd(Sym(18)), R*sind(Sym(18)), x2, r2 - R*cosd(Sym(36)), r2)
eq3 = dist2(R*sind(Sym(36)), -R*cosd(Sym(36)), R*cosd(Sym(18)), R*sind(Sym(18)), x2, r2 - R*cosd(Sym(36)), r2);

SymPy の能力的に，一度には解けないので，逐次解いていく。

まず，eq1 を解き，R を求める。

ans_R = solve(eq1, R)[1] |> simplify
ans_R |> println

   r1*((-sqrt(10) + 5*sqrt(2))*sqrt(sqrt(5) + 5) + 8*sqrt(5))/10

eq2 の R に ans_R を代入し，x2 を求める。

eq12 = eq2(R => ans_R) |> simplify
ans_x2 = solve(eq12, x2)[1]  # 1 of 2
@syms d
ans_x2 = apart(ans_x2, d) |> simplify |> sympy.sqrtdenest|> simplify
ans_x2 |> println

   (-4*sqrt(2)*r1 - 2*r1*sqrt(5 - sqrt(5)) - sqrt(10)*r2 + 5*sqrt(2)*r2)/((-3 + sqrt(5))*sqrt(5 - sqrt(5)))

eq3 の R に ans_R，x2 に ans_x2 を代入し，r2 を求める。

eq13 = eq3(R => ans_R, x2 => ans_x2) |> simplify |> numerator
ans_r2 = solve(eq13, r2)[2]  # 2 of 2
ans_r2 = apart(ans_r2/r1, d) |> x -> x*r1
ans_r2 |> println

   r1*(-sqrt(2)*sqrt(sqrt(5) + 5) - sqrt(5) + 3 + sqrt(10)*sqrt(sqrt(5) + 5)/2)

甲円の直径が 5 寸のとき，各パラメータは以下のように計算される。
乙円の直径は 2*3.032399997699489 = 6.064799995398978 である。

r1 = 5/2
t = sqrt(5 + √5)
u = sqrt(5 - √5)
r2 = r1*(3 - √5 + t*(√10/2 - √2))
x2 = (r2*(5√2 - √10) -　r1*(4√2 + 2*u))/((√5 - 3)*u)
R = r1*((5√2 - √10)*t + 8√5)/10
(r2, x2, R)

   (3.032399997699489, 8.347481064940814, 7.1007915155952475)

function draw(r1, more)
    pyplot(size=(600, 600), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   t = sqrt(5 + √5)
   u = sqrt(5 - √5)
   r2 = r1*(3 - √5 + t*(√10/2 - √2))
   x2 = (r2*(5√2 - √10) -　r1*(4√2 + 2*u))/((√5 - 3)*u)
   R = r1*((5√2 - √10)*t + 8√5)/10
   @printf("甲円の直径が %g のとき，乙円の直径は %.15g である。\n", 2r1, 2r2)
   plot([R*cosd(18), 0, -R*cosd(18), -R*sind(36), R*sind(36), R*cosd(18)],
        [R*sind(18), R, R*sind(18), -R*cosd(36), -R*cosd(36), R*sind(18)], color=:magenta, lw=0.5)
   circle(r1 - R*cosd(18), R, r1)
   circle(x2, r2 - R*cosd(36), r2, :blue)
   (x01, y01) = Float64.(intersection(0, R, R*cosd(18), R*sind(18), -R*cosd(18), 0, -R*cosd(18), -R*sind(36)))
   (x02, y02) = Float64.(intersection(0, R, R*cosd(18), R*sind(18), R*sind(36), -R*cosd(36), -R*cosd(18), -R*cosd(36)))
   plot!([x01, x02, x01, x01], [y02, y02, y01, y02], color=:green, lw=0.5)
   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 - R*cosd(18), R, "甲円:r1,(r1-R*cosd(18),R)", :red, :left, :bottom, delta=delta, deltax=-8delta)
       point(x2, r2 - R*cosd(36), "乙円:r2,(x2,r2-R*cosd(36))", :blue, :left, :bottom, delta=delta, deltax=-8delta)
       point(0, R, "(0,R)", :magenta, :center, delta=-2delta)
       point(R*cosd(18), R*sind(18), "(R*cosd(18),Rsind(18)) ", :magenta, :right, :vcenter)
       point(R*sind(36), -R*cosd(36), "(R*sind(36),-R*cosd(36))", :magenta, :right, :bottom, delta=delta)
   end  
end;

draw(5/2, true)

算額（その2096）

算額（その2096）

六十四 群馬県安中市板鼻 鷹巣神社 天保11年(1840)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円10個，外円

外円の中に甲乙丙丁戊の円を容れる。甲円，乙円の直径が 4 寸，9 寸のとき，外円の直径はいかほどか。
異版では，名前の付け方が違うが，本問に対応付けると「丙円の直径も求めよ」としている。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (r1, y1)
乙円の半径と中心座標を r2, (r2, y2)
丙円の半径と中心座標を r3, (0, y3)
丁円の半径と中心座標を r4, (x4, y4)
戊円の半径と中心座標を r5, (x5. y5)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, y1::positive,
     r2::positive, y2::negative, r3::positive, y3::positive,
     r4::positive, x4::positive, y4::positive,
     r5::positive, x5::positive, y5::positive;
eq1 = r1^2 + (y1 - y3)^2 - (r1 + r3)^2
eq2 = (x4 - r1)^2 + (y1 - y4)^2 - (r1 + r4)^2
eq3 = (x5 - r1)^2 + (y1 - y5)^2 - (r1 + r5)^2
eq4 = r2^2 + (y3 - y2)^2 - (r2 + r3)^2
eq5 = (r2 - x4)^2 + (y4 - y2)^2 - (r2 + r4)^2
eq6 = (x5 - r2)^2 + (y5 - y2)^2 - (r2 + r5)^2
eq7 = x4^2 + (y3 - y4)^2 - (r3 + r4)^2
eq8 = (x5 - x4)^2 + (y5 - y4)^2 - (r4 + r5)^2
eq9 = r1^2 + y1^2 - (R - r1)^2
eq10 = r2^2 + y2^2 - (R - r2)^2
eq11 = x5^2 + y5^2 - (R - r5)^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)
   (R, y1, y2, r3, y3, r4, x4, y4, r5, x5, y5) = u
   return [
       r1^2 - (r1 + r3)^2 + (y1 - y3)^2,  # eq1
       (-r1 + x4)^2 - (r1 + r4)^2 + (y1 - y4)^2,  # eq2
       (-r1 + x5)^2 - (r1 + r5)^2 + (y1 - y5)^2,  # eq3
       r2^2 - (r2 + r3)^2 + (-y2 + y3)^2,  # eq4
       -(r2 + r4)^2 + (r2 - x4)^2 + (-y2 + y4)^2,  # eq5
       (-r2 + x5)^2 - (r2 + r5)^2 + (-y2 + y5)^2,  # eq6
       x4^2 - (r3 + r4)^2 + (y3 - y4)^2,  # eq7
       -(r4 + r5)^2 + (-x4 + x5)^2 + (-y4 + y5)^2,  # eq8
       r1^2 + y1^2 - (R - r1)^2,  # eq9
       r2^2 + y2^2 - (R - r2)^2,  # eq10
       x5^2 + y5^2 - (R - r5)^2,  # eq11
   ]
end;

(r1, r2) = (4/2, 9/2)
iniv = BigFloat[9.1, 6.9, -1.1, 1.9, 3.5, 0.94, 2.8, 4.0, 1.8, 5.3, 5.1]

res = nls(H, ini=iniv)

   ([9.142857142857142, 6.857142857142857, -1.1428571428571428, 1.9393939393939394, 3.463203463203463, 0.9411764705882353, 2.823529411764706, 4.033613445378151, 1.7777777777777777, 5.333333333333333, 5.079365079365079], true)

甲円，乙円の直径が 4 寸，9 寸のとき，外円の直径は 128/7 = 18 + 2/7 ≒ 18.2857 である。
丙円は 128/33 = 3 + 29/33 ≒ 3.87879 である。

すべてのパラメータは以下のとおりである。

   R = 9.14286;  y1 = 6.85714;  y2 = -1.14286;  r3 = 1.93939;  y3 = 3.4632;  r4 = 0.941176;  x4 = 2.82353;  y4 = 4.03361;  r5 = 1.77778;  x5 = 5.33333;  y5 = 5.07937

function draw(r1, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, y1, y2, r3, y3, r4, x4, y4, r5, x5, y5) = [37.4/2, 12.9, -3.5, 7.9/2, 5.6, 2.3, 6.2, 5.9, 7.9/2, 12.2, 8.1]
   (R, y1, y2, r3, y3, r4, x4, y4, r5, x5, y5) = res[1]
   @printf("甲円，乙円の直径が %g, %g のとき，外円の直径は %g, 丙円の直径は %g である。\n", 2r1, 2r2, 2R, 2r3)
   @printf("R = %g;  y1 = %g;  y2 = %g;  r3 = %g;  y3 = %g;  r4 = %g;  x4 = %g;  y4 = %g;  r5 = %g;  x5 = %g;  y5 = %g\n", R, y1, y2, r3, y3, r4, x4, y4, r5, x5, y5)
   plot()
   circle(0, 0, R)
   circle2(r1, y1, r1, :blue)
   circle2(r2, y2, r2, :green)
   circle(0, y3, r3, :magenta)
   circle2(x4, y4, r4, :orange)
   circle2(x5, y5, r5, :brown)
   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, :center, :bottom, delta=delta/2)
       point(r1, y1, "甲円:r1\n(r1,y1)", :blue, :center, delta=-delta/2)
       point(r2, y2, "乙円:r2,(r2,y2)", :green, :center, delta=-delta/2)
       point(0, y3, "丙円:r3\n(0,y3)", :magenta, :center, delta=-delta/2)
       point(x4, y4, "丁円:r4,(x4,y4)", :black, :left, :bottom, delta=delta/2, deltax=-2delta)
       point(x5, y5, "戊円:r5\n(x5,y5)", :brown, :center, :bottom, delta=delta/2)
   end  
end;

draw(4/2, 9/2, true)

 

#和算