裏 RjpWiki

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

算額(その392)

2023年08月16日 | Julia

算額(その392)

埼玉県鴻巣市 新井稲荷神社 明治25年(1892)
https://yamabukiwasan.sakura.ne.jp/ymbk351.pdf

外円の中に弦を入れ,その上下に大円,小円を入れる。弦の長さが 4 寸,小円の直径が 1 寸のとき,大円の直径はいかほどか。

外円の半径と中心座標を r0, (0, 0)
大円の半径と中心座標を r1, (x1, r1 - r0 + 2r2)
小円の半径と中心座標を r2, (x2, 3r2 - r0)
弦の長さを「弦」
とおき,以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, x1::negative, r2::positive, x2::positive, 弦::positive;

eq1 = (x2 - x1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x2^2 + (3r2 - r0)^2 - (r0 - r2)^2
eq3 = x1^2 + (r1 -r0 +2r2)^2 - (r0 - r1)^2
eq4 = r0^2 - (2r2 - r0)^2 - (弦//2)^2

res = solve([eq1, eq2, eq3, eq4], (r0, r1, x1, x2))

   2-element Vector{NTuple{4, Sym}}:
    (r2 + 弦^2/(16*r2), (-16*r2^2 + 2*弦^2 - sqrt(-4*r2 + 弦)*sqrt(4*r2 + 弦)*sqrt(-16*r2^2 + 弦^2) - sqrt(-256*r2^4 + 弦^4))/(32*r2), sqrt(-16*r2^2 + 弦^2)/4 + sqrt(16*r2^2 + 弦^2)/4, sqrt(-16*r2^2 + 弦^2)/2)
    (r2 + 弦^2/(16*r2), (-16*r2^2 + 2*弦^2 - sqrt(-4*r2 + 弦)*sqrt(4*r2 + 弦)*sqrt(-16*r2^2 + 弦^2) + sqrt(-256*r2^4 + 弦^4))/(32*r2), sqrt(-16*r2^2 + 弦^2)/4 - sqrt(16*r2^2 + 弦^2)/4, sqrt(-16*r2^2 + 弦^2)/2)

2 組目の解が適解である。

大円の直径は 3.93649167310371 寸である。

2*res[2][2](r2 => 1//2, 弦 => 4).evalf() |> println

    3.93649167310371

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, 弦) = (1//2, 4)
   (r0, r1, x1, x2) = (r2 + 弦^2/(16*r2), (-16*r2^2 + 2*弦^2 - sqrt(-4*r2 + 弦)*sqrt(4*r2 + 弦)*sqrt(-16*r2^2 + 弦^2) + sqrt(-256*r2^4 + 弦^4))/(32*r2), sqrt(-16*r2^2 + 弦^2)/4 - sqrt(16*r2^2 + 弦^2)/4, sqrt(-16*r2^2 + 弦^2)/2)
   b = sqrt(r0^2 - (2r2 - r0)^2)
   @printf("r0 = %g;  r1 = %g;  x1 = %g;  x2= %g\n大円の直径 = %g;  弦 = %g;  小円の直径 = %g\n", r0, r1, x1, x2, 2r1, 2b, 2r2)
   plot()
   circle(0, 0, r0, :black)
   circle(x1, r1 - r0 +2r2, r1)
   circle(x2, 3r2 - r0, r2, :blue)
   circle(0, r2 - r0, r2, :blue)
   segment(-b, 2r2 - r0, b, 2r2 - r0)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(x1, r1 - r0 + 2r2, " 大円:r1,(x1,r1-r0+2r2) ", :red, :right, :vcenter)
       point(x2, 3r2 - r0, " 小円:r1,(x1,3r2−r0) ", :blue, :right, :vcenter)
       point(0, r2 - r0, " r2-r0", :blue)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

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

算額(その391)

2023年08月16日 | Julia

算額(その391)

埼玉県鴻巣市 新井稲荷神社 明治25年(1892)https://yamabukiwasan.sakura.ne.jp/ymbk351.pdf

大円,中円,小円,乙円,丙円が図のように配置されている。小円,乙円,丙円の直径がそれぞれ 5 寸,4 寸,3 寸のとき,大円の直径はいかほどか。

大円の半径と中心座標を r1, (0, r3 + 2r5 + 2r4 + r1)
中円の半径と中心座標を r2, (x2, y2)
小円の半径と中心座標を r3, (0, 0)
乙円の半径と中心座標を r4, (0, r3 + 2r5 + r4) 
丙円の半径と中心座標を r5, (0, r3 + r5)
とおき,以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x2::positive, y2::positive,
     r3::positive, r4::positive, r5::positive;

eq1 = x2^2 + (r3 + 2r5 + 2r4 + r1 - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + y2^2 - (r2 + r3)^2
eq3 = x2^2 + (r3 + 2r5 + r4 - y2)^2 - (r2 + r4)^2
eq4 = x2^2 + (r3 + r5 - y2)^2 - (r2 + r5)^2

res = solve([eq1, eq2, eq3, eq4], (r1, r2, x2, y2))

   2-element Vector{NTuple{4, Sym}}:
    (-r3*r4^2/(r3*r4 - r3*r5 - r5^2), r5*(r3 + r5)*(r4 + r5)/(r3*r4 - r5^2), -2*sqrt(r3)*sqrt(r4)*r5*sqrt(r3 + r5)*sqrt(r4 + r5)/(r3*r4 - r5^2), (r3^2*r4 + r3*r4*r5 - r4*r5^2 - r5^3)/(r3*r4 - r5^2))
    (-r3*r4^2/(r3*r4 - r3*r5 - r5^2), r5*(r3 + r5)*(r4 + r5)/(r3*r4 - r5^2), 2*sqrt(r3)*sqrt(r4)*r5*sqrt(r3 + r5)*sqrt(r4 + r5)/(r3*r4 - r5^2), (r3^2*r4 + r3*r4*r5 - r4*r5^2 - r5^3)/(r3*r4 - r5^2))

(r3, r4, r5) = (5, 4, 3) .// 2
(-r3*r4^2/(r3*r4 - r3*r5 - r5^2), r5*(r3 + r5)*(r4 + r5)/(r3*r4 - r5^2), -2*sqrt(r3)*sqrt(r4)*r5*sqrt(r3 + r5)*sqrt(r4 + r5)/(r3*r4 - r5^2), (r3^2*r4 + r3*r4*r5 - r4*r5^2 - r5^3)/(r3*r4 - r5^2))

   (10//1, 84//11, -9.127200289462643, 97//22)

(-r3*r4^2/(r3*r4 - r3*r5 - r5^2), r5*(r3 + r5)*(r4 + r5)/(r3*r4 - r5^2), 2*sqrt(r3)*sqrt(r4)*r5*sqrt(r3 + r5)*sqrt(r4 + r5)/(r3*r4 - r5^2), (r3^2*r4 + r3*r4*r5 - r4*r5^2 - r5^3)/(r3*r4 - r5^2))

   (10//1, 84//11, 9.127200289462643, 97//22)

解は 2 組得られるが,両者の違いは,中円の中心座標が右のものか左のものかの違いだけである。

大円の直径は 2(-r3*r4^2/(r3*r4 - r3*r5 - r5^2)) = 20 寸である。
術では 「3920 を 196 で割ればよい」とあるが,なぜそのような数値が出てくるか不明である。

   r1 = 10;  r2 = 7.63636;  x2 = 9.1272;  y2 = 4.40909
   大円の直径 = 20;  中円の直径 = 15.2727

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r3, r4, r5) = (5, 4, 3) .// 2
   (r1, r2, x2, y2) = (-r3*r4^2/(r3*r4 - r3*r5 - r5^2), r5*(r3 + r5)*(r4 + r5)/(r3*r4 - r5^2), 2*sqrt(r3)*sqrt(r4)*r5*sqrt(r3 + r5)*sqrt(r4 + r5)/(r3*r4 - r5^2), (r3^2*r4 + r3*r4*r5 - r4*r5^2 - r5^3)/(r3*r4 - r5^2))
   @printf("r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", r1, r2, x2, y2)
   @printf("大円の直径 = %g;  中円の直径 = %g\n", 2r1, 2r2)
   plot()
   circle(0, r3 + 2r5 + 2r4 + r1, r1)
   circle(x2, y2, r2, :blue)
   circle(-x2, y2, r2, :blue)
   circle(0, 0, r3, :orange)
   circle(0, r3 + 2r5 + r4, r4, :magenta)
   circle(0, r3 + r5, r5, :green)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, r3 + 2r5 + 2r4 + r1, " r3+2r5+2r4+r1\n\n      大円:r1", :red)
       point(x2, y2, "中円:r2,(x2,y2)", :blue, :center, delta=-delta)
       point(0, r3 + 2r5 + r4, "乙円:r4  r3+2r5+r4 ", :magenta, :right)
       point(0, r3 + r5, "丙円:r5  r3+r5 ", :green, :right)
       point(0, 0, "小円:r5 ", :orange, :right)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

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

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