裏 RjpWiki

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

算額(その1071)

2024年06月16日 | Julia

算額(その1071)

九十六 大船渡市立根 気仙安養寺稲荷堂 文政5年(1822)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

全円の中に甲円と心円,および乙円から始まる累円(丙円,丁円,戊円,己円...)が入っている。甲円の直径は全円の半径に等しい。心円の直径が 8 寸 9 分のとき,己円の直径を求めよ。

全円の半径と中心座標を R, (0, 0); R = 2r1
甲円の半径と中心座標を r1, (0, r1)
心円の半径と中心座標を r2, (0, -r2)
乙円の半径と中心座標を r3, (x3, y3); x3 = r3,y3 < 0
丙円の半径と中心座標を r4, (x4, y4)
とおく。

まず,甲円の半径 r1 と乙円のパラメータ r3, x3, y3 を求める。

include("julia-source.txt")

using SymPy

@syms R::positive, r1::positive,
     r3::positive, y3::negative,
     r2::positive
R = 2r1
eq1 = r3^2 + (r1 - y3)^2 - (r1 + r3)^2
eq2 = r3^2 + (r2 + y3)^2 - (r3 + r2)^2
eq3 = r3^2 + y3^2 - (R - r3)^2
res = solve([eq1, eq2, eq3], (r1, r3, y3))[1]

   (7*r2, 56*r2/9, -14*r2/3)

累円は,「ある円に外接する,半径が小さい次の円」の連続である。
乙円と丙円の関係は全円と甲円 の半径により,以下の連立方程式で得られる漸化式で計算される。

@syms R::positive, r1::positive,
     r3::positive, y3::negative,
     r2::positive
R = 2r1
@syms r4, x4, y4, x3
#x3 = r3
#(r1, r2, y3) =  (7*r2, 56*r2/9, -14*r2/3)
eq4 = x4^2 + (r1 - y4)^2 - (r1 + r4)^2
eq5 = x4^2 + y4^2 - (R - r4)^2
eq6 = (x4 - x3)^2 + (y4 - y3)^2 - (r3 + r4)^2
res = solve([eq4, eq5, eq6], (r4, x4, y4))[1]

   ((8*r1^3 + 4*r1^2*r3 - 20*r1^2*y3 - 2*r1*r3^2 - 4*r1*r3*y3 + 10*r1*x3^2 + 14*r1*y3^2 - r3^3 + 3*r3^2*y3 + r3*x3^2 + r3*y3^2 - 3*x3^2*y3 - 2*sqrt(2)*x3*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4) - 3*y3^3)/(2*(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2)), (8*r1^2*x3 - 4*r1*r3*x3 - 4*r1*x3*y3 + 2*sqrt(2)*r1*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4) - 4*r3^2*x3 + sqrt(2)*r3*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4) + 4*x3^3 + 4*x3*y3^2 - 3*sqrt(2)*y3*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4))/(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2), 3*sqrt(2)*x3*sqrt(-(-4*r1^2 - 4*r1*r3 - r3^2 + x3^2 + y3^2)*(2*r1*r3 - 2*r1*y3 - r3^2 + x3^2 + y3^2))/(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2) + (-8*r1^3 + 4*r1^2*r3 + 12*r1^2*y3 + 10*r1*r3^2 - 12*r1*r3*y3 + 2*r1*x3^2 - 6*r1*y3^2 + 3*r3^3 - 9*r3^2*y3 - 3*r3*x3^2 - 3*r3*y3^2 + 9*x3^2*y3 + 9*y3^3)/(2*(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2)))

これにより得られる式は,丙円と丁円,丁円と己円...にも同じように適用される。以下の関数は,ある円 r3, x3, y3 の次の円の 半径,中心の x, y 座標を返す。
これを次々に適用することにより,累円のパラメータを求めることができる。

nextcircle(r1, r3, x3, y3) = ((8*r1^3 + 4*r1^2*r3 - 20*r1^2*y3 - 2*r1*r3^2 - 4*r1*r3*y3 + 10*r1*x3^2 + 14*r1*y3^2 - r3^3 + 3*r3^2*y3 + r3*x3^2 + r3*y3^2 - 3*x3^2*y3 - 2*sqrt(2)*x3*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4) - 3*y3^3)/(2*(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2)), (8*r1^2*x3 - 4*r1*r3*x3 - 4*r1*x3*y3 + 2*sqrt(2)*r1*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4) - 4*r3^2*x3 + sqrt(2)*r3*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4) + 4*x3^3 + 4*x3*y3^2 - 3*sqrt(2)*y3*sqrt(8*r1^3*r3 - 8*r1^3*y3 + 4*r1^2*r3^2 - 8*r1^2*r3*y3 + 4*r1^2*x3^2 + 4*r1^2*y3^2 - 2*r1*r3^3 - 2*r1*r3^2*y3 + 2*r1*r3*x3^2 + 2*r1*r3*y3^2 + 2*r1*x3^2*y3 + 2*r1*y3^3 - r3^4 + 2*r3^2*x3^2 + 2*r3^2*y3^2 - x3^4 - 2*x3^2*y3^2 - y3^4))/(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2), 3*sqrt(2)*x3*sqrt(-(-4*r1^2 - 4*r1*r3 - r3^2 + x3^2 + y3^2)*(2*r1*r3 - 2*r1*y3 - r3^2 + x3^2 + y3^2))/(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2) + (-8*r1^3 + 4*r1^2*r3 + 12*r1^2*y3 + 10*r1*r3^2 - 12*r1*r3*y3 + 2*r1*x3^2 - 6*r1*y3^2 + 3*r3^3 - 9*r3^2*y3 - 3*r3*x3^2 - 3*r3*y3^2 + 9*x3^2*y3 + 9*y3^3)/(2*(4*r1^2 + 4*r1*r3 - 12*r1*y3 + r3^2 - 6*r3*y3 + 8*x3^2 + 9*y3^2)));

甲円の初期パラメータ (r3, x3, y3) = (27.68888888888889, 27.68888888888889, -20.76666666666667) をはじめとして,累円のパラメータを実際に求めてみる。

nextcircle(31.15, 14.658823529411768, 43.976470588235316, 18.3235294117647)

   (7.55151515151515, 37.75757575757576, 39.645454545454534)

nextcircle(31.15, 7.55151515151515, 37.75757575757576, 39.645454545454534)

   (4.371929824561412, 30.603508771929828, 49.18421052631578)

プログラムで書いてみると以下のようになる。

r2 = 8.9/2
(r1, r3, y3) = (7*r2, 56*r2/9, -14*r2/3)
R = 2r1
(r, x, y) = (r3, r3, y3)
for i = 1:10
   println((r, x, y))
   (r, x, y) = nextcircle(31.15, r, x, y)
end

   (27.68888888888889, 27.68888888888889, -20.76666666666667)
   (14.658823529411762, 43.976470588235294, 18.32352941176471)
   (7.551515151515154, 37.75757575757576, 39.64545454545455)
   (4.371929824561407, 30.60350877192981, 49.18421052631579)
   (2.7999999999999976, 25.199999999999978, 53.899999999999984)
   (1.931782945736425, 21.24961240310076, 56.5046511627907)
   (1.4079096045197808, 18.302824858757038, 58.07627118644067)
   (1.0695278969957025, 16.042918454935606, 59.09141630901287)
   (0.8390572390572382, 14.263973063973097, 59.782828282828284)
   (0.6753387533875355, 12.83143631436317, 60.2739837398374)

己円の半径は 2.7999999999999976 である(直径は 5.6 寸)。

なお,累円の半径は分数で表すことができ,分子は 1246 で,分母は g(n) = 20n^2 - 20n + 45 である。

using Printf
g(n) = 20n^2 - 20n + 45;
for n = 1:10
   @printf("n = %2d,  分母 = %4d,  半径 = %.15g\n", n, g(n), 1246/g(n))
end

   n =  1,  分母 =   45,  半径 = 27.6888888888889
   n =  2,  分母 =   85,  半径 = 14.6588235294118
   n =  3,  分母 =  165,  半径 = 7.55151515151515
   n =  4,  分母 =  285,  半径 = 4.3719298245614
   n =  5,  分母 =  445,  半径 = 2.8
   n =  6,  分母 =  645,  半径 = 1.93178294573643
   n =  7,  分母 =  885,  半径 = 1.40790960451977
   n =  8,  分母 = 1165,  半径 = 1.06952789699571
   n =  9,  分母 = 1485,  半径 = 0.839057239057239
   n = 10,  分母 = 1845,  半径 = 0.675338753387534

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   names = Char["甲乙丙丁戊己庚辛壬癸"...]
   r2 = 8.9/2
   (r1, r3, y3) = (7*r2, 56*r2/9, -14*r2/3)
   R = 2r1
   @printf("心円の直径が %g のとき,己円の直径は %g である。\n", 2r2, 2r3)
   plot(showaxis=false)
   delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
   circlef(0, 0, R, :orange)
   circlef(0, r1, r1, :green)
   circlef(0, -r2, r2, :purple)
   (rnow, xnow, ynow) = (r3, r3, y3)
   for i = 1:9
       @printf("%s: R = %.15g;  r1 = %.15g;  r = %.15g;  x = %.15g;  y = %.15g\n", names[i+1], R, r1, rnow, xnow, ynow)
       circlef(xnow, ynow, rnow, i)
       circlef(-xnow, ynow, rnow, i)
       point(xnow, ynow, names[i+1], :white, :center, :vcenter, mark=false)
       (rnow, xnow, ynow) = nextcircle(31.15, rnow, xnow, ynow)
   end
   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(0, r1, "甲", :white, :center, :vcenter, mark=false)
       point(0, -r2, "心", :white, :center, :vcenter, mark=false)
   end
end;

   心円の直径が 8.9 のとき,己円の直径は 55.3778 である。
   乙: R = 62.3;  r1 = 31.15;  r = 27.6888888888889;  x = 27.6888888888889;  y = -20.7666666666667
   丙: R = 62.3;  r1 = 31.15;  r = 14.6588235294118;  x = 43.9764705882353;  y = 18.3235294117647
   丁: R = 62.3;  r1 = 31.15;  r = 7.55151515151515;  x = 37.7575757575758;  y = 39.6454545454545
   戊: R = 62.3;  r1 = 31.15;  r = 4.37192982456141;  x = 30.6035087719298;  y = 49.1842105263158
   己: R = 62.3;  r1 = 31.15;  r = 2.8;  x = 25.2;  y = 53.9
   庚: R = 62.3;  r1 = 31.15;  r = 1.93178294573643;  x = 21.2496124031008;  y = 56.5046511627907
   辛: R = 62.3;  r1 = 31.15;  r = 1.40790960451978;  x = 18.302824858757;  y = 58.0762711864407
   壬: R = 62.3;  r1 = 31.15;  r = 1.0695278969957;  x = 16.0429184549356;  y = 59.0914163090129
   癸: R = 62.3;  r1 = 31.15;  r = 0.839057239057238;  x = 14.2639730639731;  y = 59.7828282828283

デカルトの円定理を用いる方法

円の直径(半径)だけを求めるならば,デカルトの円定理を使って,累円の半径を求めてゆくことができる。
半径 r1, prev, nxt の 3 円が内接する外円の半径を R とすると,以下の式が成り立つ。

using SymPy
@syms R, r1, prev, nxt
nxtr = -1/(1/r1 + 1/prev + 1/nxt - 2sqrt(1/(r1*prev) + 1/(prev*nxt) + 1/(nxt*r1))) - R

これを解いて nxt を求める式を得る。

result = solve(nxtr, nxt)[1]
result |> println

   (R*prev*r1*(R*prev + R*r1 - prev*r1) - 2*sqrt(-R^3*prev^3*r1^3*(-R + prev + r1)))/(R^2*prev^2 - 2*R^2*prev*r1 + R^2*r1^2 + 2*R*prev^2*r1 + 2*R*prev*r1^2 + prev^2*r1^2)

これを関数にする。この関数は,円の半径を入力すれば,その円に外接する次の円の半径を返す。

res2(prev, R, r1) = (R*prev*r1*(R*prev + R*r1 - prev*r1) - 2*sqrt(-R^3*prev^3*r1^3*(-R + prev + r1)))/(R^2*prev^2 - 2*R^2*prev*r1 + R^2*r1^2 + 2*R*prev^2*r1 + 2*R*prev*r1^2 + prev^2*r1^2)

   res2 (generic function with 1 method)

R, r1, r3 が求まった段階で,次々に関数を適用してゆく。

r2 = 8.9/2
r1 = 7*r2
R = 2r1
r3 = r2 * 56/9  # 乙円
(r2, R, r1, r3)

   (4.45, 62.300000000000004, 31.150000000000002, 27.68888888888889)

nxt = res2(r3, R, r1)  # 丙円

   14.65882352941177

nxt = res2(nxt, R, r1)  # 丁円

   7.551515151515152

nxt = res2(nxt, R, r1)  # 戊円

   4.371929824561404

nxt = res2(nxt, R, r1)  # 己円

   2.8000000000000003

nxt = res2(nxt, R, r1)    # 庚円

   1.9317829457364344

 

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

算額(その1070)

2024年06月16日 | Julia

算額(その1070)

九十六 大船渡市立根 気仙安養寺稲荷堂 文政5年(1822)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

正方形の中に直角三角形と,甲円,乙円,丙円,丁円を容れる。正方形の一辺の長さが 2 寸のとき,4 個の円の和はいかほどか。正方形の一辺の長さは弦(直角三角形の斜辺)の 4/5 である。

最後の一文は,条件なのかヒントなのか?

正方形の一辺の和を a
正方形の辺上にある正三角形の頂点座標を (b, a), (0, c)
直角三角形の 3 辺の長さを ab, bc, ca
甲円の半径と中心座標を r1, (a - r1, a - r1)
乙円の半径と中心座標を r2, (r2, r2)
丙円の半径と中心座標を r3, (x3, y3)
丁円の半径と中心座標を r4, (r4, a - r4)
とおき,以下の連立方程式を解く。

11 元連立方程式であるが,eq1〜eq9 までの条件式は簡単なものなので,SymPy でも簡単に 9 元連立方程式の解(b, c, ab, bc, ca, r1, r2, r3, r4)を求めることができる。
残るパラメータ (x3, y3) は eq10, eq11 の連立方程式を解けばよい。

include("julia-source.txt")

using SymPy

@syms a::positive, b::positive, c::positive,
     ab::positive, bc::positive, ca::positive,
     r1::positive, r2::positive, r3::positive, r4::positive,
     x3::positive, y3::positive;
eq1 = bc^2 + ca^2 - ab^2
eq2 = (a - b) + a - ab - 2r1
eq3 = c + a - ca - 2r2
eq4 = bc + ca - ab - 2r3
eq5 = b + (a - c) - bc - 2r4
eq6 = ((a - b) + a + ab)*r1 + (c + a + ca)*r2 + ( bc + ca + ab)*r3 + (b + (a - c) + bc)*r4 - 2a^2
eq6 = 4ab/5 - a
eq7 = (a - b)^2 + a^2 - ab^2
eq8 = c^2 + a^2 - ca^2
eq9 = b^2 + (a - c)^2 - bc^2
eq10 = dist2(a, 0, 0, c, x3, y3, r3)
eq11 = dist2(a, 0, b, a, x3, y3, r3)
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (b, c, ab, bc, ca, r1, r2, r3, r4))[1]

   (a/4, a/2, 5*a/4, sqrt(5)*a/4, sqrt(5)*a/2, a/4, a*(3 - sqrt(5))/4, a*(-5 + sqrt(5))/8 + sqrt(5)*a/4, a*(3/8 - sqrt(5)/8))

4 円の直径の和は 2(r1 + r2 + r3 + r4) = 3a/2 である。a = 2 寸のとき 4 円の和は 3 寸である。

2(res[6] + res[7] + res[8] + res[9]) |> simplify |> println

   3*a/2

算額の解答はここまで。

以下は,図を描くために丙円の中心座標 (x3, y3) を求める。
b, c, r3 がわかると,x3, y3 もわかる。

b = a/4
c = a/2
r3 = a*(-5 + sqrt(Sym(5)))/8 + sqrt(Sym(5))*a/4
eq10 = dist(a, 0, 0, c, x3, y3) - r3^2
eq11 = dist(a, 0, b, a, x3, y3) - r3^2
(x3, y3) = solve([eq10, eq11], (x3, y3))[2];

x3, y3 は短く簡約化できる。

@syms d
apart(x3) |> println

   a*(9/8 - 3*sqrt(5)/8)

y3 |> simplify |> sympy.sqrtdenest |> println

   a*(7 - sqrt(5))/8

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 2
   (b, c, ab, bc, ca, r1, r2, r3, r4) = a .* (1/4, 1/2, 5/4, √5/4, √5/2, 1/4, (3 - √5)/4, (√5 - 5)/8 + √5/4, (3 - √5)/8)
   (x3, y3) = a .* ((9 - 3√5)/8, (7 - √5)/8)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
   plot!([0, a, b, 0], [c, 0, a, c], color=:blue, lw=0.5)
   circle(a - r1, a - r1, r1)
   circle(r2, r2, r2, :green)
   circle(x3, y3, r3, :magenta)
   circle(r4, a - r4, r4, :purple)
   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(a, 0, " a", :green, :left, :bottom, delta=delta/2)
       point(b, a, " (b,a)", :green, :center, :bottom, delta=delta/2)
       point(0, c, "c ", :green, :right, :vcenter)
       point(a - r1, a - r1, "甲円:r1,(a-r1,a-r1)", :red, :center, delta=-delta/2)
       point(r2, r2, "乙円:r2,(r2,r2)", :green, :center, delta=-delta/2)
       point(x3, y3, "丙円:r3,(x3,y3)", :magenta, :center, delta=-delta/2)
       point(r4, a - r4, "丁円:r4\n(r4,a-r4)", :purple, :center, :bottom, delta=delta/2)
       plot!(xlims=(-4delta, a + 3delta))
   end
end;

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

算額(その1069)

2024年06月16日 | Julia

算額(その1069)

九十六 大船渡市立根 気仙安養寺稲荷堂 文政5年(1822)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

交差する2個の乙円の隙間に,甲円 3 個,乙円 4 個を容れる。丙円の直径が 5 寸のとき,乙円の直径はいかほどか。

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

include("julia-source.txt")

using SymPy

@syms r1::positive, r2::positive, x2::positive, y2::positive, r3::positive;
r1 = r3/2
eq1 = x2^2 + (r1 + y2)^2 - (r3 + r2)^2
eq2 = x2^2 + (2r1 - y2)^2 - (r1 + r2)^2
eq3 = x2^2 + (y2 - r1)^2 - (r3 - r2)^2
res = solve([eq1, eq2, eq3], (r2, x2, y2))[1]

   (3*r3/10, 2*sqrt(3)*r3/5, 3*r3/5)

乙円の半径 r2 は,丙円の半径 r3 の 3/10 倍である。
丙円の直径が 5 寸のとき,乙円の直径は 1.5 寸である。

「答」,「術」共に(山村も)「丙円径五寸」なのに,「乙円径九寸」と書いている。数値自体もおかしいが,図を見ても,乙円が丙円より大きいことはありえないことがわかるだろうに。

その他のパラメータは以下のとおりである。

   r3 = 2.5;  r1 = 1.25;  r2 = 0.75;  x2 = 1.73205;  y2 = 1.5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 5/2
   r1 = r3/2
   (r2, x2, y2) = (3*r3/10, 2*sqrt(3)*r3/5, 3*r3/5)
   @printf("丙円の直径が %g のとき,乙円の直径は %g である。\n", 2r3, 2r2)
   @printf("その他のパラメータは以下のとおりである。\n")
   @printf("r3 = %g;  r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", r3, r1, r2, x2, y2)
   plot()
   circle22(0, r1, r3)
   circle4(x2, y2, r2, :green)
   circle22(0, 2r1, r1, :blue)
   circle(0, 0, r1, :blue)
   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(0, 2r1, "甲円:r1,(0,2r1)", :blue, :center, :bottom, delta=delta)
       point(x2, y2, "乙円:r2\n(x2,y2)", :green, :center, :bottom, delta=delta)
       point(0, r1, "丙円:r3\n(0,r1)", :red, :center, :bottom, delta=delta)
   end
end;

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

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

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