裏 RjpWiki

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

算額(その1266)

2024年09月01日 | Julia

算額(その1266)

百三十 現存するが奉納場所不明 明治14年(1881)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:円5個,菱形

菱形の中に大円 1 個,中円 2 個,小円 2 個を容れる。菱長(菱形の対角線の長い方)と大円の直径が与えられたとき,小円の直径を求めるすべを述べよ。

菱長と菱平(菱形の対角線の短い方)を 2a, 2b
大円の半径と中心座標を r1, (0, 0)
中円の半径と中心座標を r2, (r2, 0)
小円の半径と中心座標を r3, (2r2 + r3, 0)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     r1::positive, r2::positive, r3::positive;
eq1 = r1/a - b /sqrt(a^2 + b^2)
eq2 = r2/(a - r2) - r1/a
eq3 = r3/(a - 2r2 - r3) - r1/a
res = solve([eq1, eq2, eq3], (b, r2, r3))[1]

   (a*r1*sqrt(1/(a - r1))/sqrt(a + r1), a*r1/(a + r1), a*r1*(a - r1)/(a^2 + 2*a*r1 + r1^2))

小円の半径 r3 はたいして簡約化できない。

@syms d
r3 = apart(res[3], d) |> factor
r3 |> println

   -a*r1*(-a + r1)/(a + r1)^2

小円の半径 r3 は,菱長 a,大円の半径 r1 の関数である。
これは,「術」の「大円径/(大円径 + 菱長)^2 * 菱長*(菱長 - 大円径)」に一致する。

たとえば,a = 5, r1 = 3 のとき,r3 = -a*r1*(-a + r1)/(a + r1)^2 = 15/32 = 0.46875 である。
菱長,大円の直径が 10, 6 のとき,小円の直径は 0.9375 である。

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

   a = 5;  r1 = 3;  b = 3.75;  r2 = 1.875;  r3 = 0.46875

function draw(a, r1, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (b, r2, r3) = (a*r1*sqrt(1/(a - r1))/sqrt(a + r1), a*r1/(a + r1), a*r1*(a - r1)/(a^2 + 2*a*r1 + r1^2))
   @printf("a = %g;  r1 = %g;  b = %g;  r2 = %g;  r3 = %g\n",
       a, r1, b, r2, r3)
   @printf("菱長,大円の直径が %g, %g のとき,小円の直径は %g である。\n", 2a, 2r1, 2r3)
   plot([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:blue, lw=0.5)
   circle(0, 0, r1)
   circle2(r2, 0, r2, :green)
   circle2(2r2 + r3, 0, r3, :magenta)
   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, 0, "大円:r1,(0,0) ", :red, :right, delta=-delta/2)
       point(r2, 0, "中円:r2,(r2,0)", :green, :center, delta=-delta/2)
       point(2r2 + r3, 0, "小円:r3\n(2r2+r3,0)", :magenta, :left, delta=-5delta)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(0, b, " b", :blue, :left, :bottom, delta=delta/2)
       xlims!(-a - 2delta, a + 8delta)
   end
end;

draw(10/2, 6/2, true)

 

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

算額(その1265)

2024年09月01日 | Julia

算額(その1265)

百二十六 群馬県倉渕村水沼 蓮華院 明治11年(1878)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:円3個,正方形,斜線

正方形を折り返した図形の中に,大円,中円,小円を容れる。大円の直径が 1.5 寸のとき,中円の直径はいかほどか。

正方形の一辺の長さを a とする。「折り返した」ことによる線分の長さの関係から,大円と小円に接する斜線と正方形の一辺の交点座標 (0, a/2), (b, 0), (c, 0) において,b = a/4, c = 2a/3 は事前に簡単に計算できる。
大円の半径と中心座標を r1, (a - r1, r1)
中円の半径と中心座標を r2, (x2, y2)
小円の半径と中心座標を r3, (r3, r3)
とおき,以下の連立方程式を解く。
なお,中円の入っている直角三角形は折る前の直角三角形と合同で,その中にある(灰色で示す)中円の直径もまた簡単に計算することはできる。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive, r1::positive,
     r2::positive, x2::positive, y2::positive, r3::positive;
b = a/4
c = 2a/3
eq1 = dist2(a/4, 0, a, a, a - r1, r1, r1)

eq2 = a + a/2 - sqrt(a^2/4 + a^2) - 2r2
eq3 = dist2(0, a/2, a, a, x2, y2, r2)
eq4 = dist2(a/4, 0, a, a, x2, y2, r2)

eq5 = dist2(a/4, 0, a, a, r3, r3, r3);

やり方は色々あるが SymPy の能力と簡便さを考慮すると,以下のように連立方程式を解くのが無難である。

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

res[1] |> println
res[2] |> println
res[3] |> println

   a/4
   a*(3 - sqrt(5))/4
   a/6

大円,中円,小円の半径は,正方形の一辺の長さの 1/4 倍,(3 - √5)/4 倍,1/6 倍である。

中円の半径は,大円の半径の (3 - √5) = 0.7639320225002102 倍である。
大円の直径が 1.5 寸のとき,中円の直径は 1.5*(3 - √5) = 1.1458980337503153 寸である。

x2, y2 については,簡約化できる。

res[4] |> sympy.sqrtdenest |> simplify |> println
res[5] |> sympy.sqrtdenest |> simplify |> println

   a*(sqrt(5) + 5)/20
   a*(25 - 7*sqrt(5))/20

function draw(d1, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1  # 任意の値
   b = a/4
   c = 2a/3
   r1 = a/4
   r2 = a*(3 - √5)/4
   r3 = a/6
   x2 = a*(√5 + 5)/20
   y2 = a*(25 - 7√5)/20
   @printf("a = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x2 = %g;  y2 = %g\n",
       a, r1, r2, r3, x2, y2)
   p = r2/r1
   @printf("中円の直径は大円の直径の %g 倍である。\n", p)
   @printf("大円の直径が %g のとき,中円の直径は %g である。\n", d1, d1*p)
   plot([a, 0, 0], [a, a, a/2], color=:blue, lw=0.5, linestyle=:dashdot)
   plot!([a, 0, 0, a, a, 0, c], [a, a/2, 0, 0, a, a/2, 0], color=:blue, lw=0.5)
   segment(b, 0, a, a, :blue)
   circle(a - r1, r1, r1)
   circle(x2, y2, r2, :green)
   circle(r3, r3, r3, :magenta)
   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, delta=-delta/2)
       point(x2, y2, "中円:r2,(x2,y2)", :green, :center, delta=-delta/2)
       circle(r2, a - r2, r2, :gray70)
       point(r2, a - r2, "中円:r2,(r2,a-r2)", :gray70, :center, delta=-delta/2)
       point(r3, r3, "小円:r3,(r3,r3)", :magenta, :center, delta=-delta/2)
       point(0, 0, "(0,0)", :blue, :center, delta=-delta)
       point(a, 0, "(a,0)", :blue, :center, delta=-delta)
       point(b, 0, "(b,0)", :blue, :center, delta=-delta)
       point(c, 0, "(c,0)", :blue, :center, delta=-delta)
       point(0, a/2, "(0,a/2) ", :blue, :right, :vcenter)
       point(0, a, "(0,a) ", :blue, :right, :vcenter)
       plot!(showaxis=false, xlims=(-10delta, a + 3delta))
   end
end;

draw(1.5, true)

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

算額(その1264)

2024年09月01日 | Julia

算額(その1264)

百二十六 群馬県倉渕村水沼 蓮華院 明治11年(1878)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:円3個,半円2個,直角三角形

直角三角形の中に互いに交わる 2 つの弧と,大円 1 個,小円 2 個を容れる。大円の直径が 16 寸のとき,鈎(直角三角形の直角を挟む 2 辺のうちの短い方の辺)はいかほどか。

この算額の図形は,弧の中心が何かを全く述べていない。問題をわかりやすく書いた「算額(その480)」と本質は同じである。
https://blog.goo.ne.jp/r-de-r/e/72ee24fcf2659d1796849823dda2fe77

下に凸の弧は鈎を直径とする円の一部,上に凸の弧は股を直径とする円の一部である。

本ブログは,図形のパラメータをすべて求めることをポリシーにしているので,「算額(その480)」の結果から,本問の答えを導くことができる。

名称・記号を対応付けると,「鈎は大円の直径の 9/4 倍」である。
大円の直径が 16 寸のとき,鈎は 16*9/4 = 36 寸である。

include("julia-source.txt");

function draw(r4, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (3*r4, 9*r4/4)
   r3 = 3*r4/4
   y31 = 2*sqrt(3)*sqrt(r3)*sqrt(r4)
   (x3, y32) = (6*r4/5, 27*r4/20)
   x4 = r1
   p = Rational(r2/r4)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  y31 = %g;  y32 = %g;  r4 = %g;  x4 = %g\n",  r1, r2, r3, x3, y31, y32, r4, x4)
   @printf("鈎は甲円の直径の %d/%d 倍である。\n", numerator(p), denominator(p))
   plot()
   circle(r1, 0, r1, beginangle=0, endangle=180)
   circle(0, r2, r2, :blue, beginangle=-90, endangle=90)
   circle(x3, y32, r3, :green)
   circle(r3, y31, r3, :green)
   circle(x4, r4, r4, :orange)
   segment(0, 2r2, 2r1, 0, :black)
   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, 0, " r1", :red, :left, :bottom, delta=delta/2)
       point(0, r2, " r2", :blue, :left, :vcenter)
       point(r3, y31, "小円:r3\n(r3,y31)", :green, :center, delta=-delta/2)
       point(x3, y32, "小円:r3\n(x3,y32)", :green, :center, delta=-delta/2)
       point(x4, r4, "大円:r4,(x4,r4)", :orange, :center, delta=-delta/2)
       point(0, 2r2, " 鈎", :black, :left, :bottom, delta=delta)
       plot!([0, 0, 2r1], [2r2, 0, 0], color=:black, lw=0.5)
   end
end;

draw(16/2, true)

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

算額(その480)

2024年09月01日 | Julia

算額(その480)

2024/09/01 改訂
解析解を求めるようにした

宮城県丸森町小斎日向 鹿島神社 明治13年
徳竹亜紀子,谷垣美保: 2021年度の算額調査,仙台高等専門学校名取キャンパス 研究紀要,第 58 号, p.7-28, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2022/03/kiyo2022-2.pdf
キーワード:円3個,半円2個,直角三角形

直角三角形の 2 つの辺それぞれを直径とする大半円と中半円があり,甲円 1 個,乙円 2 個が入っている。甲円の直径が 12 寸のとき,乙円の直径はいかほどか。

この問題(図)は,「算額(その373)」に1個の円(甲円)を加えたものである。

大半円の半径と中心座標を r1, (r1, 0)
中半円の半径と中心座標を r2, (0, r2)
乙円の半径と中心座標を r3, (r3, y31) および (x3, y32)
甲円の半径と中心座標を (x4, y4)
とおき,以下の連立方程式を(SymPy では一度に解けないので)逐次的に解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive, y31::positive,
     x3::positive, y32::positive, r4::positive, x4::positive;
x4 = r1
eq1 = (r1 - r3)^2 + y31^2 - (r1 + r3)^2
eq2 = (r1 - x3)^2 + y32^2 - (r1 - r3)^2
eq3 = x3^2 + (r2 - y32)^2 - (r2 - r3)^2
eq4 = dist2(0, 2r2, 2r1, 0, r3, y31, r3)
eq5 = dist2(0, 2r2, 2r1, 0, x4, r4, r4)
eq6 = x4^2 + (r2 - r4)^2 - (r2 + r4)^2;

eq5, eq6 を解き,r1, r2 を求める。

res1 = solve([eq5, eq6], (r1, r2))[1]

   (3*r4, 9*r4/4)

eq1 に r1, r2 を代入し,y31 を求める。

eq11 = eq1(r1 => 3r4, r2 => 9r4/4) |> simplify
ans_y31 = solve(eq11, y31)[1]
ans_y31 |> println

   2*sqrt(3)*sqrt(r3)*sqrt(r4)

eq4 に r1, r2, y31 を代入し,r3 を求める。

eq14 = eq4(r1 => 3r4, r2 => 9r4/4, y31 => ans_y31)/(9r4^2/4) |> simplify
ans_r3 = solve(eq14, r3)[1]  # 1 of 3
ans_r3 |> println

   3*r4/4

乙円の半径 r3 は,甲円の半径 r4 の 3/4 倍である。
甲円の直径が 12 寸のとき,乙円の直径は 9 寸である。

なお,このとき,直角三角形の辺の長さの比は 3:4:5 になる。

問に答えるためにはここまでで十分であるが,図を描くために残りのパラメータも求める。

eq2, eq3 に r1, r2, r3 を代入し,x3, y32 を求める。

eq12 = eq2(r1 => 3r4, r2 => 9r4/4, r3 => ans_r3) |> expand
eq13 = eq3(r1 => 3r4, r2 => 9r4/4, r3 => ans_r3) |> expand
res2 = solve([eq12, eq13], (x3, y32))[1]

   (6*r4/5, 27*r4/20)

function draw(r4, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, r2) = (3*r4, 9*r4/4)
    r3 = 3*r4/4
    y31 = 2*sqrt(3)*sqrt(r3)*sqrt(r4)
    (x3, y32) = (6*r4/5, 27*r4/20)
    x4 = r1
    p = Rational(r3/r4)
    @printf("r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  y31 = %g;  y32 = %g;  r4 = %g;  x4 = %g\n",  r1, r2, r3, x3, y31, y32, r4, x4)
    @printf("乙円の直径は甲円の直径の %d/%d 倍である。\n", numerator(p), denominator(p))
    plot()
    circle(r1, 0, r1, beginangle=0, endangle=180)
    circle(0, r2, r2, :blue, beginangle=-90, endangle=90)
    circle(x3, y32, r3, :green)
    circle(r3, y31, r3, :green)
    circle(x4, r4, r4, :orange)
    segment(0, 2r2, 2r1, 0, :black)
    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, 0, " r1", :red, :left, :bottom, delta=delta/2)
        point(0, r2, " r2", :blue, :left, :vcenter)
        point(r3, y31, "乙円:r3\n(r3,y31)", :green, :center, delta=-delta/2)
        point(x3, y32, "乙円:r3\n(x3,y32)", :green, :center, delta=-delta/2)
        point(x4, r4, "甲円:r4,(x4,r4)", :orange, :center, delta=-delta/2)
        plot!([0, 0, 2r1], [2r2, 0, 0], color=:black, lw=0.5)
    end
end;

draw(12/2, true)

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

算額(その1263)

2024年09月01日 | Julia

算額(その1263)

百二十六 群馬県倉渕村水沼 蓮華院 明治11年(1878)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:円3個,外円,正方形2個

外円の中に,正方形 2 個,円 2 個を容れる。外円の直径が 10 寸のとき,等円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (R - r, 0)
とおき,以下の方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive
eq1 = R - r ⩵ √Sym(2)r
res = solve(eq1, r)[1] |> simplify
res |> println

   R*(-1 + sqrt(2))

等円の半径 r は,外円の半径 R の (√2 - 1) 倍である。
外円の直径が 10 寸のとき,等円の直径は 10(√2 - 1) = 4.142135623730951 寸である。

function draw(R, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = (√2 - 1)R
   @printf("外円の直径が %g のとき,等円の直径は %g である。\n", 2R, 2r)
   plot([0, R/2, 0, -R/2, 0], [0, R/2, R, R/2, 0], color=:green, lw=0.5)
   plot!([0, R/2, 0, -R/2, 0], [0, -R/2, -R, -R/2, 0], color=:green, lw=0.5)
   circle(0, 0, R, :blue)
   circle2(R - r, 0, r)
   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(R/2, R/2, " (R/2,R/2)", :green, :left, :vcenter)
       point(R - r, 0, "(R-r,0)", :red, :center, delta=-delta/2)
       point(R, 0, " R", :blue, :left, delta=-delta/2)
   end
end;

draw(10/2, true)

 

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

算額(その1262)

2024年09月01日 | Julia

算額(その1262)

百二十六 群馬県倉渕村水沼 蓮華院 明治11年(1878)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:円5個,正三角形

正三角形の中に,等円 5 個を容れる。正三角形の一辺の長さが与えられとき,等円の直径を求めよ。

注:5 個の円の中心は正五角形の頂点である。

正三角形の一辺の長さが 2a,等円の半径を r,右上の等円の中心座標を (r + 2r*cos(2π/5), r + 2r*sin(2π/5)) とおき,「斜辺と等円の距離が等円の半径に等しい」とする方程式を解く。

include("julia-source.txt");

using SymPy
@syms a, r
eq1 = dist2(a, 0, 0,√Sym(3)a, r + 2r*cos(2PI/5), r + 2r*sin(2PI/5), r)
res = solve(eq1, r)[2]  # 2 of 2
res |> println

   a*(-2*sqrt(3) + 3 + sqrt(6*sqrt(5) + 30) + 3*sqrt(5))/(2*sqrt(3) + sqrt(6)*sqrt(sqrt(5) + 5) + 2*sqrt(2)*sqrt(sqrt(5) + 5) + 2*sqrt(15) + 8 + 4*sqrt(5) + sqrt(30)*sqrt(sqrt(5) + 5))

等円の半径は a の定数倍であるが,その正確な値は長い式で表される。式と言っても定数式なので,まとめると a の 0.2248064842773 倍ということである。

res.evalf() |> println

   0.2248064842773*a

正三角形の一辺の長さが 3 寸のとき,等円の直径は 0.674419452831899 寸である。

2res(a => 3/2).evalf() |> println
2(3/2)*0.2248064842773 |> println

   0.674419452831899
   0.6744194528319

function draw(a, more=false)
   pyplot(size=(650, 650), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = a*(-2√3 + 3 + sqrt(6√5 + 30) + 3√5)/(2√3 + √6sqrt(√5 + 5) + 2√2sqrt(√5 + 5) + 2√15 + 8 + 4√5 + √30*sqrt(√5 + 5))
   @printf("正三角形の一辺の長さが %g のとき,等円の直径は %g である。\n", 2a, 2r)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:blue, lw=0.5)
   circle2(r, r, r)
   circle2(r + 2r*cos(2π/5), r + 2r*sin(2π/5), r)
   circle(0, r + 2r*sin(2π/5) + 2r*sin(π/5), r)
   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, √3a, " √3a", :blue, :left, :vcenter)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(r + 2r*cos(2π/5), r + 2r*sin(2π/5), "(r+2r*cos(2π/5),\nr+2r*sin(2π/5))", :red, :center)
   end
end;

draw(3/2, true)

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

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

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