裏 RjpWiki

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

算額(その415)

2023年09月02日 | Julia

算額(その415)

東京都府中市 大國魂神社 明治18年(1885)

早坂平四郎: 算額の一考察,苫小牧工業専門学校紀要
https://www.tomakomai-ct.ac.jp/wp01/wp-content/uploads/2014/06/kiyou5-8.pdf

甲円と乙円が交わっており,共通部分に丙円1こと丁円2個を入れる。それぞれの円は互いに内接・外接している。
甲円,乙円,丙円の直径が与えられたとき,丁円の直径を求めよ。

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

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::negative,
     r4::positive, x4::negative, y4::positive;
eq1 = (r1 - r3 - x4)^2 + y4^2 - (r1 - r4)^2
eq2 = x4^2 + y4^2 - (r3 + r4)^2
eq3 = (r3 - r2 - x4)^2 + y4^2 - (r2 - r4)^2
res = solve([eq1, eq2, eq3], (r4, x4, y4))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (r3*(r1 - r3)*(r2 - r3)/(r1*r2 - r3^2), -r3^2*(r1 - r2)/(r1*r2 - r3^2), -2*sqrt(r1)*sqrt(r2)*r3*sqrt(r1 - r3)*sqrt(r2 - r3)/(r1*r2 - r3^2))
    (r3*(r1 - r3)*(r2 - r3)/(r1*r2 - r3^2), -r3^2*(r1 - r2)/(r1*r2 - r3^2), 2*sqrt(r1)*sqrt(r2)*r3*sqrt(r1 - r3)*sqrt(r2 - r3)/(r1*r2 - r3^2))

最初の組が適解である。

res[1][1] |> simplify

甲円,乙円,丙円の半径 r1, r2, r3 がそれぞれ 55, 45, 17 のとき,丁円の半径は
(r3*(r1 - r3)*(r2 - r3)/(r1*r2 - r3^2) = 8.274473924977126 である。

   r1 = 55;  r2 = 45;  r3 = 17
   r4 = 8.27447;  x4 = -1.32205;  y4 = 25.2399

(r1, r2, r3) = (55, 45, 17)
r3*(r1 - r3)*(r2 - r3)/(r1*r2 - r3^2)

   8.274473924977126

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = (55, 45, 17)
   (r4, x4, y4) = (r3*(r1 - r3)*(r2 - r3)/(r1*r2 - r3^2), -r3^2*(r1 - r2)/(r1*r2 - r3^2), 2*sqrt(r1)*sqrt(r2)*r3*sqrt(r1 - r3)*sqrt(r2 - r3)/(r1*r2 - r3^2))
   @printf("r1 = %g;  r2 = %g;  r3 = %g\n", r1, r2, r3)
   @printf("r4 = %g;  x4 = %g;  y4 = %g\n", r4, x4, y4)
   plot()
   circle(r1 - r3, 0, r1, :blue)
   circle(r3 - r2, 0, r2, :red)
   circle(0, 0, r3, :green)
   circle(x4, y4, r4, :orange)
   circle(x4, -y4, r4, :orange)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r1 - r3, 0, "r1-r3  甲円:r1", :blue, delta=-delta)
       point(r3 - r2, 0, "乙円:r2  r3-r2 ", :red, :right, delta=-delta)
       point(0, 0, "丙円:r3", :green, :center, delta=-delta)
       point(x4, y4, "    丁円:r4,(x4,y4)", :orange, :left, :vcenter)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その414)

2023年09月02日 | Julia

算額(その414)

愛媛県松山市桜谷町 伊佐爾波神社 嘉永3年(1850)

早坂平四郎: 算額の一考察,苫小牧工業専門学校紀要
https://www.tomakomai-ct.ac.jp/wp01/wp-content/uploads/2014/06/kiyou5-8.pdf

直線の上に等円 3 個と大円が互いに接している。さらに,大円の中に4本の弦を引き分割された領域に,等円と同じ半径を持ち弦と接する 2 個の円と,中央の菱形に内接する中円を入れる(菱形の長い方の対角線は,上述の直線と並行である)。

大円の半径と中心座標を a, (0, a)
等円の半径と中心座標を r1, (x1, r1), ( x1 - r1, y1) ただし,y1 = r1(1 + √3)
菱形の中の円の半径と中心座標を r2, (0, a)
弦の1つの端点の座標を (bx, by)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, r1::positive, x1::negative, y1::positive,
     r2::positive, bx::positive;
a = 10
eq1 = x1^2 + (a - r1)^2 - (a + r1)^2
eq2 = (x1 - r1)^2 + (a - r1*(1+sqrt(Sym(3))))^2 - (a + r1)^2
eq3 = distance(-bx, a - sqrt(a^2 - bx^2), a, a, 0, a) - r2^2
eq4 = distance(-bx, a - sqrt(a^2 - bx^2), a, a, 0, r1) - r1^2;
res = solve([eq1, eq2, eq3, eq4], (r1, x1, bx, r2))

   1-element Vector{NTuple{4, Sym}}:
    (210 - 120*sqrt(3), 60 - 40*sqrt(3), -5891176044722613657600*sqrt(-154227 + 89043*sqrt(3))/52283326179841 - 3401272075258036224000*sqrt(3)*sqrt(-154227 + 89043*sqrt(3))/52283326179841 - 18210720*sqrt(3)/7230721 + 47757270/7230721 + 5891187146667514675200*sqrt(3)*sqrt(-51409 + 29681*sqrt(3))/52283326179841 + 10203835454942977476480*sqrt(-51409 + 29681*sqrt(3))/52283326179841, 10*sqrt(-185511737623 - 132968544*sqrt(-51409 + 29681*sqrt(3)) + 76769280*sqrt(3)*sqrt(-51409 + 29681*sqrt(3)) + 107105251656*sqrt(3))/sqrt(3836935845121 - 2215255943040*sqrt(3)))

   r1 = 2.1539;  x1 = -9.28203
   bx = 5.48676;  by = 1.63965;  r2 = 4.75039

菱形の中の円の半径は 4.75038687805787 である。

res[1][4] |> simplify

res[1][4].evalf()

   4.75038687805787

術では以下の式が示されているが,計算される値は不適切なもののようだ。

47/(16*sqrt(57√3 - 73) - 141(7 - 4√3)) * 10

   6.616796128691446

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 10
   (r1, x1, bx, r2) = res[1]
   by = a - sqrt(a^2 - bx^2)
   @printf("r1 = %g;  x1 = %g\n", r1, x1)
   @printf("bx = %g;  by = %g;  r2 = %g\n", bx, by, r2)
   plot()
   circle(0, a, a, :blue)
   circle(x1, r1, r1)
   circle(x1 - 2r1, r1, r1)
   circle(x1 - r1, (1 + √3)r1, r1)
   circle(0, r1, r1)
   circle(0, 2a - r1, r1)
   circle(0, a, r2, :green)
   segment(a, a, -bx, by)
   segment(-a, a, bx, by)
   segment(a, a, -bx, 2a-by)
   segment(-a, a, bx, 2a-by)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, r1, "(x1,r1)", :red, :center, delta=-delta)
       point(x1 - r1, r1*(1 + √3), "(x1-r1,(1+√3)r1)", :red, :center, delta=-delta)
       point(0, a, " a", :green)
       point(float(bx), float(by), " (bx,by)", :black)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その413)

2023年09月02日 | Julia

算額(その413)

石川県七尾市中島町宮前(旧熊木村) 久麻加夫都阿良加志比古神社 文政6年(1823)

早坂平四郎: 算額の一考察,苫小牧工業専門学校紀要
https://www.tomakomai-ct.ac.jp/wp01/wp-content/uploads/2014/06/kiyou5-8.pdf

半径 a の外円の中に,交わる 2 個の大円と,半径の等しい 4 個の小円が入っている。a を知って,大円と小円の半径を求めよ。

外円の半径と中心座標を a, (0, 0)
大円の半径と中心座標を r1, (a - r1, 0)
小円の半径と中心座標を r2, (0, a-r2), (r2, 0)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, r1::positive, r2::positive;
eq1 = a - 2r1 + 2r2
eq2 = (a - r1)^2 + (a - r2)^2 - (r1 + r2)^2
res = solve([eq1, eq2], (r1, r2))

   1-element Vector{Tuple{Sym, Sym}}:
    (a*(-5/4 + sqrt(33)/4) + a/2, a*(-5/4 + sqrt(33)/4))

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

   a*(-3 + sqrt(33))/4
   a*(-5 + sqrt(33))/4

大円の半径は (√33 - 3)a/4,小円の半径は (√33 - 5)a/4 である。

a = 10 のとき,6.8614066163450715 と 1.8614066163450715。

a = 10
a*(-3 + sqrt(33))/4, a*(-5 + sqrt(33))/4

   (6.8614066163450715, 1.8614066163450715)

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 10
   r1 = a*(-3 + sqrt(33))/4
   r2 = a*(-5 + sqrt(33))/4
   @printf("r1 = %g;  r2 = %g\n", r1, r2)
   plot()
   circle(0, 0, a, :black)
   circle(a - r1, 0, r1)
   circle(r1 - a, 0, r1)
   circle(0, a - r2, r2, :blue)
   circle(0, r2 - a, r2, :blue)
   circle(r2, 0, r2, :blue)
   circle(-r2, 0, r2, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, a - r2, " a-r2", :blue)
       point(r2, 0, "r2", :blue, delta=-delta)
       point(a - r1, 0, "a-r1", :red, :right, :bottom, delta=delta)
       point(a, 0, "a ", :black, :right, :bottom, delta=delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その412)

2023年09月02日 | Julia

算額(その412)

百十五 群馬県富岡市一ノ宮 貫前神社 明治4年(1871)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.

早坂平四郎: 算額の一考察
https://www.tomakomai-ct.ac.jp/wp01/wp-content/uploads/2014/06/kiyou5-8.pdf

キーワード:円6個

図のように,相接する等円 3 個はその中心が直線上にある。中央の等円の上に大円が乗り,大円と等円に外接する小円が 2 個ある。
大円と小円の半径がわかっているとき,等円の半径を求める式を作れ。

大円の半径と中心座標を r2, (0, r1 + r2)
小円の半径と中心座標を r3, (x3, y3) ただし x3 = r1
等円の半径と中心座標を r1, (0, 0), (2r1, 0)
とおき以下の連立方程式を解く。

「大円と小円の半径がわかっているとき,等円の半径を求める式」は,sympy では簡略化できない長い式なので,大円と小円の半径を与えて等円の半径を求める。」

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive,
      r3::positive, x3::positive, y3::positive;

(r2, r3) = (5, 1)
x3 = r1
eq1 = (2r1 - x3)^2 + y3^2 - (r1 + r3)^2
eq2 = x3^2 + (r1 + r2 - y3)^2 - (r2 + r3)^2
res = solve([eq1, eq2], (r1, y3))

   3-element Vector{Tuple{Sym, Sym}}:
    (((-1/2 - sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3) - 10/(3*(-1/2 - sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3)))*(2 + (-1/2 - sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3) - 10/(3*(-1/2 - sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3)))/2, 1 + (-1/2 - sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3) - 10/(3*(-1/2 - sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3)))
    ((-10/(3*(-1/2 + sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3)) + (-1/2 + sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3))*(2 - 10/(3*(-1/2 + sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3)) + (-1/2 + sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3))/2, 1 - 10/(3*(-1/2 + sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3)) + (-1/2 + sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3))
    ((-10/(3*(10 + 10*sqrt(111)/9)^(1/3)) + (10 + 10*sqrt(111)/9)^(1/3))*(-10/(3*(10 + 10*sqrt(111)/9)^(1/3)) + 2 + (10 + 10*sqrt(111)/9)^(1/3))/2, -10/(3*(10 + 10*sqrt(111)/9)^(1/3)) + 1 + (10 + 10*sqrt(111)/9)^(1/3))

3 番目の組が適解である。
大円,小円の半径を 5 寸, 1 寸としたとき,
等円の半径は 2.865876288528763 寸である。

(-10/(3*(10 + 10*sqrt(111)/9)^(1/3)) + (10 + 10*sqrt(111)/9)^(1/3))*(-10/(3*(10 + 10*sqrt(111)/9)^(1/3)) + 2 + (10 + 10*sqrt(111)/9)^(1/3))/2

   2.865876288528763

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, y3) = res[3]
   r2 = 5
   r3 = 1
   x3 = r1
   @printf("r1 = %g;  x3 = %g;  y3 = %g\n", r1, x3, y3)
   plot()
   circle(0, 0, r1)
   circle(2r1, 0, r1)
   circle(-2r1, 0, r1)
   circle(x3, y3, r3, :blue)
   circle(-x3, y3, r3, :blue)
   circle(0, r1 + r2, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r1 + r2, " r1+r2  大円")
       point(r1, 0, " r1", :red, delta=-delta)
       point(2r1, 0, " 2r1 等円", :red, delta=-delta)
       point(x3, y3, "小円:r3,(x3,y3)", :blue, :left, :bottom, delta=2delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

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

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