裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

算額(その1012)

2024年05月29日 | Julia

算額(その1012)

八七 加須市大越六間 天神社 明治14年(1881)

埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

5 個の等円が隣同士外接し,かつ,内側にある正五角形の辺にも外接している。等円の直径が 49.0165 寸のとき,正五角形の一辺の長さを求めよ。

等円の中心は,原点 O を中心とする半径 OA の円周上にある。OA = r1 + r2*cos(36°)
等円の半径 AD = r1
AB/2 = r1 = OA*sin(36°)
正五角形の頂点は,原点 O を中心とする半径 OC の円周上にある。 OC = r2 = 4r1*(-sqrt(5 - √5) + 2√2)/((1 + √5)*sqrt(5 - √5))
∠COD は 36°
正五角形の辺と等円の接点は,原点 O を中心とする半径 OD の円周上にある。
正五角形の一辺の長さは 2CD = 2OC*sin(36°)

等円の直径 2r1 = 49.0165 のとき,五角形の一辺の長さは 24.97515419514368 である。

術では五角形の一辺の長さは 25 寸とあるが,正確ではない。
正確に25 寸であるためには等円の直径が 49.0652626376286 でなければならない。
問で「等円径四十九寸令壱リ六毛五糸」とあるのが誤記で,「等円径四十九寸令分六釐五毛」なのかもしれない。

以下により,r1, r2 を求め,正五角形の一辺の長さ 2r2*sin(36°) = 24.97515419514368 を求める。

include("julia-source.txt");

using SymPy

@syms r1::poitive, r2::poitive

eq1 = (r1 + r2*cosd(Sym(36)))*sind(Sym(36)) - r1;
solve(eq1, r2)[1] |> println

   4*r1*(-sqrt(5 - sqrt(5)) + 2*sqrt(2))/((1 + sqrt(5))*sqrt(5 - sqrt(5)))

r1 = 49.0165/2
r2 = 4*r1*(-sqrt(5 - sqrt(5)) + 2*sqrt(2))/((1 + sqrt(5))*sqrt(5 - sqrt(5)))
2r2*sind(36)

   24.97515419514368

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 49.0165/2
   r2 = 4r1*(-sqrt(5 - √5) + 2√2)/((1 + √5)*sqrt(5 - √5))
   l = 2r2*sind(36)
   println("i = $l")
   θ = 54:72:414
   x = cosd.(θ)
   y = sind.(θ)
   plot()
   r0 = r1 + r2*cosd(36)
   for i = 1:5
       segment(r2*x[i], r2*y[i], r2*x[i + 1], r2*y[i + 1], :blue, lw=0.5)
       circle(r0*cosd(θ[i] + 36), r0*sind(θ[i] + 36), r1)
   end
   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)
       circle(0, 0, r0, :gray70)
       #circle(r0*cosd(θ[1] + 36), r0*sind(θ[1] + 36), r1, :magenta)
       point(r0*cosd(θ[1] + 36), r0*sind(θ[1] + 36), " A", :green, :left, :bottom, delta=delta/2)
       segment(0, 0, r0*cosd(θ[1] + 36), r0*sind(θ[1] + 36), :green)
       point(r0*cosd(θ[5] + 36), r0*sind(θ[5] + 36), " B")
       segment(r0*cosd(θ[1] + 36), r0*sind(θ[1] + 36), r0*cosd(θ[5] + 36), r0*sind(θ[5] + 36), :green)
       point(0, 0, " O", :green, :left, delta=-delta/2)
       circle(0, 0, r2, :gray80)
       point(r2*cosd(θ[1]), r2*sind(θ[1]), "C", :green, :left, :bottom, delta=0.1delta, deltax=0.3delta)
       segment(0, 0, r2*cosd(θ[1]), r2*sind(θ[1]), :green)
       point(0, r2*sind(θ[1]), " D", :green, :left, delta=-delta/2)
       circle(0, 0, r2*sind(θ[1]), :gray80)
   end
end;

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1011)

2024年05月29日 | Julia

算額(その1011)

八六 加須市多聞寺 愛宕神社 明治13年(1880)

埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

長方形の中に甲円 4 個,乙円 2 個,丙円 8 個と菱形が入っている。甲円の直径が 1 寸のとき,乙円の直径はいかほどか。

長方形の長辺と短辺の長さを 2a, 2b; b = 2r1
菱形の短い方の対角線の長さを 2c; 長い方の対角線は 2b
甲円の半径と中心座標を r1, (a - r1, r1)
乙円の半径と中心座標を r2, (0, y2)
丙円の半径と中心座標を r3, (x3, b - r3), (r3, r3)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::poitive, b::poitive, c::poitive,
     r1::poitive, r2::poitive, y2::poitive,
     r3::poitive, x3::poitive
b = 2r1
eq1 = (a - r1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2 |> expand
eq2 = r3^2 + (y2 - r3)^2 - (r2 + r3)^2 |>expand
eq3 = dist2(0, b, c, 0, x3, b - r3, r3)/4r1
eq4 = dist2(0, b, c, 0, a - r1, r1, r1)/(4r1^2)
eq5 = dist2(0, b, c, 0, 0, y2, r2)
eq6 = dist2(0, b, c, 0, r3, r3, r3)/(4c*r1);

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)
   (a, c, r2, y2, r3, x3) = u
   return [
       a^2 - 2*a*r1 - 2*a*x3 + r1^2 - 4*r1*r3 + 2*r1*x3 + x3^2,  # eq1
       -r2^2 - 2*r2*r3 + r3^2 - 2*r3*y2 + y2^2,  # eq2
       -c*r3*x3 - r1*r3^2 + r1*x3^2,  # eq3
       a^2 - a*c - 2*a*r1 + c*r1,  # eq4
       4*c^2*r1^2 - 4*c^2*r1*y2 - c^2*r2^2 + c^2*y2^2 - 4*r1^2*r2^2,  # eq5
       c*r1 - c*r3 - 2*r1*r3 + r3^2,  # eq6
   ]
end;
r1 = 1/2
iniv = BigFloat[1.17, 0.29, 0.17, 0.39, 0.13, 0.17]
res = nls(H, ini=iniv)

   ([1.1666666666666667, 0.2916666666666667, 0.1701388888888889, 0.3923611111111111, 0.125, 0.16666666666666666], true)

甲円の直径が 1 寸のとき,乙円の直径は 0.340278 寸(3 分 4 厘有奇)である。
注1:「答」の乙円径□□四釐有奇」の欠損文字は「三分」であろう。
注2:「術」「甲円の直径を 45 倍し 144 で割る」というのは数値解を求めたということであろうが,45/144=0.3125 は答えとも不一致である。17/50 = 0.34 は良い近似であろう。

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

   a = 1.16667;  b = 1;  c = 0.291667;  r2 = 0.170139;  y2 = 0.392361;  r3 = 0.125;  x3 = 0.166667

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, c, r2, y2, r3, x3) = res[1]
   b = 2r1
   @printf("甲円の直径が %g のとき,乙円の直径は %g である。\n", 2r1, 2r2)
   @printf("a = %g;  b = %g;  c = %g;  r2 = %g;  y2 = %g;  r3 = %g;  x3 = %g\n", a, b, c, r2, y2, r3, x3)
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:black, lw=0.5)
   plot!([0, c, 0, -c, 0], [-b, 0, b, 0, -b], color=:orange, lw=0.5)
   circle4(a - r1, b - r1, r1)
   circle22(0, y2, r2, :green)
   circle4(x3, b - r3, r3, :blue)
   circle4(r3, r3, r3, :blue)
   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(a - r1, r1, "甲円:r1,(a-r1,r1)", :red, :center, :bottom, delta=delta/2)
       point(0, y2, " 乙円:r2,(0,y2)", :green, :left, :vcenter)
       point(r3, b - r3, " 丙円:r3,(x3,b-r3)", :blue, :left, :vcenter)
       point(r3, r3, " 丙円:r3,(r3,r3)", :blue, :left, :vcenter)
       point(a, b, "(a,b)", :black, :right, :bottom, delta=delta/2)
   end
end;

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村