裏 RjpWiki

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

算額(その920)

2024年05月06日 | ブログラミング

算額(その920)

香川県善通寺市与北町 皇美屋神社 明治11年(1878)
本田益夫:算額随想-香川県内の算額について-,私家版,昭和45年(1970).

直線上に大円 2 個,小円 3 個,内円 1 個が積み上がっている。大円,小円の直径がそれぞれ 14 寸,7 寸のとき,内円の直径はいかほどか。

大円の半径と中心座標を r1, (r1, r1)
小円の半径と中心座標を r2, (r2, y2), (0, y22)
内円の半径と中心座標を r3, (0, y3)
とおき,以下の連立方程式を解く。

大円,小円の直径が特定の値のときには,直接数値を代入して連立方程式を解けばよい。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, y2::positive, y22::positive,
     r3::positive, y3::positive
(r1, r2) = (14, 7) .// 2
eq1 = (r1 - r2)^2 + (y2 - r1)^2 - (r2 + r1)^2
eq2 = r1^2 + (y3 - r1)^2 - (r3 + r1)^2
eq3 = r2^2 + (y2 - y3)^2 - (r2 + r3)^2
eq4 = r2^2 + (y22 - y2)^2 - 4r2^2;
res = solve([eq1, eq2, eq3, eq4], (y2, y22, r3, y3))[2]  # 2 of 2

   (7 + 7*sqrt(2), 7*sqrt(3)/2 + 7 + 7*sqrt(2), 2, 4*sqrt(2) + 7)

大円,小円の直径がそれぞれ 14 寸,7 寸のとき,内円の直径は 4 寸である。

---

以下では,一般解を求める。

@syms r1::positive, r2::positive, y2::positive, y22::positive,
     r3::positive, y3::positive
eq1 = (r1 - r2)^2 + (y2 - r1)^2 - (r2 + r1)^2 |> expand
eq2 = r1^2 + (y3 - r1)^2 - (r3 + r1)^2 |> expand
eq3 = r2^2 + (y2 - y3)^2 - (r2 + r3)^2 |> expand
eq4 = r2^2 + (y22 - y2)^2 - 4r2^2 |> expand;

y2, y3 は簡単に求めることができる。

println(eq1, " から y2 を求める")
ans_y2 = solve(eq1, y2)[2]  # 2 of 2
ans_y2 |> println

   r1^2 - 4*r1*r2 - 2*r1*y2 + y2^2 から y2 を求める
   2*sqrt(r1)*sqrt(r2) + r1

println(eq2, " から y3 を求める")
ans_y3 = solve(eq2, y3)[2]  # 2 of 2
ans_y3 |> println

   r1^2 - 2*r1*r3 - 2*r1*y3 - r3^2 + y3^2 から y3 を求める
   r1 + sqrt(r3)*sqrt(2*r1 + r3)

r3 は,上で求めた y2, y3 を eq3 に代入し,方程式を解けばよい。もう少し簡約化することもできるが,大差ない。

eq13 = eq3(y2 => ans_y2, y3 => ans_y3) |> expand
println(eq13, " から r3 を求める")
ans_r3 = solve(eq13, r3)[1]
ans_r3 |> println

   -4*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(2*r1 + r3) + 4*r1*r2 + 2*r1*r3 - 2*r2*r3 から r3 を求める
   2*(-2*sqrt(2)*r1^(3/2)*r2^(3/2) + r1*r2*(r1 + r2))/(r1^2 - 6*r1*r2 + r2^2)

y22 は eq4 を解き,式中の y2 に,すでに求めた ans_y2 を代入すればよい。

println(eq4, " から y22 を求める")
ans_y22 = solve(eq4, y22)[2]
ans_y22(y2 => ans_y2) |> simplify |> println

   -3*r2^2 + y2^2 - 2*y2*y22 + y22^2 から y22 を求める
   2*sqrt(r1)*sqrt(r2) + r1 + sqrt(3)*r2

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (14, 7) .// 2
   r3 = 2*(-2√(2(r1^3*r2^3)) + r1*r2*(r1 + r2))/(r1^2 - 6*r1*r2 + r2^2)
   y3 = r1 + √(2r1*r3 + r3^2)
   y2 = r1 + 2√(r1*r2)
   y22 = 2√(r1*r2) + r1 + √3r2
   plot()
   circle2(r1, r1, r1)
   circle2(r2, y2, r2, :blue)
   circle(0, y22, r2, :blue)
   circle(0, y3, r3, :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(r1, r1, "大円:r1,(r1,r1)", :red, :center, delta=-delta/2)
       point(r2, y2, "小円:r2,(r2,y2)", :blue, :center, delta=-delta/2)
       point(0, y22, "小円:r2,(0,y22)", :blue, :center, delta=-delta/2)
       point(0, y3, "内円:r3,(0,y3)", :black, :center, delta=-delta/2)
       segment(-2r1, 0, 2r1, 0, :magenta, lw=0.5)
   end
end;

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

コメントを投稿

ブログラミング」カテゴリの最新記事