裏 RjpWiki

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

算額(その305)

2023年06月29日 | Julia

算額(その305)

「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県伊賀市 菅原神社 嘉永7年(1854)

問題文4 

外大円の中に長方形を起き,菱形を切り取った残りの領域に円を置く。円の直径が 3 尺,長方形の残りの長さ BE が 4.5 尺である。このとき,菱形の一辺(菱面=BD),長軸(同立=AD),短軸(同横= BC),外大円の円周,矢(r0 - y)を求めよ。

外大円の半径を r0 とおく。r0 = 同立/2
長方形の右上の頂点の座標を (x, y) とする。
以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, x::positive, y::positive,
   菱面::positive, 同立::positive, 同横::positive,
   外大円円周::positive, 矢::positive;

r1 = 3/2
eq1 = 4.5 + 2y - sqrt(4.5^2 + (2y)^2) - 2r1
eq2 = 4.5^2 + (2y)^2 - 菱面^2  # ⊿BDE において
eq3 = (菱面 + 4.5)^2 + (2y)^2 - 同立^2  # ⊿ADE において
eq4 = (同立/2)^2 + (同横/2)^2 - 菱面^2  # ⊿OBE において
eq5 = 外大円円周 - 同立*PI
eq6 = 矢 - (同立/2 - y);

eq1 〜 eq6 までの6元連立方程式にすると解けなくなってしまうので,eq1 〜 eq5 の5元連立方程式を解く。

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (y, 菱面, 同立, 同横, 外大円円周, 矢))

res = solve([eq1, eq2, eq3, eq4, eq5], (y, 菱面, 同立, 同横, 外大円円周))

   1-element Vector{NTuple{5, Sym}}:
    (3.00000000000000, 7.50000000000000, 13.4164078649987, 6.70820393249937, 42.1488883862444)

res[1][3]/2 - res[1][1] |> println  # 矢

   3.70820393249937

なお,solve(eq1) で y だけを求めれば,ドミノ倒しで残りの変数は確定する。

res = solve(eq1)[1] |> println  # y

   3.00000000000000

sqrt(6^2 + 4.5^2)  # 菱面 = sqrt((2y)^2 + 4.5^2)

   7.5

sqrt((7.5 + 4.5)^2 + (2*3)^2)  # 同立 = sqrt((菱面 + 4.5) + (2y)^2)

   13.416407864998739

sqrt(7.5^2 - (13.416407864998739/2)^2)*2  # 同横 = sqrt(菱面^2-(同立/2)^2)*2

   6.708203932499367

13.4164078649987*pi  # 外大円円周 = 同立*π

   42.14888838624424

(13.416407864998739/2) - 3  # 矢 = 同立/2 - y

   3.7082039324993694

   y = 3.000000
   菱面 = 7.500000
   同立 = 13.416408
   同横 = = 6.708204
   外大円円周 = = 42.148888
   矢 = 3.708204

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (y, 菱面, 同立, 同横, 外大円円周) = (3.00000000000000, 7.50000000000000, 13.4164078649987, 6.70820393249937, 42.1488883862444)
   矢 = 3.70820393249937
   @printf("y = %.6f\n", y)
   @printf("菱面 = %.6f\n", 菱面)
   @printf("同立 = %.6f\n", 同立)
   @printf("同横 = = %.6f\n", 同横)
   @printf("外大円円周 = = %.6f\n", 外大円円周)
   @printf("矢 = %.6f\n", 矢)
   plot([6, 6, -6, -6, 6], [-y, y, y, -y, -y], color=:black, lw=0.5)
   plot!([6, 1.5, -6, -1.5, 6], [-y, y, y, -y, -y], color=:red, lw=0.5)
   circle(4.5, 1.5, 1.5, :green)
   segment(1.5, y, -1.5, -y, :blue)
   segment(-6, y, 6, -y, :magenta)
   if more
       point(4.5, 1.5, "(4.5,1.5)", :green, :center)
       point(1.5, y, "B:(1.5,y) ", :red, :left, :bottom)
       point(6, y, "E:(6,y)", :black, :left, :bottom)
       point(-6, y, " A", :green, :left, :bottom)
       point(-1.5, -y, " C", :green, :left, :bottom)
       point(6, -y, " D", :green, :left, :bottom)
       point(0, 0, " O", :green, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その304)

2023年06月29日 | Julia

算額(その304)

「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県伊賀市 菅原神社 嘉永7年(1854)

問題文1 

正方形に内接する円を描き,円の中に図のように菱形と2種類の円を2つずつ入れる。菱形の一辺と短軸は等しい。小円の直径が 3.55 尺のとき,正方形の一変,正方形に内接する円の円周,菱形の一変,菱形内の大円の直径,矢(菱形の上の頂点から正方形の上辺までの距離)を求めよ。

正方形の中心を原点とする。正方形の一辺の長さを 2r0 とする(内接する円の半径が r0)。菱形内の小円の半径,中心座標を r1, (r1, y1),大円の半径,中心座標を r2, (r2, 0) とする。
以下の連立方程式をとき,y1, r2, r0 を求める。

include("julia-source.txt");

using SymPy
@syms r1::positive, y1::positive, r2::positive, r0::positive, y::positive;

r1 = 3.55/2
sqrt3 = sqrt(Sym(3))
y = r0/sqrt3
y = r2*sqrt3
eq1 = (r2 - r1)^2 + y1^2 - (r2 + r1)^2
eq2 = y - y1 - sqrt3*r1
eq3 = 2y - sqrt(y^2 + r0^2);

res = solve([eq1, eq2, eq3], (y1, r2, r0))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (2*sqrt(3)*r1, 3*r1, 9*r1)

y1 = 2*sqrt(3)*r1
r2 = 3r1
r0 = 9r1
である。他の変数の値は,これらから求められる。

   中円の径 = 2r2 = 10.650000
   菱面 = 2y = 18.446341
   方面 = 2r0 = 31.950000
   円廻り = π*2r0 = 100.373885
   矢 = r0 - y = 6.751829

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 3.55/2
   y1 = 2*sqrt(3)*r1
   r2 = 3r1
   r0 = 9r1
   y = r0/sqrt(3)
   @printf("中円の径 = 2r2 = %.6f\n", 2r2)
   @printf("菱面 = 2y = %.6f\n", 2y)
   @printf("方面 = 2r0 = %.6f\n", 2r0)
   @printf("円廻り = π*2r0 = %.6f\n", pi*2r0)
   @printf("矢 = r0 - y = %.6f\n", r0 - y)
   plot([r0, r0, -r0, -r0, r0], [-r0, r0, r0, -r0, -r0], color=:black, lw=0.5)
   circle(0, 0, r0, :green)
   circle(r1, y1, r1, :blue)
   circle(-r1, y1, r1, :blue)
   circle(r2, 0, r2, :red)
   circle(-r2, 0, r2, :red)
   plot!([r0, 0, -r0, 0, r0], [0, y, 0, -y, 0], color=:magenta, lw=0.5)
   x = sqrt(r0^2 - y^2)
   segment(-x, y, x, y)
   if more
       point(0, y, " y", :magenta, :left, :bottom)
       point(0, r0, " r0", :green, :left, :bottom)
       point(r0, 0, " r0", :green, :left, :bottom)
       point(r2, 0, " 中円\n r2", :red, :left, :bottom)
       point(r1, y1, " 小円(r1,y1)", :blue, :left, :bottom)
       point(0, (r0 + y)/2, " 矢=r0-y", :black, mark=false)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その303)

2023年06月29日 | Julia

算額(その303)

「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県伊賀市 菅原神社 嘉永7年(1854)

問題文3 

外側の正三角形の一辺の長さを 2a,一番大きい正方形,二番目,三番目の正方形の一辺の長さをそれぞれ 2b, 2c, 2d とし,以下の連立方程式を解く。2d = 7.8179 が与えられているが,d も未知数として方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, c::positive, d ::positive

d = 7.8179/2
sqrt3 = sqrt(Sym(3))
eq1 = sqrt3*a - 2(b + c + d) - sqrt3*d
eq2 = sqrt3*a - 2(b + c) - sqrt3*c
eq3 = sqrt3*a - 2b - sqrt3*b

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

   Dict{Any, Any} with 3 entries:
     b => 4*sqrt(3)*d/3 + 7*d/3
     a => 5*d + 26*sqrt(3)*d/9
     c => d + 2*sqrt(3)*d/3

res[a] |> simplify |> println
res[b] |> simplify |> println
res[c] |> simplify |> println

   d*(45 + 26*sqrt(3))/9
   d*(4*sqrt(3) + 7)/3
   d*(3 + 2*sqrt(3))/3

二方(2 番目に大きい正方形の一辺の長さ)= 2c

res[c](d => 7.8179/2).evalf()*2

   16.8452333389952

一方(1 番目大きい正方形の一辺の長さ)= 2b

res[b](d => 7.8179/2).evalf()*2

   36.296433344657

外三角(外側の正三角形の一辺の長さ)

res[a](d => 7.8179/2).evalf()*2

   78.207944468979

三角面(内側の正三角形の一辺の長さ)
内側の正三角形は正方形に内接する円(半径 b)に内接する

(res[b](d => 7.8179/2) * sqrt3/2 * 2).evalf()

   31.4336333432415

小円径(内側の正三角形に内接する円の直径)b/2 * 2 である

(res[b](d => 7.8179/2) / 2 * 2).evalf()

   18.1482166723285

   a = 39.103972;  b = 18.148217;  c = 8.422617; d = 3.908950
   二方 = 2c = 16.845233
   一方 = 2b = 36.296433
   外三角 = 2a = 78.207944
   三角面 = √3b = 31.433633
   小円径 = b = 18.148217

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   d = 7.8179/2
   a = d*(45 + 26*sqrt(3))/9
   b = d*(4*sqrt(3) + 7)/3
   c = d*(3 + 2*sqrt(3))/3
   @printf("a = %.6f;  b = %.6f;  c = %6f; d = %.6f\n", a, b, c, d)
   @printf("二方 = 2c = %.6f\n", 2c)
   @printf("一方 = 2b = %.6f\n", 2b)
   @printf("外三角 = 2a = %.6f\n", 2a)
   @printf("三角面 = √3b = %.6f\n", √3b)
   @printf("小円径 = b = %.6f\n", b)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:black, lw=0.5)
   rect(-d, 2(b + c), d, 2(b + c + d), :red)
   rect(-c, 2b, c, 2(b + c), :blue)
   rect(-b, 0, b, 2b, :green)
   circle(0, b, b, :brown)
   circle(0, b, b/2, :orange)
   plot!([b√3/2, 0, -b√3/2, b√3/2], [b/2, 2b, b/2, b/2], color=:magenta, lw=0.5)
   if more
       point(b, 2b, " (b,2b)")
       point(c, 2(b + c), " (c,2b+2c)", :blue)
       point(d, 2(b + c + d), " (c,2b+2c+2d)", :red)
       point(a, 0, " a", :black, :left, :bottom)
       point(0, √3a, " √3a", :brown, :left, :bottom)
       point(0, b, " b", :brown)
       point(0, b/2, " b/2", :orange)
       point(√3b/2, b/2, "  (√3b/2,b/2)", :magenta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その302)

2023年06月29日 | Julia

算額(その302)

「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
三重県伊賀市 菅原神社 嘉永7年(1854)

問題文5 面積が 163350 歩平方の鈎股弦がある。鈎股弦の長さの差は分からないが,弦が 82.5 尺であるとき,鈎,股,中勾,方面,小面を求めよ。
注: 1 尺 = 10 歩

正方形の一辺の長さを a として,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive, a ::positive

eq1 = 鈎 * 股 / 2 - 1633.50
eq2 = 鈎^2 + 股^2 - 82.5^2
eq3 = (鈎 - a)/a - a/(股 - a)

solve([eq1, eq2, eq3], (鈎, 股, a))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (49.5000000000000, 66.0000000000000, 28.2857142857143)
    (66.0000000000000, 49.5000000000000, 28.2857142857143)

鈎 < 股 ゆえ,最初の組が適解。

鈎 = 49.500000
股 = 66.000000
中勾(正方形の対角線) = 40.002041
方面(正方形の一辺の長さ) = 28.285714
小面(小さな正方形の一辺の長さ) = 20.001020

using Plots

function draw(zoomin=false, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, a) = (49.5000000000000, 66.0000000000000, 28.2857142857143)
   @printf("鈎 = %.6f\n", 鈎)
   @printf("股 = %.6f\n", 股)
   @printf("中勾(正方形の対角線) = %.6f\n", √2a)
   @printf("方面(正方形の一辺の長さ) = %.6f\n", a)
   @printf("小面(小さな正方形の一辺の長さ) = %.6f\n", a/√2)
   plot([0, 股 , 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   rect(0, 0, a, a, :red)
   plot!([a/2, a, a/2, 0, a/2], [0, a/2, a, a/2, 0], color=:blue, lw=0.5)
   if more
       point(0, 鈎, " 鈎", :black, :left, :bottom)
       point(股, 0, " 股", :black, :left, :bottom)
       point(a, 0, " a", :red, :left, :bottom)
       point(a, a/2, " (a,a/2)", :blue, :left, :bottom)
       point(a/2, a, " (a/2,a)", :blue, :left, :bottom)
       point(a, a, " (a,a)", :red, :left, :bottom)
       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アクセスランキング にほんブログ村