裏 RjpWiki

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

算額(その476)

2023年10月25日 | Julia

算額(その476)

宮城県丸森町小斎日向 鹿島神社 弘化5年

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

鈎股弦があり,鈎,股,弦を直径としその中点を中心とする 3 個の円,天円,地円,人円がある。
鈎股弦に内接する全円と,同じ直径を持つもう 1 つの全円は人円,地円,甲円に外接し,天円に内接している。
甲円は全円,地円と外接し,天円と内接している。
乙円は人円と地円の交わる部分にあり,人円と地円に内接する。

乙円の直径が与えられたとき,甲円の直径を求めよ。

鈎股弦は算額の写真を私が測定したところ,各辺が 3:4:5 の直角三角形と推測できる。

鈎, 股, 弦 の長さを 6, 8, 10 とする(偶数になるように設定する)
甲円の半径と中心座標を (r1, x1, y1)
乙円の半径と中心座標を (r2, x2, y2)
天円の半径と中心座標を (r3, x3, y3) = (弦/2, 股/2, 鈎/2)
地円の半径と中心座標を (r4, x4, |y4) = (股/2, 股/2, 0)
人円の半径と中心座標を (r5, x5, y5) = (鈎/2, 0, 鈎/2)
全円の半径と中心座標を (r6, x6, y6) = ((鈎 + 股 - 弦)/2, 鈎 + 股 - 弦)/2, 鈎 + 股 - 弦)/2)
とおいて以下の連立方程式を解く。

以下の連立方程式は条件が 7 個,未知数が 8 個なので,乙円の直径は一意に定まらない。
そこで,乙円の直径が最大になるときの解を求めるため,r1, x1, y1, x, y, r2, x2 について求めるが,それぞれは y3 を変数として含む式で表される。

include("julia-source.txt")

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
     r1::positive, x1::positive, y1::positive,
     r2::positive, x2::positive, y2::positive,
     r3::positive, x3::positive, y3::positive,
     r4::positive, x4::positive, y4::positive,
     r5::positive, x5::positive, y5::positive,
     r6::positive, x6::positive, y6::positive,
     x::positive, y::positive;
(鈎, 股, 弦) = (3, 4, 5).*2
(r3, x3, y3) = (弦//2, 股//2, 鈎//2)  # 天円
(r4, x4, y4) = (股//2, 股//2, 0)  # 地円
(r5, x5, y5) = (鈎//2, 0, 鈎//2)  # 人円
(r6, x6, y6) = fill((鈎 + 股 - 弦)//2, 3)  # 全円
eq1 = (x1 - x3)^2 + (y1 - y3)^2 - (r3 - r1)^2
eq2 = (x1 - x4)^2 + y1^2 - (r1 + r4)^2
eq3 = (x1 - x)^2 + (y1 - y)^2 - (r1 + r6)^2
eq4 = (x - x4)^2 + y^2 - (r6 + r4)^2
eq5 = (x - x3)^2 + (y - y3)^2 - (r3 - r6)^2
eq6 = (x4 - x2)^2 + y2^2 - (r4 - r2)^2
eq7 = x2^2 + (y5 - y2)^2  - (r5 - r2)^2

res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r1, x1, y1, x, y, r2, x2))

   4-element Vector{NTuple{7, Sym}}:
    (3/2, 4 - sqrt(10), 9/2, 4, 6, y2/5 - 8*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 16/5, 4*y2/5 - 2*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 4/5)
    (3/2, 4 - sqrt(10), 9/2, 4, 6, y2/5 + 8*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 16/5, 4*y2/5 + 2*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 4/5)
    (3/2, sqrt(10) + 4, 9/2, 4, 6, y2/5 - 8*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 16/5, 4*y2/5 - 2*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 4/5)
    (3/2, sqrt(10) + 4, 9/2, 4, 6, y2/5 + 8*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 16/5, 4*y2/5 + 2*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 4/5)

4 組の解が求まるが,3 番目のものが適解である。r2 はその 6 番目で,y2/5 - 8*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 16/5 である。

res[3][6] |> println

   y2/5 - 8*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 16/5

この式は,y2 についての微分の式が 0 のとき,r2 が最大値を取ることを意味している。

pyplot(size=(400, 400), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(res[1][6], xlabel="y2", ylabel="r2", xlims=(0, 4))

r2 が最大値を取るときの y2 を求める。

solve(diff(res[3][6]))[1] |> println

   9/5

y2 = 9/5 のとき,r2 が最大値を取る。そのときの実際の r2 は,以下で求められるように r2 = 1 である。

つまり,「鈎, 股, 弦 の長さを 6, 8, 10 とすると」乙円の半径は 1 である。

res[3][6](y2 => 9//5) |> println

   1

甲円の半径は 3/2 である。つまり,甲円の直径は乙円の直径の 1.5 倍である。

res[3][1](y2 => 9//5) |> println

   3/2

   天円: r3 = 5;  x3 = 4;  y3 = 3
   地円: r4 = 4;  x4 = 4;  y4 = 0
   人円: r5 = 3;  x5 = 0;  y5 = 3
   全円: r6 = 2;  x6 = 2;  y6 = 2  x = 4;  y = 6
   甲円: r1 = 1.5;  x1 = 7.16228;  y1 = 4.5
   乙円: r2 = 1;  x2 = 1.6;  y2 = 1.8

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, 弦) = (3, 4, 5).*2
   (r3, x3, y3) = (弦/2, 股/2, 鈎/2)  # 天円
   (r4, x4, y4) = (股/2, 股/2, 0)  # 地円
   (r5, x5, y5) = (鈎/2, 0, 鈎/2)  # 人円
   (r6, x6, y6) = fill((鈎 + 股 - 弦)/2, 3)  # 全円
   y2 = 9/5
   (r1, x1, y1, x, y, r2, x2) = (3/2, sqrt(10) + 4, 9/2, 4, 6, y2/5 - 8*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 16/5, 4*y2/5 - 2*sqrt(6)*sqrt(y2^2 - 3*y2 + 6)/15 + 4/5)
   @printf("天円: r3 = %g;  x3 = %g;  y3 = %g\n", r3, x3, y3)
   @printf("地円: r4 = %g;  x4 = %g;  y4 = %g\n", r4, x4, y4)
   @printf("人円: r5 = %g;  x5 = %g;  y5 = %g\n", r5, x5, y5)
   @printf("全円: r6 = %g;  x6 = %g;  y6 = %g  x = %g;  y = %g\n", r6, x6, y6, x, y)
   @printf("甲円: r1 = %g;  x1 = %g;  y1 = %g\n", r1, x1, y1)
   @printf("乙円: r2 = %g;  x2 = %g;  y2 = %g\n", r2, x2, y2)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
   circle(x3, y3, r3)
   circle(x4, y4, r4, :blue)
   circle(x5, y5, r5, :green)
   circle(x6, y6, r6, :magenta)
   circle(x1, y1, r1, :orange)
   circle(x, y, r6, :magenta)
   circle(x2, y2, r2, :black)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, y1, " 甲円:r1,(x1,y1)", :orange, :center, :bottom, delta=delta)
       #point(x2, y2, "乙円:r2,(x2,y2)", :orange, :center, :bottom, delta=delta)
       point(x3, y3, " 天円:r3,(x3, y3)", :red, :left, :bottom, delta=delta)
       point(x4, 0, "地円:r4,(x4, 0)", :blue, :center, :top, delta=-delta)
       point(0, y5, "人円:r5,(0,y5) ", :green, :right, :bottom, delta=delta)
       point(x6, y6, " 全円:r6,(x6,y6)", :magenta, :left, :bottom)
       point(x, y, "全円:r6,(x,y)", :magenta, :center, :bottom, delta=delta)
       point(x2, y2, " 乙円:r2,(x2,y2)", :black, :left, :top)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;


コメント    この記事についてブログを書く
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする
« 算額(その475) | トップ | 算額(その477) »
最新の画像もっと見る

コメントを投稿

Julia」カテゴリの最新記事