裏 RjpWiki

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

算額(その523)

2023年11月30日 | Julia

算額(その523)

和算図形問題あれこれ - 令和4年11月の問題-No.1
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html

2 個の甲円が交差しており,乙円 2 個と丙円 3 個が入っている。甲円,丙円の直径が 100 寸,36 寸であるとき,乙円の直径はいくらか。

甲円の半径と中心座標を r1, (0, r1 - r3), (0, r3 - r1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (x3, y3), (0, 0)
とおき,以下の連立方程式を解く。
なお,r1, r3 を変数のままとして解くことは,solve() の性能上無理なようで,問に書かれた通りの数値を代入して解を求める。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, x2::negative, y2::positive,
   r3::positive, x3::positive, y3::positive;
(r1, r3) = (100, 36) .// 2
eq1 = x2^2 + (y2 - (r3 - r1))^2 - (r1 + r2)^2
eq2 = x3^2 + (y3 - (r3 - r1))^2 - (r1 + r3)^2
eq3 = x2^2 + (y2 - (r1 - r3))^2 - (r1 - r2)^2
eq4 = x3^2 + (y3 - (r1 - r3))^2 - (r1 - r3)^2
eq5 = (x3 - x2)^2 + (y3 - y2)^2 - (r2 + r3)^2;
res = solve([eq1, eq2, eq3, eq4, eq5], (r2, x2, y2, x3, y3))

   1-element Vector{NTuple{5, Sym}}:
    (20500*sqrt(7)/4181 + 72242/4181, -66075*sqrt(41)/16724 + 25215*sqrt(287)/33448, 128125*sqrt(7)/16724 + 903025/33448, 15*sqrt(287)/8, 225/8)

乙円の半径 r2 は 20500√7/4181 + 72242/4181 = 82*(250√7 + 881)/4181,直径は 60.502225246029234 である。

2*82*(250√7 + 881)/4181

   60.502225246029234

   r2 = 30.2511;  x2 = -12.527;  y2 = 47.2674;  x3 = 31.7645;  y3 = 28.125
   乙円の直径 = 60.5022

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r3) = (100, 36) .// 2
   (r2, x2, y2, x3, y3) = res[1]
   #r2 = 4.30617;  x2 = 38.0693;  y2 = 6.72839;  x3 = 31.7645;  y3 = 28.125
   @printf("r2 = %g;  x2 = %g;  y2 = %g;  x3 = %g;  y3 = %g\n", r2, x2, y2, x3, y3)
   @printf("乙円の直径 = %g\n", 2r2)
   plot()
   circle(0, r1 - r3, r1)
   circle(0, r3 - r1, r1)
   circle(x2, y2, r2, :blue)
   circle(-x2, -y2, r2, :blue)
   circle(x3, y3, r3, :green)
   circle(-x3, -y3, r3, :green)
   circle(0, 0, r3, :green)
   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, r1 - r3, "甲円:r1 \n(0,r1-r3) ", :red, :right, :vcenter)
       point(0, r3 - r1, "(0,r3-r1)")
       point(x2, y2, "乙円:r2 \n(x2,y2) ", :blue, :right, :vcenter)
       point(x3, y3, " 丙円:r3(x3,y3)", :green, :center, :top, delta=-delta/2)
   end
end;

 

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

算額(その522)

2023年11月30日 | Julia

算額(その522)

岩手県山田町 武内大明神 文政3年(1820)
一関市博物館>>和算に挑戦>>平成28年度出題問題&解答例>>平成28年度出題問題(2)[中級問題] https://www.city.ichinoseki.iwate.jp/museum/wasan/h28/normal.html

外円内に長方形,甲円,乙円が入っている,乙円の直径が 1,長方形の短辺の長さが 6 のとき,外円の直径を求めよ。

長方形の長辺と短辺の長さをそれぞれ 2b, 2a とする。b = r0 - 2r1 である。
外円の半径と中心座標を r0, (0,  0)
甲円の半径と中心座標を r1, (b + r1, 0)
乙円の半径と中心座標を r2, (b + r2, y2)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r0::positive, r1::positive, r2::positive, y2::positive, b::positive;
a = r0 - 2r1
eq1 = (a + r2)^2 + y2^2 - (r0 - r2)^2
eq2 = (r1 - r2)^2 + y2^2 - (r1 + r2)^2
eq3 = a^2 - (r0^2 - b^2)
res = solve([eq1, eq2, eq3], (r0, r1, y2))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (b^2/(4*r2), b*(b - sqrt(b^2 - 16*r2^2))/(8*r2), sqrt(2)*sqrt(b)*sqrt(b - sqrt(b - 4*r2)*sqrt(b + 4*r2))/2)
    (b^2/(4*r2), b*(b + sqrt(b^2 - 16*r2^2))/(8*r2), sqrt(2)*sqrt(b)*sqrt(b + sqrt(b - 4*r2)*sqrt(b + 4*r2))/2)

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

   外円の直径 = 9;  r0 = 4.5;  r1 = 0.572949;  y2 = 1.07047

外円の直径は,長方形の短辺の長さの二乗を乙円の直径の 4 倍で割って得られる。

乙円の直径が 1,長方形の短辺の長さが 6 のとき,外円の直径は 36/4 = 9 である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, b) = (1//2, 6//2)
   (r0, r1, y2) = (b^2/(4*r2), b*(b - sqrt(b^2 - 16*r2^2))/(8*r2), sqrt(2)*sqrt(b)*sqrt(b - sqrt(b - 4*r2)*sqrt(b + 4*r2))/2)
   @printf("外円の直径 = %g;  r0 = %g;  r1 = %g;  y2 = %g\n", 2r0, r0, r1, y2)
   (r2, b) = (1//2, 6//2)
   a = r0 - 2r1
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:black, lw=0.5)
   circle(0, 0, r0)
   circle(a + r1, 0, r1, :blue)
   circle(-a - r1, 0, r1, :blue)
   circle4(a+r2, y2, r2, :magenta)
   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(a, 0, "a ", :black, :right, :bottom, delta=delta/2)
       point(0, b, " b", :black, :left, :bottom, delta=delta/2)
       point(a + r1, 0, "a+r1", :blue, :center, :bottom, delta=delta/2)
       point(a + r2, y2, "乙円:r2,(a+r2,y2) ", :magenta, :right, :vcenter)
   end
end;

 

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

算額(その521)

2023年11月29日 | Julia

算額(その521)

岩手県一関市 菅原神社 嘉永3年(1850)
一関市博物館>>和算に挑戦>>平成19年度出題問題&解答例>>平成19年度出題問題(2)[中級問題]
https://www.city.ichinoseki.iwate.jp/museum/wasan/h19/normal.html

正方形の中に四分円 3 個と等円 2 個が入っている。等円の直径が 1 のとき,正方形の一辺の長さを求めよ。また,斜線部分の面積を求めよ。

等円の半径と正方形の一辺の長さを r, a とおき,以下の方程式を解く。

include("julia-source.txt");

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

   16*r

r = 1/2 のとき,正方形の一辺の長さは 8 である。

求めるべき面積は A, B の和の 2 倍である。
A は半径 a の面積の 15/360 である。
B は半径 a の面積の 1/6 から,一辺の長さが a の正三角形の面積を差し引いたものである。

a = 8
A = PI*a^2*(15//360)
B = PI*a^2/6 - a*(a//2)√Sym(3)/2;

求める面積は (80π - 96√3)/3 = 28.3501782535237 である。

2(A + B) |> simplify |> println

   -32*sqrt(3) + 80*pi/3

2(A + B).evalf()

   28.3501782535237

function fill(a)
   x = vcat(a.*cosd.(270:330), [0])
   y = vcat(a.*sind.(270:330) .+ a, [0])
   plot!(x, y, color=:black, lw=0.5, seriestype=:shape, fillcolor=:blue, alpha=0.5)
   x = vcat(a.*cosd.(30:45), [0, cosd(30)])
   y = vcat(a.*sind.(30:45), [0, sind(30)])
   plot!(x, y, color=:black, lw=0.5, seriestype=:shape, fillcolor=:red, alpha=0.5)
end;

function draw(r, n, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r) = (8, 1//2)
   @printf("正方形の一辺の長さ = %g;  等円の半径 = %g\n", a, r)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
   circle(0, 0, a, beginangle=0, endangle=90)
   circle(a, 0, a, beginangle=90, endangle=180)
   circle(0, a, a, beginangle=270, endangle=360)
   circle(a/2, a-r, r, :brown)
   circle(a - r, a/2, r, :brown)
   fill(a)
   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.5a, 0.24a, "A", :black, mark=false)
       point(0.5a, 0.4a, "B", :black, mark=false)
       point(a/2, a - r, " (a/2,a-r)", :black, :center, :bottom, delta=delta/2)
       point(a - r, a/2, " (a-r,a/2)", :black, :center, :bottom, delta=delta/2)
       point(a, 0, "a ", :black, :right, :bottom, delta=delta/2)
   end
end;

 

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

算額(その520)

2023年11月29日 | Julia

算額(その520)

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

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

大きさの同じ球を下段に7個,中段に3個,上段に1個,互いに接するように積み上げる。球の直径が 1 のとき,盤上から一番上の球の天辺までの高さを求めよ。

球の半径を r とする。最上段の 4 個の球の中心を結んでできる正四面体を考える。



x^2 + y^2 + z^2 = BA^2, y = r, x = y/√3, BA = 2r である。正四面体の高さは sqrt(x^2 + y^2 + z^2) - 2r を z について解けば求まる。
z = r*2√6/3 である。
下段の3個の球の中心を通る水平面と上の球の中心を通る水平面の距離は z である。
設問ではこの3個の球のそれぞれが上段になりその下にある3個が構成する同様の構造体がある。最下段の 7 個の球の中心とその上の3個の球の中心を通る2つお水平面の距離も z である。
最下段の中心を通る水平面と原点の距離は r,最上段の球の中心と天辺との距離も r である。
したがって,原点と最上段の球の天辺との距離は r + z + z + r = 2r(1 + 2√6/3)
r = 1/2 のとき,1 + 2√6/3 = 2.6329931618554516 である。

include("julia-source.txt");

using SymPy
@syms r::positive, x, y, z
y = r
x = y/sqrt(Sym(3))
eq = sqrt(x^2 + y^2 + z^2) - 2r
z = solve(eq, z)[2] |> simplify
z |> println

   2*sqrt(6)*r/3

2r + 2z |> simplify |> println

   2*r*(3 + 2*sqrt(6))/3

(2z + 2r)(r => 1//2) |> println

   1 + 2*sqrt(6)/3

1 + 2*sqrt(6)/3

   2.6329931618554516

 

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

算額(その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でシェアする

算額(その516)

2023年11月27日 | Julia

算額(その516)

『算法新書』 文政13年(1830)
一関市博物館>>和算に挑戦>>平成26年度出題問題>>平成26年度出題問題(2)[中級問題](中学・高校生向き)

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

外円の中に弦 2 本と小円が 3 個入っている。左右の小円は弦にその中点で接している。外円の直径が 1 寸のとき,小円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
小円の半径と中心座標を r, (x, y), (0, r - R)
弦と円の右下の交点座標を (a, -sqrt(R^2 - a^2))
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R, r, a, x, y

eq1 = r/(2R - r) - (R - 2r)/R;

eq1 だけで,R が与えられたときに r を求めることができる。

solve(eq1, r) |> println

   Sym[R*(3 - sqrt(5))/2, R*(sqrt(5) + 3)/2]

r = R*(3 - sqrt(5))/2 が適解である。
すなわち,外円の直径が 1 のとき,小円の直径は (3 - sqrt(5))/2 ≒ 0.3819660112501051 である。

(3 - sqrt(5))/2

   0.3819660112501051

図を描くために必要な a, x, y を求めるには以下の eq2, eq3, eq4 を加えて,4元連立方程式を解けばよい。

eq2 = 2sqrt(R^2 - (R - 2r)^2) - sqrt((R + sqrt(R^2 - a^2))^2 + a^2)
eq3 = x^2 + y^2 - (R - r)^2
eq4 = (y/x)*((R + sqrt(R^2 - a^2))/a) + 1
res = solve([eq1, eq2, eq3, eq4], (r, a, x, y));

16 組の解が得られるが,6 番目のものが適解である。

res[6][1] |> println  # r

   R*(3 - sqrt(5))/2

res[6][2] |> println  # a

   4*R*sqrt(-38 + 17*sqrt(5))

res[6][3] |> println  # x

   sqrt(2)*sqrt(R^2*(-8*R*(-199 + 89*sqrt(5))/(R + sqrt(609 - 272*sqrt(5))*sqrt(R^2)) - sqrt(5) + 3))/2

res[6][4] |> println  # y

   2*sqrt(R^3/(R + sqrt(609 - 272*sqrt(5))*sqrt(R^2)))*sqrt(-199 + 89*sqrt(5))

res[6]

   (R*(3 - sqrt(5))/2, 4*R*sqrt(-38 + 17*sqrt(5)), sqrt(2)*sqrt(R^2*(-8*R*(-199 + 89*sqrt(5))/(R + sqrt(609 - 272*sqrt(5))*sqrt(R^2)) - sqrt(5) + 3))/2, 2*sqrt(R^3/(R + sqrt(609 - 272*sqrt(5))*sqrt(R^2)))*sqrt(-199 + 89*sqrt(5)))

   R = 1;  r = 0.381966;  a = 0.458792;  x = 0.600566;  y = 0.145898

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 1
   (r, a, x, y) = (
       R*(3 - sqrt(5))/2,
       4*R*sqrt(-38 + 17*sqrt(5)),
       sqrt(2)*sqrt(R^2*(-8*R*(-199 + 89*sqrt(5))/(R + sqrt(609 - 272*sqrt(5))*sqrt(R^2)) - sqrt(5) + 3))/2,
       2*sqrt(R^3/(R + sqrt(609 - 272*sqrt(5))*sqrt(R^2)))*sqrt(-199 + 89*sqrt(5)))
   @printf("R = %g;  r = %g;  a = %g;  x = %g;  y = %g\n", R, r, a, x, y)
   plot()
   circle(0, 0, R)
   circle(x, y, r, :blue)
   circle(-x, y, r, :blue)
   circle(0, r - R, r, :blue)
   segment(0, R, a, -sqrt(R^2 - a^2), :green)
   segment(0, R, -a, -sqrt(R^2 - a^2), :green)
   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(x, y, " (x,y)", :blue, :left, :vcenter)
       point(a, -√(R^2 - a^2), "(a, -√(R^2-a^2))", :blue)
       point(0, r - R, " r-R", :blue, :left, :vcenter)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta)
   end
end;

 

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

算額(その515)

2023年11月26日 | Julia

算額(その515)

横川良助直胤: 神壁算法追加,文化4年(1807)
一関市博物館>>和算に挑戦>>平成27年度出題問題>>平成27年度出題問題(2)[中級問題]&解答例
https://www.city.ichinoseki.iwate.jp/museum/wasan/h27/normal.html

上底と下底が 4寸5分,8寸の等脚台形の中に円が内接しており,その円に正方形が内接しており,さらにその正方形に円が内接しており...と繰り返されている。
。4 番目の円に内接する正方形の一辺の長さを求めよ。なお,最後の正方形の辺は座標軸に水平・垂直である。

上底,仮定を 2a, 2b とする。
台形の高さ h は sqrt((a + b)^2 - (b - a)^2) である。
よって,第1円の半径は h/2 である。
第1正方形の対角線の長さが第1円の半径に等しい。
第n円の直径は第n-1円の直径の 1/√2 である。
以下,この繰り返し。
最後の正方形は45度回転して描く。

   a = 2.25;  b = 4;  h = 6;  r1 = 3
   1 番目の内接円の半径 = 3
   2 番目の内接円の半径 = 2.12132
   3 番目の内接円の半径 = 1.5
   4 番目の内接円の半径 = 1.06066
   最後の正方形の一辺の長さは 1.5(1 寸 5 分)

include("julia-source.txt");

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (45, 80) .// 20
   h = sqrt((a + b)^2 - (b - a)^2)
   r1 = h/2
   r = zeros(10)
   @printf("a = %g;  b = %g;  h = %g;  r1 = %g\n", a, b, h, r1)
   plot([b, a, -a, -b, b], [0, h, h, 0, 0], color=:black, lw=0.5)
   r[1] = h/2
   n = 4
   for i = 1:n - 1
       @printf("%d 番目の内接円の半径 = %g\n", i, r[i])
       circle(0, h/2, r[i])
       plot!([r[i], 0, -r[i], 0, r[i]], [h/2, h/2 + r[i], h/2, h/2 - r[i], h/2], color=:blue, lw=0.5)
       r[i + 1] = r[i]/√2
   end
   @printf("%d 番目の内接円の半径 = %g\n", n, r[n])
   circle(0, h/2, r[n])
   t = r[n]*√2/2
   plot!([t, -t, -t, t], [h/2 + t, h/2 + t, h/2 - t, h/2 - t, h/2 + t], color=:green, lw=0.5)
   @printf("最後の正方形の一辺の長さは %g\n", 2t)
   if more == true
       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(a, h, " (a,h)", :black, :left, :vcenter)
       point(b, 0, "b ", :black, :right, :bottom, delta=delta/2)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その514)

2023年11月26日 | Julia

算額(その514)

横川良助直胤: 神壁算法追加,文化4年(1807)
一関市博物館>>和算に挑戦>>平成27年度出題問題>>平成27年度出題問題(2)[中級問題]&解答例
https://www.city.ichinoseki.iwate.jp/museum/wasan/h27/normal.html

大円と小円が交わっており,上部の小円部分に甲円 1 個,乙円 2 個が入っている。

下図は説明のためのものなので,問の条件のときのものとは見た目が違う。


小円,甲円,乙円の直径がそれぞれ,8 寸,5 寸,3 寸のとき,大円の直径はいかほどか。

大円の直径と中心座標を r0, (0, 0)
小円の直径と中心座標を r1, (0, r0 + 2r2 - r1)
甲円の直径と中心座標を r2, (0, r0 + r2)
乙円の直径と中心座標を r3, (x3, y3)
として,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r0::positive, r1::positive, r2::positive, r3::positive, x3::positive, y3::positive

(r1, r2, r3) = (8, 3, 2) .// 2
eq1 = x3^2 + y3^2 - (r0 + r3)^2
eq2 = x3^2 + (y3 - (r0+ 2r2 - r1))^2 - (r1 - r3)^2
eq3 = x3^2 + (r0 + r2 - y3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (r0, x3, y3))

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

2 組の解が得られるが,2 番目のものが適解である。
小円,甲円,乙円の直径がそれぞれ,8 寸,5 寸,3 寸のとき,大円の直径は 27 寸である。

(r1, r2, r3) = (8, 3, 2) .// 2
(
   -r2^2*(r1 - r2 - r3)/(r1*r2 - r1*r3 - r2^2),
   2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 - r2 - r3)/(r1 - r2),
   (r1^2*r2^2 - 3*r1^2*r2*r3 + r1^2*r3^2 - 2*r1*r2^3 + 3*r1*r2^2*r3 + r1*r2*r3^2 + r2^4)/((r1 - r2)*(r1*r2 - r1*r3 - r2^2))
)

   (27//2, 2.3999999999999995, 143//10)

using Plots
function draw(r1, r2, r3, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, x3, y3) = (-r2^2*(r1 - r2 - r3)/(r1*r2 - r1*r3 - r2^2), 2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 - r2 - r3)/(r1 - r2), (r1^2*r2^2 - 3*r1^2*r2*r3 + r1^2*r3^2 - 2*r1*r2^3 + 3*r1*r2^2*r3 + r1*r2*r3^2 + r2^4)/((r1 - r2)*(r1*r2 - r1*r3 - r2^2)))
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  r0 = %g;  x3 = %g;  y3 = %g\n", r1, r2, r3, r0, x3, y3)
   @printf("大円の直径 = %g\n", 2r0)
   plot()
   circle(0, 0, r0)
   circle(0, r0 + 2r2 - r1, r1, :blue)
   circle(0, r0 + r2, r2, :magenta)
   circle(x3, y3, r3, :green)
   circle(-x3, y3, r3, :green)
   if more == true
       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+r2, " r0+r2", :magenta, :left, :vcenter)
       point(0, r0 + 2r2 - r1, " r0+2r2-r1")
       point(x3, y3, "(x3,y3)")
       point(0, r0, "r0 ", :red, :right, :top, delta=-delta)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その513)

2023年11月26日 | Julia

算額(その513)

岩手県盛岡市 八幡宮 文化6年(1869)
一関市博物館>>和算に挑戦>>平成27年度出題問題>>平成27年度出題問題(2)[中級問題]&解答例
https://www.city.ichinoseki.iwate.jp/museum/wasan/h27/normal.html

大円内に 8 個の小円が入っている。大円の直径が 25 寸のとき,小円の直径を求めよ。

大円の直径と中心座標を R, (0, 0)
小円の直径と中心座標を r, (r, 2√3r), (0, √3r), (r, 0)
として,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive

eq = r^2 + (2sqrt(Sym(3))r)^2 - (R - r)^2
solve(eq, r)[1] |> println

   R*(-1 + sqrt(13))/12

小円の半径は大円の半径の (√13 - 1)/12 倍である。

   R = 12.5;  r = 2.71412;  小円の直径 = 5.42823

大円の直径が 25 寸のとき小円の直径は約 5.42823 寸である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 25/2
   r = R*(√13 - 1)/12
   @printf("R = %g;  r = %g;  小円の直径 = %g\n", R, r, 2r)
   plot()
   circle(0, 0, R)
   circle4(r, 2√3r, r, :blue)
   circle(0, √3r, r, :blue)
   circle(0, -√3r, r, :blue)
   circle(r, 0, r, :blue)
   circle(-r, 0, r, :blue)
   if more == true
       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, 0, "r", :blue, :center, :bottom, delta=delta/2)
       point(0, √3r, " √3r", :blue, :left, :vcenter)
       point(r, 2√3r, " 小円:r\n (r,2√3r)", :blue, :center, :top, delta=-delta)
       point(R, 0, "R ", :red, :right, :bottom, delta=delta/2)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その512)

2023年11月26日 | Julia

算額(その512)

兵庫県上郡町 上郡天満神社 文化3年(1806)
https://kamigori-town-museum.jp/cultural-asetts/559/
藤井康生: 上郡町天満神社算額見学記,和算,第93号,2002年1月31日
http://www.wasan.jp/kinkikaisi/wasan093.pdf

外円内に弦と矢および甲円,乙円が入っている。甲円と乙円の直径の和が 1尺4分,外円の直径と矢の差が 1尺4寸4分のとき矢を求めよ。

弦と y 軸の交点座標を (0, a) とおく。
外円の直径と中心座標を r0, (0, 0)
甲円の直径と中心座標を r1, (0, a + r1)
乙円の直径と中心座標を r2, (0, a + r2)
として,以下の連立方程式を解く。
矢(の長さ)は r0 - a である。

include("julia-source.txt");

using SymPy
@syms a::positive, r0::positive, r1::positive, r2::positive, x2::positive

eq1 = 2r1 + 2r2 - 104//10
eq2 = 2r0 - (r0 - a) - 144//10
eq3 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq4 = r1^2 + (a + r1)^2 - (r0 - r1)^2
eq5 = x2^2 + (a + r2)^2 - (r0 - r2)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (a, r0, r1, r2, x2))

   1-element Vector{NTuple{5, Sym}}:
    (63/20, 45/4, 18/5, 8/5, 42/5)

   a = 3.15;  r0 = 11.25;  r1 = 3.6; r2 = 1.6;  x2 = 8.4

矢 = res[1][2] - res[1][1]
矢 |> println

   81/10

矢は 8寸1分 である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r0, r1, r2, x2) = res[1]
   @printf("a = %g;  r0 = %g;  r1 = %g; r2 = %g;  x2 = %g\n",
       a, r0, r1, r2, x2)
   @printf("矢 = %g\n", r0 - a)
   plot()
   circle(0, 0, r0)
   circle(r1, a + r1, r1, :blue)
   circle(x2, a + r2, r2, :green)
   segment(-sqrt(r0^2 - a^2), a, sqrt(r0^2 - a^2), a, :brown)
   if more == true
       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(r1, a + r1, " 甲円:r1\n (r1,a+r1)", :blue, :left, :vcenter)
       point(x2, a + r2, "乙円:r2,(x2,a+r2) ", :black, :right, :vcenter)
       point(0, a, "a ", :brown, :right, :bottom, delta=delta/2)
       point(0, r0, "r0 ", :red, :right, :bottom, delta=delta/2)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その511)

2023年11月26日 | Julia

算額(その511)

長野県長野市安茂里正覚院 久保寺観音堂 文政13年(1830)
中村信弥「改訂増補 長野県の算額」

http://www.wasan.jp/zoho/zoho.html

半円内に2本の斜線を引き,甲円 2 個と乙円を入れる。乙円の径が 1 寸のとき,甲円の径はいかほどか。

各円の半径,中心座標を以下のようにする。
外円: R, (0, 0)
甲円: r1, (x1, y1) および (x2, r1)
乙円: r3, (x3, y3)
外円の円周上にある2本の斜線の交点座標を (x, y)
とおき,以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms x1, y1, r1, x2, r3, x3, y3, x, R

r3 = 1//2
y = sqrt(R^2 - x^2)  # (x, y) は外円の円周上にある
eq1 = x1^2 + y1^2 - (R - r1)^2  # 右の甲円が外円に内接する
eq2 = x3^2 + y3^2 - (R - r3)^2  # 乙円が外円に内接する
eq3 = sqrt((R + x)^2 + y^2) + sqrt((R - x)^2 + y^2) - 2R - 2r1  # 鉤股弦に内接する円の直径
eq4 = (y3/x3)*(y/(x + R)) + 1
eq5 = numerator(distance(-R, 0, x, y, x2, r1) - r1^2) |> simplify
eq6 = numerator(distance(R, 0, x, y, x1, y1) - r1^2) |> simplify
eq7 = ((x - R)/2 - x3)^2 + (y/2 - y3)^2 - r3^2
eq8 = (y1/x1) * (y/(x - R)) + 1;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (x1, y1, r1, x2, x3, y3, x, R))

using Printf
using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (x1, y1, r1, x2, x3, y3, x, R) = u
   t = sqrt(R^2 - x^2)
   return [
       x1^2 + y1^2 - (R - r1)^2,  # eq1
       x3^2 + y3^2 - (R - 1/2)^2,  # eq2
       -2*R - 2*r1 + sqrt(R^2 - x^2 + (R - x)^2) + sqrt(R^2 - x^2 + (R + x)^2),  # eq3
       1 + y3*sqrt(R^2 - x^2)/(x3*(R + x)),  # eq4
       (-4*R^2*r1^2 + (R*r1 - R*t + r1*x - x2*t)^2 + (R^2 - R*x + R*x2 - r1*t - x*x2)^2)/(4*R^2),  # eq5
       (-4*R^2*r1^2 + (R*y1 - R*t - x*y1 + x1*t)^2 + (-R^2 - R*x + R*x1 + x*x1 + y1*t)^2)/(4*R^2),  # eq6
       (-y3 + t/2)^2 + (-R/2 + x/2 - x3)^2 - 1/4,  # eq7
       1 + y1*t/(x1*(-R + x)),  # eq8
   ]
end;

iniv = BigFloat[1.7, 4.2, 2.1, -3.5, -5.8, 2.4, -4.6, 6.4]
res = nls(H, ini=iniv)
println(res)

   (BigFloat[3.461538461538461538461538509169459561367959236632964526959234465030510845593508, 8.307692307692307692307692429999766197094115669634049393814154165178838427207838, 4.000000000000000000000000058695639811160390953981958412762403947086964336791015, -7.000000000000000000000000103675105612121327124114588925459978536567327493943071, -11.53846153846153846153846171778119379511106203957672690108302039508364569817559, 4.807692307692307692307692375277671136103845336879392455246078657573883943485054, -9.153846153846153846153846290689476553117442078831806765600936112963486841120747, 13.00000000000000000000000019241303418103091110455672021661459202029874280726291], true)

   r3 = 0.5;  x1 = 3.46154;  y1 = 8.30769; r1 = 4;  x2 = -7;  x3 = -11.5385;  y3 = 4.80769;  x = -9.15385;  y = 9.23077;  R = 13
   甲円の径 = 2r1 = 8

甲円の半径は 4 である。元の単位では甲円の直径は 8 寸である。

using Plots
function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 1//2
   (x1, y1, r1, x2, x3, y3, x, R) = res[1]
   y = sqrt(R^2 - x^2)
   @printf("r3 = %g;  x1 = %g;  y1 = %g; r1 = %g;  x2 = %g;  x3 = %g;  y3 = %g;  x = %g;  y = %g;  R = %g\n",
       r3, x1, y1, r1, x2, x3, y3, x, y, R)
   @printf("甲円の径 = 2r1 = %g\n", 2r1)
   plot([R, x, -R, R], [0, y, 0, 0], color=:brown, lw=0.5)
   circle(0, 0, R, beginangle=0, endangle=180)
   circle(x1, y1, r1, :blue)
   circle(x2, r1, r1, :green)
   circle(x3, y3, r3, :magenta)
   if more == true
       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(x1, y1, " 甲円:r1\n (x1,y1)", :blue, :left, :vcenter)
       point(x2, r1, " 甲円:r1\n (x2,r1)", :blue, :left, :top, delta=-delta)
       point(x3, y3, "  乙円:r3,(x3,y3)", :magenta, :left, :vcenter)
       point(x, y, "(x,y) ", :brown, :right, :vcenter)
   else
       plot!(showaxis=false)
   end
end;

 

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

算額(その510)

2023年11月25日 | Julia

算額(その510)

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション

https://dl.ndl.go.jp/pid/3508431/1/11

(16) 京都府京都市右京区山ノ内宮脇町 山王神社 明治23年(1890)
近畿数学史学会:近畿の算額「数学の絵馬を訪ねて」,平成4年5月16日 初版第一刷,大阪教育図書株式会社,大阪市.

正方形内に甲円 1 個,乙円 2 個が入っている。甲円の直径が 9 寸のとき,正方形の一辺の長さはいかほどか。

正方形の一辺の長さを a = 4r2
甲円の半径と中心座標を r1, (0, 4r2 - r1)
乙円の半径と中心座標を r2, (r2, r2)
として方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive
eq = r2^2 + (3r2 - r1)^2 - (r1 + r2)^2
res = solve(eq, r2)[1]
res |> println

   8*r1/9

乙円の半径は甲円の半径の 8/9 である。
甲円の直径が 9 寸ならば乙円の直径は 8 寸,正方形の一辺の長さはその 2 倍,16 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 9/2
   r2 = 8r1/9
   a = 4r2
   @printf("r1 = %g;  r2 = %g;  正方形の一辺の長さ = %g\n", r1, r2, a)
   plot([a/2, a/2, -a/2, -a/2, a/2], [0, a, a, 0, 0], color=:blue, lw=0.5)
   circle(r2, r2, r2, :green)
   circle(-r2, r2, r2, :green)
   circle(0, a - r1, r1)
   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, a - r1, " 甲円:r1,(0,a-r1)", :red, :left, :vcenter)
       point(r2, r2, " 乙円:r2(r2,r2)", :green, :left, :vcenter)
       point(a/2, 0, "a/2=2r2 ", :blue, :right, :bottom, delta=delta)
   else
       plot!(showaxis=false)
   end
end;

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

算額(その509)

2023年11月25日 | Julia

算額(その509)

岩手県一関市上大槻街路 一関神明社 天保5年(1834)
http://www.wasan.jp/iwate/sinmeisha.html

正方形内に,四分円 2 個,半円 1 個,天円,地円が入っている。
地円の直径が 1 寸のとき,天円の直径はいかほどか。

正方形の一辺の長さを a
天円の半径と中心座標を r1, (r1, a - r1)
地円の半径と中心座標を r2, (x2, r2)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive, x2::positive

eq1 = x2^2 + (a - r2)^2 - (a + r2)^2
eq2 = 2(a - r1)^2 - (a + r1)^2
eq3 = (a - x2)^2 + (a//2 -r2)^2 - (a//2 + r2)^2
res = solve([eq1, eq2, eq3], (r1, r2, x2))

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

4 組の解が得られるが,1 番目のものが適解である。
r1 は r2 の 2 倍である。
地円の直径が 1 寸のとき,天円の直径は 2 寸である。

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

   2

ちなみに,地円の直径が 1 寸となるのは,正方形の一辺の長さが 4*sqrt(2) + 6 のときである。

   a = 11.6569;  r1 = 2;  r2 = 1;  x2 = 6.82843

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 4*sqrt(2) + 6
   (r1, r2, x2) = (a*(3 - 2*sqrt(2)), -a/2 + a*(2 - sqrt(2)), a*(2 - sqrt(2)))
   @printf("a = %g;  r1 = %g;  r2 = %g;  x2 = %g\n", a, r1, r2, x2)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:blue, lw=0.5)
   circle(0, a, a, :green, beginangle=270, endangle=360)
   circle(a, 0, a, :green, beginangle=90, endangle=180)
   circle(a, a/2, a/2, :magenta, beginangle=90, endangle=270)
   circle(r1, a - r1, r1)
   circle(x2, r2, r2, :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(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
       point(0, a, " a", :blue, :left, :bottom, delta=delta/2)
       point(a, a/2, "(a,a/2) ", :magenta, :right, :vcenter)
       point(r1, a - r1, " 天円:r1\n (r1,a-r1)", :red, :left, :vcenter)
       point(x2, r2, "地円:r2,(x2,r2)", :black, :center, :top, delta=-delta/2)
   else
       plot!(showaxis=false)
   end
end;

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

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

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