裏 RjpWiki

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

算額(その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) | トップ | 算額(その1264) »
最新の画像もっと見る

コメントを投稿

Julia」カテゴリの最新記事