裏 RjpWiki

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

算額(その519)

2023年11月28日 | Julia

算額(その519)

岩手県一関市 観福寺 明治34年
一関市博物館>>和算に挑戦>>平成22年度出題問題&解答例>>平成22年度出題問題(2)[中級問題]

https://www.city.ichinoseki.iwate.jp/museum/wasan/h22/normal.html

外円内に n 個の菱形が入っている。外円の直径を r としたとき,菱形の長い方の対角線の長さを求めよ。

菱形の短い方の対角線の長さを,長い方の対角線の長さを 2b,2a とおき,以下の方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b, r, n
b = r/n
eq = sqrt(a^2 + (r - b)^2) - r
a = solve(eq, a)[2]
a |> simplify |> println

   sqrt(r^2*(2*n - 1))/n

菱形の長い方の対角線の長さは 2r*sqrt(2n - 1)/n。

外円の直径が 10,n = 2 のとき,2r*sqrt(2n - 1)/n = 8.660254037844386 である。

r = 10/2
n = 2
2r*sqrt(2n - 1)/n

   8.660254037844386

draw(5, 2, true)
   a = 4.33013;  b = 2.5;  長い方の対角線 = 8.66025

draw(5, 5, true)
   a = 3;  b = 1;  長い方の対角線 = 6

function draw(r, n, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = r//n
   a = r*sqrt(2n - 1)/n
   @printf("a = %g;  b = %g;  長い方の対角線 = %g\n", a, b, 2a)
   plot()
   circle(0, 0, r)
   for i = 1:n
       plot!([a, 0, -a, 0, a], r - (2i - 1)b .+ [0, b, 0, -b, 0], color=:blue, lw=0.5)
       point(0, r - (2i - 1)b, "r-$(2i-1)b ", :blue, :right, :vcenter)
   end
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(r*sqrt(2*n - 1)/n, r - b, "(r√(2n-1)/n,r-b)", :black, :right, :top, delta=-delta)
       point(r, 0, " r", :red, :left, :bottom, delta=delta/2)
   end
end;

 

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

算額(その518)

2023年11月28日 | Julia

算額(その518)

村松茂清: 『算俎』,寛文3年(1663)
一関市博物館>>和算に挑戦>>平成23年度出題問題&解答例>>平成23年度出題問題(2)[中級問題]

https://www.city.ichinoseki.iwate.jp/museum/wasan/h23/normal.html

正三角形内に等円 3 個が内接している。等円の直径が 10 のとき,黒積 A, B の面積を求めよ。

等円の半径が r,正三角形の一辺の長さが 2a とする。

include("julia-source.txt");

using SymPy
@syms a, r
eq = r - (a - r)/sqrt(Sym(3))
a = solve(eq, a)[1]
a |> println

   r*(1 + sqrt(3))

正三角形の一辺の長さは等円の半径の 1 + √3 倍である。

   r = 5;  a = 13.6603;  等円の直径 = 10;  正三角形の一辺の長さ = 27.3205

黒積 A は,「底辺 (a - r),高さ r の直角三角形から,半径 r の円の面積の 1/6 を引いたもの」の 2 倍である。
-r^2*(pi - 3*sqrt(3))/3 = (√3 - π/3)r^2

2((a - r)r/2 - PI*r^2/6) |> factor |> println

   -r^2*(pi - 3*sqrt(3))/3

黒積 B は,3個の等円の中心を結んでできる正三角形の面積から,半径が r の円の面積の半分(円/6 の 3 個分)を引いたものである。
-r^2*(pi - 2*sqrt(3))/2 = (√3 - π/2)r^2

2r * √Sym(3)r/2 - PI*r^2/2 |> factor |> println

   -r^2*(pi - 2*sqrt(3))/2

function black(r, a)
   x = vcat([cosd(θ)r for θ in 240:300], [cosd(θ)r + r for θ in 120:180], [cosd(θ)r - r for θ in 0:60])
   y = vcat([sind(θ)r + (1 + √3)r for θ in 240:300], [sind(θ)r + r for θ in 120:180], [sind(θ)r + r for θ in 0:60])
   plot!(x, y, color=:black, lw=0.5, seriestype=:shape, fillcolor=:black, alpha=0.5)
   x = vcat([cosd(θ)r + r for θ in 270:390], [a, r])
   y = vcat([sind(θ)r + r for θ in 270:390], [0, 0])
   plot!(x, y, color=:black, lw=0.5, seriestype=:shape, fillcolor=:black, alpha=0.5)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 10//2
   a = (1 + √3)r
   @printf("r = %g;  a = %g;  等円の直径 = %g;  正三角形の一辺の長さ = %g\n", r, a, 2r, 2a)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:blue, lw=0.5)
   circle(r, r, r)
   circle(-r, r, r)
   circle(0, (1 + √3)r, r)
   black(r, a)
   plot!([r, 0, -r, r], r .+ [0, √3r, 0, 0], color=:green, lw=0.5)
   plot!([r, r, a, r], [0, r, 0, 0], color=:magenta, lw=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0.15r, 1.5r, "B", :black, mark=false)
       point(1.6r, 0.5r, "A", :black, mark=false)
       point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(0, √3a, " √3a", :blue, :left, :vcenter)
       point(0, (1 + √3)r, " (1+√3)r", :red, :left, :vcenter)
       point(r, r, " (r,r)", :red, :left, :vcenter)
   end
end;

 

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

算額(その517)

2023年11月28日 | Julia

算額(その517)

千葉胤秀: 『算法新書』巻の五,文政13年(1830)
一関市博物館>>和算に挑戦>>平成23年度出題問題&解答例>>平成23年度出題問題(3)[上級問題]

https://www.city.ichinoseki.iwate.jp/museum/wasan/h23/hard.html

外円内に5個の円が内接・外接している。外円の直径は 6 寸,甲円の直径が 2 寸のとき,丙円の直径を求めよ。

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

include("julia-source.txt");

using SymPy

@syms r0, r1, r2, x2, y2, r3, r4

eq1 = x2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2
eq2 = x2^2 + (y2 - r3 + r0)^2 - (r2 + r3)^2
eq3 = x2^2 + y2^2 - (r0 - r2)^2
eq4 = x2^2 + (r0 - 2r1 - r4 - y2)^2 - (r2 + r4)^2
eq5 = r1 + r3 + r4 - r0

res = solve([eq1, eq2, eq3, eq4, eq5], (r2, x2, y2, r3, r4))

   3-element Vector{NTuple{5, Sym}}:
    (0, 0, r0, r0, -r1)
    (2*r0*r1*(r0 - r1)/(r0^2 + r1^2), -2*sqrt(2)*r0*r1*(r0 - r1)/(r0^2 + r1^2), r0*(r0^2 - 2*r0*r1 - r1^2)/(r0^2 + r1^2), r0*(r0 - r1)/(r0 + r1), (r0*r1 - r1^2)/(r0 + r1))
    (2*r0*r1*(r0 - r1)/(r0^2 + r1^2), 2*sqrt(2)*r0*r1*(r0 - r1)/(r0^2 + r1^2), r0*(r0^2 - 2*r0*r1 - r1^2)/(r0^2 + r1^2), r0*(r0 - r1)/(r0 + r1), (r0*r1 - r1^2)/(r0 + r1))

3 組の解が得られるが,3 番目の解が適解である。

乙円の半径(r2)は r0*(r0 - r1)/(r0 + r1) で,外円の直径(2r0)が 6 寸,甲円の直径(2r1)が 2 寸のとき,丙円の直径は 3.0 である。

r0 = 3
r1 = 1
2(r0*(r0 - r1)/(r0 + r1))

   3.0

   r0 = 3;  r1 = 1;  r2 = 1.2;  x2 = 1.69706;  y2 = 0.6;  r3 = 1.5;  r4 = 0.5

「甲円の直径(半径)の逆数」と「丙円の直径(半径)の逆数」の和は,「丙円の逆数の2倍」と等しい。
今回の図では左右にある丙円は同じ大きさであるが大きさが異なる場合も「丙円1の逆数と丙円2の逆数の和」に等しい。

1/r1 + 1/res[3][4] |> simplify |> println

   (1.0*r0*(r0 - r1) + r0 + r1)/(r0*(r0 - r1))

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

   (r0^2 + r1^2)/(r0*r1*(r0 - r1))

1/r1 + 1/r3 = 1.66667;  2/r2 = 1.66667

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

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

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

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