算額（その475）

算額（その475）

宮城県石巻市雄勝町 葉山神社 明治17年(1884)

徳竹亜紀子，谷垣美保: 2021年度の算額調査，仙台高等専門学校名取キャンパス 研究紀要，第 58 号, p.7-28, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2022/03/kiyo2022-2.pdf

外円内に同じ大きさの菱形 2 個と円が入っている。円の直径が 1 寸のとき外円の直径はいかほどか。

外円の半径と中心座標を r0, (0, 0)
菱形の長い方の対角線と短い方の対角線の長さを 2b, 2a
円の半径と中心座標を r1, (x1, y1)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms r0::positive, a::positive, b::positive,
     r1::positive, x1::positive, y1::positive;

eq1 = 2a + 2b - 2r0
eq2 = (r0 - a)^2 + b^2 - r0^2
eq3 = x1^2 + y1^2 - (r0 - r1)^2
eq4 = distance(0, r0 - 2a, b, r0 - a, x1, y1) - r1^2
eq5 = distance(0, r0 - 2a, a, b - r0, x1, y1) - r1^2;

a, b は r0 が分かれば決まる。しかし，r0, x1, y1 は SymPy では有限の時間内に計算できないようなので，a, b も含めて数値解を求める。

solve([eq1, eq2], (a, b))

   1-element Vector{Tuple{Sym, Sym}}:
    (-sqrt(2)*r0/2 + r0, sqrt(2)*r0/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 v, r.f_converged
end;

function H(u)
   (r0, a, b, x1, y1) = u
   return [
       2*a + 2*b - 2*r0,  # eq1
       b^2 - r0^2 + (-a + r0)^2,  # eq2
       x1^2 + y1^2 - (r0 - r1)^2,  # eq3
       -r1^2 + (y1 - (a^2*y1 - 2*a*b^2 + a*b*x1 + b^2*r0)/(a^2 + b^2))^2 + (-b*(2*a^2 - a*r0 + a*y1 + b*x1)/(a^2 + b^2) + x1)^2,  # eq4
       -r1^2 + (y1 - (-a*(2*a^2 - a*r0 - 2*a*x1 + a*y1 - b*x1 + 2*r0*x1) + y1*(5*a^2 + 4*a*b - 8*a*r0 + b^2 - 4*b*r0 + 4*r0^2))/(5*a^2 + 4*a*b - 8*a*r0 + b^2 - 4*b*r0 + 4*r0^2))^2 + (-a*(4*a^2 + 2*a*b - 6*a*r0 + a*x1 + 2*a*y1 - b*r0 + b*y1 + 2*r0^2 - 2*r0*y1)/(5*a^2 + 4*a*b - 8*a*r0 + b^2 - 4*b*r0 + 4*r0^2) + x1)^2,  # eq5
   ]
end;

r1 = 1//2
iniv = [big"85", 25, 60, 48, 34] ./ 37
iniv = [big"2.297", 0.676, 1.622, 1.297, 0.919]
res = nls(H, ini=iniv);

   外円の直径 = 2.38017;  r0 = 1.19008;    a = 0.348568;  b = 0.841517;  x1= 0.653281;  y1= 0.222351

外円の直径は 2寸3分8毛有奇である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1//2
   (r0, a, b, x1, y1) = res[1]
   @printf("外円の直径 = %g;  r0 = %g;    a = %g;  b = %g;  x1= %g;  y1= %g\n", 2r0, r0, a, b, x1, y1)
   plot([0, b, 0, -b, 0], [r0 - 2a, r0 - a, r0, r0 - a, r0 - 2a], color=:blue, lw=0.5)
   plot!([0, a, 0, -a, 0], [-r0, b - r0, r0 - 2a, b - r0, -r0], color=:blue, lw=0.5)
   circle(0, 0, r0, :magenta)
   circle(x1, y1, r1)
   circle(-x1, y1, r1)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, r0, " r0", :magenta, :left, :bottom, delta=delta/2)
       point(0, -r0, " -r0", :magenta, :left, :top)
       point(0, r0 - a, " r0-a", :blue, :left, :vcenter)
       point(b, r0 - a, "(b,r0-a)  ", :blue, :right, :vcenter)
       point(0, r0 - 2a, " r0-2a", :blue, :left, :vcenter)
       point(0, b - r0, " b-r0", :blue, :left, :vcenter)
       point(a, b - r0, " (a,b-r0)", :blue, :left, :vcenter)
       point(x1, y1, "等円:r1,(x1,y1)", :red, :center, :top, delta=-delta)
   end
end;

 

算額（その474）

算額（その474）

宮城県石巻市雄勝町 葉山神社 明治17年(1884)

徳竹亜紀子，谷垣美保: 2021年度の算額調査，仙台高等専門学校名取キャンパス 研究紀要，第 58 号, p.7-28, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2022/03/kiyo2022-2.pdf

正方形の中に三角形（圭），大円，天円，小円が入っている。小円の直径が 1 寸のとき，大円の直径を求めよ。

正方形の一辺の長さを 2a
大円の半径と中心座標を r1, (r1, 0), (0, r1)
小円の半径と中心座標を r2, (0, 0)
天円の半径と中心座標を r3, (a - r3, y3)
とおき，以下の連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms a::positive, r1::positive, r2::positive,
     r3::positive, y3::positive;

eq1 = (r1 - r3)^2 + y3^2 - (r1 + r3)^2
eq2 = distance(0, 2r1 - a, a, a, a - r1, 0) - r1^2
eq3 = (a - 2r1)*a/sqrt((2a - 2r1)^2 + a^2) - r2
eq4 = (a - y3)*r1 - r3*a;

どういうわけか，一度に全部の解を求めようとすると y3 が求まらない（すべての解の式に y3 が含まれ，y3 は未知数のまま）。
そこで，まず a, r1 を求め，次にそれを既知として r3, y3 を求める。
問では r2 が与えられているが，算額の常套手段で，r2 は未知数のまま解の式に残す。

res = solve([eq2, eq3], (a, r1))

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

求まった解 a, r1 を既知として r3, y3 を求める。

@syms a::positive, r1::positive, r2::positive,
     r3::positive, y3::positive;

(a, r1) = (5*r2, 5*r2/3)
eq1 = (r1 - r3)^2 + y3^2 - (r1 + r3)^2
eq4 = (a - y3)*r1 - r3*a
solve([eq1, eq4], (r3, y3))

   1-element Vector{Tuple{Sym, Sym}}:
    (-10*r2*(-1/9 + sqrt(10)/9)/3 + 5*r2/3, 10*r2*(-1/9 + sqrt(10)/9))

天円の半径は r2 * 5(11 - 2*sqrt(10))/27 となり，小円の半径の 5(11 - 2*sqrt(10))/27 = 0.865823088826526 倍である。

-10*r2*(-1//9 + sqrt(Sym(10))//9)//3 + 5*r2//3 |> simplify |> println

   5*r2*(11 - 2*sqrt(10))/27

5(11 - 2*sqrt(10))/27

   0.865823088826526

小円の直径が 1 寸のとき，天円の直径は 0.865823088826526 = 8分6厘5毛8糸 有奇（あまりあり）。

   天円の直径 = 0.865823;  a = 2.5;  r1 = 0.833333;  r3= 0.432912;  y3= 1.20127

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1//2
   (a, r1) = (5*r2, 5*r2/3)
   (r3, y3) = (-10*r2*(-1/9 + sqrt(10)/9)/3 + 5*r2/3, 10*r2*(-1/9 + sqrt(10)/9))
   @printf("天円の直径 = %g;  a = %g;  r1 = %g;  r3= %g;  y3= %g\n", 2r3, a, r1, r3, y3)
   plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:black, lw=0.5)
   circle42(0, a - r1, r1)
   circle4(a - r3, y3, r3, :blue)
   circle(0, 0, r2, :magenta)
   plot!([-a, 0, a], [a, 2r1 - a, a], color=:green, lw=0.5)
   plot!([-a, 0, a], [-a, a - 2r1, -a], color=:green, lw=0.5)
   x = r1*a/(a - 3)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(a - 2r1, 0, " a-2r1")
       point(0, a - 2r1, " a-2r1")
       point(0, a - r1, " 大円:r1\n (0,a-r1)", :red, :left, :vcenter)
       point(a - r1, 0, "r1", :red, :center, :bottom, delta=delta)
       point(a, 0, "a ", :red, :right, :bottom, delta=delta)
       point(0, 0, " 小円\n r2", :magenta, :left, :top, delta=-delta/2)
       point(0, r2, " r2", :magenta, :left, :top, delta=-delta)
       point(a - r3, y3, "天円:r3 \n(a-r3,y3) ", :blue, :right, :vcenter)
       
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;