裏 RjpWiki

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

算額(その301)

2023年06月28日 | Julia

算額(その301)

「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県伊賀市 永保寺 天保15年(1844)

正三角形の中に図のように円,正方形が入っている。正三角形の面積が 2625 平方寸のとき,それぞれの大きさを求めよ。

正三角形の一辺の長さを 2a とする。
下の同心円の中心座標は (0, a/2),大きい円の半径 r1 は a/2 である。
上にある円の半径を r2 とし,以下の方程式を解き,a, r1, r2 を求める。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive
eq1 = a*sqrt(Sym(3))a - 2625
eq2 = a*tan(PI/6) - r1
eq3 = (sqrt(Sym(3))a - r2 - 2r1)sin(PI/6) - r2;

今回は,正三角形の面積が 2625 平方寸ということで,a が特定されている。

res = solve(eq1, a)
res[1] |> println

   5*3^(1/4)*sqrt(35)  # a

一般解を求める場合には以下の2元連立方程式を解く。

res = solve([eq2, eq3], (r1, r2))

   Dict{Any, Any} with 2 entries:
     r2 => sqrt(3)*a/9
     r1 => sqrt(3)*a/3

ついで,図中の b 及び小円の半径 r3 を求める。

   a = 38.929994;  r1 = 22.476243;  r2 = 7.492081
   b = 15.893104;  r3 = 11.238121
   三角面(正三角形の一辺の長さ)= 77.85998861
   中勾(三角形の高さ)= 67.42872808
   下円径(大きい円の直径)= 44.95248538
   上円径(上の円の直径)= 14.98416179
   大方面(大きい正方形の一辺の長さ) = 31.78620725
   小方面(小さい正方形の一辺の長さ) = 22.47624269
   小円径(小円の直径) = 22.47624269
   小円廻り(小円の円周) = 70.61119892

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 5*3^(1/4)*sqrt(35)
   r1 = sqrt(3)*a/3
   r2 = sqrt(3)*a/9
   b = r1*cos(PI/4)
   r3 = b*cos(PI/4)
   @printf("a = %.6f;  r1 = %.6f;  r2 = %.6f\nb = %.6f;  r3 = %.6f\n", a, r1, r2, b, r3)
   @printf("三角面(正三角形の一辺の長さ)= %.8f\n", 2a)  # 77.85998861
   @printf("中句(三角形の高さ)= %.8f\n", a * √3)  # 67.42872808
   @printf("下円径(大きい円の直径)= %.8f\n", sqrt(3)*a/3 * 2)  # 44.95248538
   @printf("上円径(上の円の直径)= %.8f\n", sqrt(3)*a/9 * 2)  # 14.98416179
   @printf("大方面(大きい正方形の一辺の長さ) = %.8f\n", 2r1*cos(pi/4))  # 31.78620725
   @printf("小方面(小さい正方形の一辺の長さ) = %.8f\n", 2(b/√2))  # 22.47624269
   @printf("小円径(小円の直径) = %.8f\n", 2r3)  # 70.61119892087622
   @printf("小円廻り(小円の円周) = %.8f\n", 2r3*pi)  # 70.61119892087622
   plot([a, 0 , -a, a], [0, √3a, 0, 0], color=:black, lw=0.5)
   circle(0, r1, r1)
   circle(0, r2 + 2r1, r2, :blue)
   rect(-b, r1 - b, b, r1 + b, :magenta)
   plot!([b, 0, -b, 0, b], [r1, r1 + b, r1, r1 - b, r1], color=:green, lw=0.5)
   circle(0, r1, r3, :orange)
   if more
       point(0, √3a, " √3a", :black, :left, :bottom)
       point(a, 0, " a", :black, :left, :bottom)
       point(0, r1, " r1", :red, :left, :bottom)
       point(0, r1+r3, "\n r1+r3", :orange, :left, :top)
       point(0, r1 + b, " r1+b", :green, :left, :bottom)
       point(0, r1 - b, " r1-b", :green, :left, :top)
       point(b, r1, "\n(b,r1)", :green, :left, :top)
       point(0, 2r1 + r2, " 2r1+r2", :blue)
       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アクセスランキング にほんブログ村