裏 RjpWiki

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

算額(その755)

2024年03月05日 | Julia

算額(その755)

宮城県白石市小原字馬頭山 三瀧神社奉納算額(萬蔵稲荷神社所蔵) 明治8年
徳竹亜紀子,谷垣美保:宮城県白石市小原地区の算額調査,仙台高等専門学校名取キャンパス研究紀要,第57号,2021.

https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2021/04/No57_2.pdf

外円の中に,甲円,乙円,大円,中円,小円が入っている。外円,甲円,乙円の直径がそれぞれ 64 寸,52 寸 9 分,36 寸 1 分のとき,中円の直径はいかほどか。

中円は甲円に内接し,乙円と外接している。小円は甲円と外接し,乙円に内接している。
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, r1 - R)
乙円の半径と中心座標を r2, (0, R - r2)
大円の半径と中心座標を r3, (0, r1 - r2)
中円の半径と中心座標を r4, (r4, y4)
小円の半径と中心座標を r5, (r5, y5)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive, r3::positive,
     r4::positive, y4::negative, r5::positive, y5::positive
#(R, r1, r2) = (640, 529, 361) .// 20
eq1 = r4^2 + (r1 - R - y4)^2 - (r1 - r4)^2
eq2 = r5^2 + (y5 - R + r2)^2 - (r2 - r5)^2
eq3 = r4^2 + (R - r2 - y4)^2 - (r2 + r4)^2
eq4 = r5^2 + (y5 - r1 + R)^2 - (r1 + r5)^2
eq5 = 2r1 + 2r2 - 2r3 - 2R
res = solve([eq1, eq2, eq3, eq4, eq5], (r3, r4, y4, r5, y5))

   4-element Vector{NTuple{5, Sym{PyCall.PyObject}}}:
    (-R + r1 + r2, (R*r1 - R*r2 + (-2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), -2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2), -(R*r1 - R*r2 + (-2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), -2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))
    (-R + r1 + r2, (R*r1 - R*r2 + (-2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), -2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2), -(R*r1 - R*r2 + (2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), 2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))
    (-R + r1 + r2, (R*r1 - R*r2 + (2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), 2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2), -(R*r1 - R*r2 + (-2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), -2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))
    (-R + r1 + r2, (R*r1 - R*r2 + (2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), 2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2), -(R*r1 - R*r2 + (2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2), 2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))

for i = 1:4
   print("#$i: ")
   for j = 1:5
       print((res[i][j](R => 64.0/2, r1 => 52.9/2, r2 => 36.1/2)).evalf(), "; ")
   end
   println()
end

   #1: 12.5000000000000; 12.0000000000000; -13.6000000000000; -12.0000000000000; -13.6000000000000; 
   #2: 12.5000000000000; 12.0000000000000; -13.6000000000000; 5.21297815932332; 25.6808988764045; 
   #3: 12.5000000000000; -5.21297815932332; 25.6808988764045; -12.0000000000000; -13.6000000000000; 
   #4: 12.5000000000000; -5.21297815932332; 25.6808988764045; 5.21297815932332; 25.6808988764045; 

4 組の解が得られるが,2 番目のものが適解である。
中円の半径を表す数式はかなり長いものであり,SymPy ではうまく簡約化できない。

res[2][2] |> println

   (R*r1 - R*r2 + (-2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(-R + r1 + r2)/(r1 + r2) + R*(r1 - r2)/(r1 + r2))*(-2*R + r1 + r2))/(r1 + r2)

外円,甲円,乙円の直径がそれぞれ 64 寸,52 寸 9 分,36 寸 1 分のとき,中円の直径は 24 寸である。

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

外円の直径 = 64;  甲円の直径 = 52.9;  乙円の直径 = 36.1
大円の直径 = 25;  小円の直径 = 10.426;  r3 = 12.5;  r4 = 12;  y4 = -13.6;  r5 = 5.21298;  y5 = 25.6809

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1, r2) = (640, 529, 361) .// 20
   # (r3, r4, y4, r5, y5) = (25/2, 12, -68/5, 41292/7921, 11428/445)
   t0 = r1 + r2
   t1 = 2sqrt(R*r1*r2*(-R + t0))/t0
   t2 = R*(r1 - r2)/t0
   (r3, r4, y4, r5, y5) = (
       t0 - R,
       (R*(r1 - r2) + (t2 - t1)*(t0 - 2R))/t0,
       t2 - t1,
       -(R*(r1 - r2) + (t2 + t1)*(t0 - 2R))/t0,
       t2 + t1)
   @printf("外円の直径 = %g;  甲円の直径 = %g;  乙円の直径 = %g\n", 2R, 2r1, 2r2)
   @printf("中円の直径 = %g;  大円の直径 = %g;  小円の直径 = %g;  r3 = %g;  r4 = %g;  y4 = %g;  r5 = %g;  y5 = %g\n",
       2r4, 2r3, 2r5, r3, r4, y4, r5, y5)
   plot()
   circle(0, 0, R, :black)
   circle(0, r1 - R, r1)
   circle(0, R - r2, r2, :blue)
   circle(0, r1 - r2, r3, :green)
   circle2(r4, y4, r4, :magenta)
   circle2(r5, y5, r5, :orange)
   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, r1 - R, " 甲円:r1,(0,r1-R)", :red, :left, :vcenter)
       point(0, R - r2, "乙円:r2 \n(0,R-r2) ", :blue, :right, :vcenter)
       point(0, r1 - r2, " 大円:r3\n (0,r1-r2)", :green, :left, :vcenter)
       point(r4, y4, "中円:r4,(r4,y4)", :magenta, :center, delta=-delta)
       point(r5, y5, "小円:r5\n(r5,y5)", :orange, :center)
       point(R, 0, " R", :black, :left, :bottom, delta=delta/2)
   end
end;

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

算額(その754)

2024年03月05日 | Julia

算額(その754)

草加市金明町 旭神社 寛政11年(1799)
山口正義:やまぶき3,第45号

https://yamabukiwasan.sakura.ne.jp/ymbk45.pdf

鈎股弦(直角三角形)において,弦の二乗と中鉤の和が 27 寸 4 分,股の二乗と内接円の直径の和が 18 寸のとき,弦はいくつか。

直角三角形の直角を挟む二辺の短い方を「鈎」,長い方を「股」,直角の対辺を「弦」と呼ぶ。
中鈎とは直角の頂点から弦へおろした垂線である。
直角三角形に内接する円の直径と鈎・股・弦には eq3 で示す関係がある。鈎・股・弦の間の関係式 eq4 はピタゴラスの定理。未知数が 5 個なので,条件式 eq5 を付け加える。
図を描くために,中鈎と弦の交点(中鈎の脚)の座標を (x, y) として以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r::positive, 中鈎::positive, x::positive, y::positive,
     鈎::positive, 股::positive, 弦::positive
eq1 = 弦^2 + 中鈎 - 274//10
eq2 = 股^2 + 2r - 18
eq3 = 鈎 + 股 - 弦 - 2r
eq4 = 鈎^2 + 股^2 - 弦^2
eq5 = 弦*中鈎 - 鈎*股
res = solve([eq1, eq2, eq3, eq4, eq5], (鈎, 股, 弦, 中鈎, r))

   1-element Vector{NTuple{5, Sym{PyCall.PyObject}}}:
    (3, 4, 5, 12/5, 1)

(鈎, 股, 弦, 中鈎, r) = res[1]
eq6 = x^2 + y^2 - 中鈎^2
eq7 = (鈎 - y)/x - 鈎/股
res2 = solve([eq6, eq7], (x, y))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (36/25, 48/25)

甲円の直径は 18 寸である。

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

   弦 = 5, 鈎 = 3;  股 = 4;  中鈎 = 2.4;  r = 1;  x = 1.44;  y = 1.92

なお,引用元にも書いてあるが,和算的には以下のような12次方程式を解くことになるそうだが,これでは「弦」しか求まらない。

eq = x^12 - 109.6*x^10 -6x^9 + 4472.56*x^8 + 512.4x^7 - 79881.696*x^6 - 14547.04*x^5 + 504392.2576*x^4 + 137662.816*x^3 + 475181.92*x^2 - 12182.4*x+ 104976
solve(eq)

   2-element Vector{Sym{PyCall.PyObject}}:
    5.00000000000000
    5.46220825640114

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, 弦, 中鈎, r) = (3, 4, 5, 12/5, 1)
   (x, y) = (36/25, 48/25)
   @printf("弦 = %g, 鈎 = %g;  股 = %g;  中鈎 = %g;  r = %g;  x = %g;  y = %g\n", 弦, 鈎, 股, 中鈎, r, x, y)
   plot([0, 股, 0, 0 ], [0, 0, 鈎, 0], color=:black, lw=0.5)
   segment(0, 0, x, y, :green)
   circle(r, r, 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, r, "内接円:r,(r,r)", :red)
       point(0, 鈎, " 鈎", :black, :left, :bottom)
       point(股, 0, " 股", :black, :left, :bottom)
       point(x, y, "(x,y)", :black, :left, :bottom)
       point(0.7, 1.5, "中鈎", :green, mark=false)
       point(2.5, 1.3, "弦", :black, mark=false)
   end
end;

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

算額(その753)

2024年03月05日 | Julia

算額(その753)

四二 浦和市西堀(現さいたま市桜区西堀) 氷川神社 嘉永5年(1852)
埼玉県立図書館:埼玉県史料集 第二集『埼玉の算額』,昭和44年,誠美堂印刷所,埼玉県与野市.

さいたま市桜区 西堀氷川神社 嘉永5年(1852)
山口正義:やまぶき3,第45号

https://yamabukiwasan.sakura.ne.jp/ymbk45.pdf

直線上に甲円と乙円が載っており,互いに接している。甲円と乙円の中心を結ぶ直線を弦として,座標軸に平行な二辺(短い方を鈎,長い方を股と呼ぶ)とする直角三角形を考える。乙円の直径と股の和が 20 寸,直角三角形に内接する円(黒円と名付ける)の直径が 4 寸,のとき,甲円の直径はいかほどか。

甲円の半径と中心座標を r1, (0, r1)
乙円の半径と中心座標を r2, (股, r2)
黒円の半径と中心座標を r3, (股 - r3, r1 - r3)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive,
     鈎::positive, 股::positive, 弦::positive
r3 = 4//2
鈎 = r1 - r2 
股 = 2sqrt(r1*r2)
弦 = r1 + r2
eq1 = 2r2 + 股 - 20
eq2 = 鈎 + 股 - 弦 - 2r3
res = solve([eq1, eq2], (r1, r2))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (9, 4)

甲円の直径は 18 寸である。

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

r1 = 9;  r2 = 4;  r3 = 2;  股 = 12

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (9, 4)
   r3 = 4//2
   鈎 = r1 - r2 
   股 = 2sqrt(r1*r2)
   @printf("甲円の直径 = %g, r1 = %g;  r2 = %g;  r3 = %g;  股 = %g\n", 2r1, r1, r2, r3, 股)
   plot()
   circlef(股 - r3, r1 - r3, r3, :gray50)
   plot!([0, 股, 股, 0], [r1, r2, r1, r1], color=:blue, lw=0.5)
   circle(0, r1, r1)
   circle(股, r2, r2, :green)
   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, r1, " 甲円:r1,(0,r1)", :red, :center, :bottom, delta=delta)
       point(股, r2, " 乙円:r2,(股,r2)", :green, :center, delta=-delta)
       point(股 - r3, r1 - r3, "黒円:r3", :yellow, :center, :bottom, delta=delta)
       point(6, r1, "股", :blue, :center, :bottom, delta=delta, mark=false)
       point(6, 3r1/4, "弦", :blue, :center, delta=-3delta, mark=false)
       point(股, 3r1/4, " 鈎", :blue, :left, :vcenter, mark=false)
   end
end;

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

算額(その752)

2024年03月05日 | Julia

算額(その752)

東京都稲城市 穴沢天神社 明治10年
山口正義:やまぶき,第15号

https://yamabukiwasan.sakura.ne.jp/ymbk15.pdf

重なった等円と長方形で区切られた領域に大円と小円が入っている。大円と小円の直径がそれぞれ 5 分 7 厘,2 分  1厘のとき,長方形の長辺の長さはいかほどか。

長方形の長辺と短辺の長さを 2b, a
大円の半径と中心座標を r1, (0, r0)
小円の半径と中心座標を r2, (0, r2), (0, a + r2)
等円の半径と中心座標を r0, (r0 - r1, r0)
とおいて以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r0::positive, r1:: poaitive,
     r2::positive, a::positive, b::positive
eq1 = (r0 - r1)^2 + (r0 - r2)^2 - (r0 + r2)^2
eq2 = (r0 - r1)^2 + (a - r0 + r2)^2 - (r0 - r2)^2
eq3 = b - (2r0 - r1)

res = solve([eq1, eq2, eq3], (a, b, r0))

   4-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (r1 - 2*sqrt(r2)*sqrt(r1 + r2) + r2 - sqrt(-4*r1*r2 + 8*r2^(3/2)*sqrt(r1 + r2) - 8*r2^2 + (r1 - 2*sqrt(r2)*sqrt(r1 + r2) + r2)^2), r1 - 4*sqrt(r2)*sqrt(r1 + r2) + 4*r2, r1 - 2*sqrt(r2)*sqrt(r1 + r2) + 2*r2)
    (r1 - 2*sqrt(r2)*sqrt(r1 + r2) + r2 + sqrt(-4*r1*r2 + 8*r2^(3/2)*sqrt(r1 + r2) - 8*r2^2 + (r1 - 2*sqrt(r2)*sqrt(r1 + r2) + r2)^2), r1 - 4*sqrt(r2)*sqrt(r1 + r2) + 4*r2, r1 - 2*sqrt(r2)*sqrt(r1 + r2) + 2*r2)
    (r1 + 2*sqrt(r2)*sqrt(r1 + r2) + r2 - sqrt(-4*r1*r2 - 8*r2^(3/2)*sqrt(r1 + r2) - 8*r2^2 + (r1 + 2*sqrt(r2)*sqrt(r1 + r2) + r2)^2), r1 + 4*sqrt(r2)*sqrt(r1 + r2) + 4*r2, r1 + 2*sqrt(r2)*sqrt(r1 + r2) + 2*r2)
    (r1 + 2*sqrt(r2)*sqrt(r1 + r2) + r2 + sqrt(-4*r1*r2 - 8*r2^(3/2)*sqrt(r1 + r2) - 8*r2^2 + (r1 + 2*sqrt(r2)*sqrt(r1 + r2) + r2)^2), r1 + 4*sqrt(r2)*sqrt(r1 + r2) + 4*r2, r1 + 2*sqrt(r2)*sqrt(r1 + r2) + 2*r2)

4 組の解が得られるが,最後のものが適解である。

i = 4
@printf("a = %g;  b = %g;  r0 = %g\n",
       res[i][1](r1 => 0.57/2, r2 => 0.21/2),
       res[i][2](r1 => 0.57/2, r2 => 0.21/2),
       res[i][3](r1 => 0.57/2, r2 => 0.21/2))

   a = 1.29841;  b = 1.51444;  r0 = 0.899722

長方形の長辺の長さは 2(r1 + 4*sqrt(r2)*sqrt(r1 + r2) + 4*r2) = 3.028888507587845である。
算額の答では「直長三寸〇九厘二九有奇」とあるが,写し間違いであろう。

r1 = 0.57/2
r2 = 0.21/2
2(r1 + 4*sqrt(r2)*sqrt(r1 + r2) + 4*r2)

   3.028888507587845

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

r1 = 0.285;  r2 = 0.105;  a = 1.29841;  b = 1.51444;  r0 = 0.899722

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (57, 21) .// 200
   a = 1.29841;  b = 1.51444;  r0 = 0.899722
   @printf("長方形の長辺の長さ = %g;  r1 = %g;  r2 = %g;  a = %g;  b = %g;  r0 = %g\n", 2b, r1, r2, a, b, r0)
   plot([b, b, -b, -b, b], [0, a, a, 0, 0], color=:black, lw=0.5)
   circle2(r0 -  r1, r0, r0)
   circle(0, r0, r1, :blue)
   circle(0, r2, r2, :green)
   circle(0, a + r2, r2, :green)
   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(r0 - r1, r0, "等円:r0,(r0-r1,r0)", :red, :left, delta=-delta)
       point(0, r0, "大円:r1\n(0,r1)", :blue, :center, delta=-delta)
       point(0, r2, "   小円:r2,(0,r2)", :green, :left, :bottom, delta=delta)
       point(0, a + r2, "   小円:r2,(0,a+r2)", :green, :left, :bottom, delta=delta)
       point(0, a, "a", :black, :center, delta=-delta)
       point(b, 0, " b", :black, :left, :bottom, delta=delta)
   end
end;

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

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

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