裏 RjpWiki

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

算額(その581)

2023年12月22日 | Julia

算額(その581)

長崎市 鎮西大社諏訪神社 明治20年(1887)
米光丁: 長崎県の和算の概説

http://hyonemitsu.web.fc2.com/Nagasakiwasan.pdf

問題 4. 外円の中に大小の円を入れる。外円と大円の直径が 32 寸,18 寸のとき,小円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (x1, y - r1)
小円の半径と中心座標を r2,(x2, y - r2), (0, y + r2)
弦と y 軸の交点座標を (0, y)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, y::positive, r1::positive, x1::negative,
     r2::positive, x2::positive

eq1 = (x2 - x1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x1^2 + (y - r1)^2 - (R - r1)^2
eq3 = x2^2 + (y - r2)^2 - (R - r2)^2
eq4 = y + 2r2 - R
res = solve([eq1, eq2, eq3, eq4], (y, x1, r2, x2))

   1-element Vector{NTuple{4, Sym}}:
    ((-4*R^(11/2)*sqrt(r1) + 32*R^(9/2)*r1^(3/2) + 128*R^(7/2)*r1^(5/2) - 2048*R^(5/2)*r1^(7/2) + 7168*R^(3/2)*r1^(9/2) - 8192*sqrt(R)*r1^(11/2) + R^6 - 4*R^5*r1 - 64*R^4*r1^2 + 384*R^3*r1^3 + 256*R^2*r1^4 - 5120*R*r1^5 + 8192*r1^6)/(R^5 - 8*R^4*r1 - 32*R^3*r1^2 + 512*R^2*r1^3 - 1792*R*r1^4 + 2048*r1^5), 2*sqrt(2)*(-R^(11/2)*sqrt(r1) - R^(9/2)*r1^(3/2) + 68*R^(7/2)*r1^(5/2) - 80*R^(5/2)*r1^(7/2) - 1088*R^(3/2)*r1^(9/2) + 2560*sqrt(R)*r1^(11/2) + 4*R^5*r1 - 14*R^4*r1^2 - 200*R^3*r1^3 + 1184*R^2*r1^4 - 1408*R*r1^5 - 1024*r1^6)/(sqrt(R^(3/2)*sqrt(r1) + 8*sqrt(R)*r1^(3/2) - 5*R*r1 - 4*r1^2)*(-R^4 + 4*R^3*r1 + 48*R^2*r1^2 - 320*R*r1^3 + 512*r1^4)), 2*(R^(11/2)*sqrt(r1) - 8*R^(9/2)*r1^(3/2) - 32*R^(7/2)*r1^(5/2) + 512*R^(5/2)*r1^(7/2) - 1792*R^(3/2)*r1^(9/2) + 2048*sqrt(R)*r1^(11/2) - R^5*r1 + 8*R^4*r1^2 + 32*R^3*r1^3 - 512*R^2*r1^4 + 1792*R*r1^5 - 2048*r1^6)/(R^5 - 8*R^4*r1 - 32*R^3*r1^2 + 512*R^2*r1^3 - 1792*R*r1^4 + 2048*r1^5), sqrt(8*R^(3/2)*sqrt(r1) + 64*sqrt(R)*r1^(3/2) - 40*R*r1 - 32*r1^2))

小円の半径を表す式は長いが,問で与えられた数数値を代入すると,小円の半径は6,直径は 12 寸である。

res[1][3](R => 32//2, r1 => 18//2) |> println

   6

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

R = 16;  r1 = 9;  x1 = -4.89898;  r2 = 6;  x2 = 9.79796;  y = 4

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1) = (32, 18) .// 2
   (y, x1, r2, x2) = ((-4*R^(11/2)*sqrt(r1) + 32*R^(9/2)*r1^(3/2) + 128*R^(7/2)*r1^(5/2) - 2048*R^(5/2)*r1^(7/2) + 7168*R^(3/2)*r1^(9/2) - 8192*sqrt(R)*r1^(11/2) + R^6 - 4*R^5*r1 - 64*R^4*r1^2 + 384*R^3*r1^3 + 256*R^2*r1^4 - 5120*R*r1^5 + 8192*r1^6)/(R^5 - 8*R^4*r1 - 32*R^3*r1^2 + 512*R^2*r1^3 - 1792*R*r1^4 + 2048*r1^5), 2*sqrt(2)*(-R^(11/2)*sqrt(r1) - R^(9/2)*r1^(3/2) + 68*R^(7/2)*r1^(5/2) - 80*R^(5/2)*r1^(7/2) - 1088*R^(3/2)*r1^(9/2) + 2560*sqrt(R)*r1^(11/2) + 4*R^5*r1 - 14*R^4*r1^2 - 200*R^3*r1^3 + 1184*R^2*r1^4 - 1408*R*r1^5 - 1024*r1^6)/(sqrt(R^(3/2)*sqrt(r1) + 8*sqrt(R)*r1^(3/2) - 5*R*r1 - 4*r1^2)*(-R^4 + 4*R^3*r1 + 48*R^2*r1^2 - 320*R*r1^3 + 512*r1^4)), 2*(R^(11/2)*sqrt(r1) - 8*R^(9/2)*r1^(3/2) - 32*R^(7/2)*r1^(5/2) + 512*R^(5/2)*r1^(7/2) - 1792*R^(3/2)*r1^(9/2) + 2048*sqrt(R)*r1^(11/2) - R^5*r1 + 8*R^4*r1^2 + 32*R^3*r1^3 - 512*R^2*r1^4 + 1792*R*r1^5 - 2048*r1^6)/(R^5 - 8*R^4*r1 - 32*R^3*r1^2 + 512*R^2*r1^3 - 1792*R*r1^4 + 2048*r1^5), sqrt(8*R^(3/2)*sqrt(r1) + 64*sqrt(R)*r1^(3/2) - 40*R*r1 - 32*r1^2))
   @printf("小円の直径 = %g;  R = %g;  r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g;  y = %g\n",
      2r2, R, r1, x1, r2, x2, y)
   plot()
   circle(0, 0, R, :magenta)
   circle(x1, y - r1, r1)
   circle(x2, y - r2, r2, :blue)
   circle(0, y + r2, r2, :blue)
   x = sqrt(R^2 - y^2)
   segment(-x, y, x, y)
   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(x1, y - r1, "大円:r1,(x1,y-r1)", :red, :center, :top, delta=-delta/2)
       point(x2, y - r2, "小円:r2,(x2,y-r2)", :blue, :center, :top, delta=-delta/2)
       point(0, y + r2, " 小円:r2\n (x2,y+r2)", :blue, :left, :vcenter)
       point(0, y, " y", :black, :left, :top, delta=-delta/2)
       point(0, R, " R", :magenta, :left, :bottom, delta=delta/2)
   end
end;

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

算額(その580)

2023年12月22日 | Julia

算額(その580)

群馬の算額 19-4 榛名神社 文化8年
http://takasakiwasan.web.fc2.com/gunnsann/g019-4.html

鈎股弦(直角三角形)の中に正方形,大円,小円を入れる。中鈎とは直角の頂点から弦へおろした垂線である。大円,小円の直径が 12 寸,5 寸のとき,正方形の一辺の長さを求めよ。

「鈎」,「股」の長さを,鈎,股とする。
中鈎と弦の交点座標を (x, y) とする。
正方形の一辺の長さを a とする。
大円の半径と中心座標を r1, (股 - a + r1, r1)
小円の半径と中心座標を r2, (股 - r2, a - r2)
として,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, a::positive, r1::positive, x1::positive,
     r2::positive, x::positive, y::positive
(r1, r2) = (12, 5) .// 2
eq1 = a/(股 - a) - 鈎/股
eq2 = distance(股, 0, x, y, 股 - a + r1, r1) - r1^2
eq3 = distance(股, 0, x, y, 股 - r2, a - r2) - r2^2
eq4 = 股*y + 鈎*(股 - x) - 股*鈎
eq5 = a*(鈎 - a) + a*(股 - a) + 2*a^2 - 股*鈎;
# res = solve([eq1, eq2, eq3, eq4, eq5], (鈎, 股, a, x, y))

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=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (鈎, 股, a, x, y) = u
   return [
       a/(-a + 股) - 鈎/股,  # eq1
       (-y*(-a*x + a*股 + 6*x + 6*y - 6*股)/(x^2 - 2*x*股 + y^2 + 股^2) + 6)^2 + (-a + 股 - (股*(x^2 - 2*x*股 + y^2 + 股^2) + (x - 股)*(-a*x + a*股 + 6*x + 6*y - 6*股))/(x^2 - 2*x*股 + y^2 + 股^2) + 6)^2 - 36,  # eq2
       (a - y*(2*a*y - 5*x - 5*y + 5*股)/(2*(x^2 - 2*x*股 + y^2 + 股^2)) - 5/2)^2 + (股 - (股*(x^2 - 2*x*股 + y^2 + 股^2) + (x - 股)*(2*a*y - 5*x - 5*y + 5*股)/2)/(x^2 - 2*x*股 + y^2 + 股^2) - 5/2)^2 - 25/4,  # eq3
       y*股 - 股*鈎 + 鈎*(-x + 股),  # eq4
       2*a^2 + a*(-a + 股) + a*(-a + 鈎) - 股*鈎,  # eq5
   ]
end;

(r1, r2) = (12, 5) .// 2
iniv = BigFloat[20.6, 50.1, 14.5, 42.1, 17.9]
res = nls(H, ini=iniv)

   (BigFloat[21.35733319823884343090635521344017153294802876982109796830915235365518705818788, 50.39219873866784319300311036791102016301992183128791550477294030211998940687325, 15.00000000000000000000000000000000000000000000000000000000000000000000000000028, 42.82892754606824024613300750612840676283894513302086625119773230436335822945181, 18.15185086223904707248824686827827216043434407584091820858049919461591482581305], true)

正方形の一辺の長さは 15 寸である。

他のパラメータは以下の通り。
鈎 = 21.3573;  股 = 50.3922;  a = 15;  x1 = 41.3922;  x = 42.8289;  y = 18.1519

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (12, 5) .// 2
   (鈎, 股, a, x1, x, y) = [98, 238, 69, 195, 200, 85] .* (12/57)
   (鈎, 股, a, x, y) = res[1]
   @printf("鈎 = %g;  股 = %g;  a = %g;  x1 = %g;  x = %g;  y = %g\n",
       鈎, 股, a, x1, x, y)
   plot([0, 股, 股, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(股 - a + r1, r1, r1)
   circle(股 - r2, a - r2, r2, :blue)
   segment(股, 0, x, y)
   rect(股 - a, 0, 股, a, :green)
   plot!(xlims=(-1, 57), ylims=(-2.5, 25))
   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 + r1, r1, "大円:r1\n(股-a+r1,r1)", :red, :center, :bottom, delta=delta)
       point(股 - r2, a - r2, "小円:r2,(股-r2,a-r2) ", :blue, :right, :vcenter)
       point(股, a, " (股,a)", :green, :left, :vcenter)
       point(x, y, "(x,y)", :black, :right, :bottom, delta=delta)
       point(股 - a, 0, "股-a", :green, :center, :top, delta=-delta/2)
       point(股, 0, "股", :green, :center, :top, delta=-delta/2)
       point(股, 鈎, "(股,鈎) ", :black, :right, :bottom, delta=delta/2)
   end
end;

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

算額(その579)

2023年12月22日 | Julia

算額(その579)

群馬の算額 19-6 榛名神社 文化8年
http://takasakiwasan.web.fc2.com/gunnsann/g019-6.html

3 本の線分で 3 個の円を挟んでいる。大円,小円の直径が 4 寸,1 寸 4 分 4 厘 のとき,界斜(図の記号では点 (x,0) と (x0, y0) を結ぶ線分)の長さを求めよ。

界斜の両端の座標を (x,0) と (x0, y0) とする。
大円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r1, (x11, r1), (x12, r1)
として,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, x11::positive, x12::positive, r2::positive, x2::positive,
     x::positive, x0::positive, y0::positive
(r1, r2) = (144//200, 4//2)
eq1 = (x2 - x12)^2 + (r2 - r1)^2 - (r1 + r2)^2
eq2 = distance(0, 0, x0, y0, x11, r1) - r1^2
eq3 = distance(0, 0, x0, y0, x2, r2) - r2^2    
eq4 = distance(x, 0, x0, y0, x11, r1) - r1^2
eq5 = distance(x, 0, x0, y0, x12, r1) - r1^2
eq6 = distance(x, 0, x0, y0, x2, r2) - r2^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (x11, x12, x2, x, x0, y0))

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=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (x11, x12, x2, x, x0, y0) = u
   return [
       (-x12 + x2)^2 - 144/25,  # eq1
       (-x0*(25*x0*x11 + 18*y0)/(25*(x0^2 + y0^2)) + x11)^2 + (-y0*(25*x0*x11 + 18*y0)/(25*(x0^2 + y0^2)) + 18/25)^2 - 324/625,  # eq2
       (-x0*(x0*x2 + 2*y0)/(x0^2 + y0^2) + x2)^2 + (-y0*(x0*x2 + 2*y0)/(x0^2 + y0^2) + 2)^2 - 4,  # eq3
       (x11 - (x11*(x^2 - 2*x*x0 + x0^2 + y0^2) + y0*(25*x*y0 - 18*x + 18*x0 - 25*x11*y0)/25)/(x^2 - 2*x*x0 + x0^2 + y0^2))^2 + (-y0*(25*x^2 - 25*x*x0 - 25*x*x11 + 25*x0*x11 + 18*y0)/(25*(x^2 - 2*x*x0 + x0^2 + y0^2)) + 18/25)^2 - 324/625,  # eq4
       (x12 - (x^2*x12 - 2*x*x0*x12 + x*y0^2 - 18*x*y0/25 + x0^2*x12 + 18*x0*y0/25)/(x^2 - 2*x*x0 + x0^2 + y0^2))^2 + (-y0*(25*x^2 - 25*x*x0 - 25*x*x12 + 25*x0*x12 + 18*y0)/(25*(x^2 - 2*x*x0 + x0^2 + y0^2)) + 18/25)^2 - 324/625,  # eq5
       (x2 - (x^2*x2 - 2*x*x0*x2 + x*y0^2 - 2*x*y0 + x0^2*x2 + 2*x0*y0)/(x^2 - 2*x*x0 + x0^2 + y0^2))^2 + (-y0*(x^2 - x*x0 - x*x2 + x0*x2 + 2*y0)/(x^2 - 2*x*x0 + x0^2 + y0^2) + 2)^2 - 4,  # eq6
   ]
end;

(r1, r2) = (144//200, 4//2)
iniv = BigFloat[2.2, 4, 6.4, 2.6, 5.412, 3.8]
res = nls(H, ini=iniv)

   (BigFloat[2.32537499999999981452205370806705695806346869674534523102666679819954313421592, 4.05937499999999973306648353848901437850698923197181920887229922564784525145764, 6.459374999999999688657562553482752350696667705728683787732940545284447986400679, 2.709374999999999822647282592843187726735573892402228819331695263097710712289562, 5.012399221453286844773215617274192837885357156437953123416192597373948544218003, 3.433079584775086464352099968927128171012730897483676505879678311456937136579126], true)

界斜の長さは sqrt((x0 - x)^2 + y0^2) ≒ 4.134である。

x  = 2.7093749999999998
x0 = 5.0123992214532868
y0 = 3.4330795847750864
sqrt((x0 - x)^2 + y0^2)

   4.134

他のパラメータは以下の通り。
r1 = 0.72;  r2 = 2;  x11 = 2.32537;  x12 = 4.05937;  x2 = 6.45937

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (144//200, 4//2)
   (x11, x12, x2, x, x0, y0) = res[1]
   l = sqrt((x0 - x)^2 + y0^2)
   @printf("界斜の長さ = %g;  r1 = %g;  r2 = %g;  x11 = %g;  x12 = %g;  x2 = %g;  x = %g;  x0 = %g;  y0 = %g\n",
       l, r1, r2, x11, x12, x2, x, x0, y0)
   plot()
   circle(x11, r1, r1)
   circle(x12, r1, r1)
   circle(x2, r2, r2, :blue)
   abline(0, 0, y0/x0, 0, 6.5)
   segment(x, 0, x0, y0)
   plot!(ylims=(-0.3, 4.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(x11, r1, "小円:r1\n(x11,r2)", :red, :center, :bottom, delta=delta)
       point(x12, r1, "小円:r1\n(x12,r2)", :red, :center, :bottom, delta=delta)
       point(x2, r2, "大円:r2,(x2,r2)", :magenta, :center, :top, delta=-delta)
       point(x, 0, "x", :black, :center, :top, delta=-delta/2)
       point(x0, y0, "(x0,y0)", :black, :right, :bottom, delta=delta/2)
   end
end;

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

算額(その578)

2023年12月22日 | Julia

算額(その578)

群馬の算額 17-5 八幡宮
http://takasakiwasan.web.fc2.com/gunnsann/g017-5.html

直線,垂線と斜線に大円,中円,小円が接している。大円,中円の直径が 5 寸,3 寸のとき,小円の直径を求めよ。

斜線の x 切片,y 切片を x, y とする。
大円の半径と中心座標を r1, (r1, r1)
中円の半径と中心座標を r2, (-r2, r2)
小円の半径と中心座標を r3, (r3, y3)
として,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive, x::positive, y::positive

eq1 = distance(x, 0, 0, y, r1, r1) - r1^2
eq2 = distance(x, 0, 0, y, -r2, r2) - r2^2
eq3 = distance(x, 0, 0, y, r3, r3) - r3^2
res = solve([eq1, eq2, eq3], (r3, x, y))

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

4 番目の解が適解である。

小円の直径は r2*(r1 - r2)/(r1 + r2) で,r1 = 5//2, r2 = 3//2 のときは,3/4 = 0.75 = 7 分 5 厘である。

その他のパラメータは以下の通り。
小円の直径 = 0.75;  r1 = 2.5;  r2 = 1.5;  r3 = 0.375;  x = 1;  y = 1.875

(r1, r2) = (5, 3) .// 2
r2*(r1 - r2)/(r1 + r2)

   3//8

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (5, 3) .// 2
   (r3, x, y) = (r2*(r1 - r2)/(r1 + r2), r1 - r2, 2*r1*r2/(r1 + r2))
   @printf("小円の直径 = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x = %g;  y = %g\n", 2r3, r1, r2, r3, x, y)
   plot()
   circle(r1, r1, r1)
   circle(-r2, r2, r2, :blue)
   circle(r3, r3, r3, :magenta)
   abline(0, y, -y/x, -5, 5)
   plot!(xlims=(-3.5, 5.2), ylims=(-0.5, 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(r1, r1, "大円:r1,(r1,r1)", :red, :center, :top, delta=-delta)
       point(-r2, r2, "中円:r2,(-r2,r2)", :blue, :center, :top, delta=-delta)
       point(r3, r3, "  小円:r3,(r3,r3)", :magenta, :left, :bottom, delta=5delta)
       point(x, 0, " x", :black, :left, :bottom, delta=delta/2)
       point(0, y, "  y", :black, :left, :vcenter)
   end
end;

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

算額(その577)

2023年12月22日 | Julia

算額(その577)

群馬の算額 18-5 八幡宮
http://takasakiwasan.web.fc2.com/gunnsann/g018-5.html

2 本の直線と 1 本の垂線に接している甲円,乙円,丙円がある。乙円,丙円の直径が 12 寸,3 寸のとき,甲円の直径を求めよ。

甲円の半径と中心座標を r1, (-r1, r1)
乙円の半径と中心座標を r2, (r2, r2)
丙円の半径と中心座標を r3, (r3, y3)
とする。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive, y3::positive,
     a::negative, b::positive, x::positive, y::positive

(r2, r3) = (12, 3) .// 2

   (6//1, 3//2)

乙円と丙円が外接することから,丙円の中心座標 y3 が求まる。

eq1 = (r2 - r3)^2 + (y3 - r2)^2 - (r2 + r3)^2
y3 = solve(eq1)[1]
y3 |> println

   12

乙円と丙円の中心を通る直線の傾き a と y 切片 b を求める。

eq2 = r3*a + b - y3
eq3 = r2*a + b - r2
res = solve([eq2, eq3], (a, b))

   Dict{Any, Any} with 2 entries:
     b => 14
     a => -4/3

甲円と丙円の共通接線と乙円と丙円の中心を通る直線の y 切片は同じ 14 である。
甲円と乙円の中心座標と共通接線の距離についての式を立てる。共通接線は (0, 14) と (x, y) を通る。x は任意で(
例えば 10 とか), y は x に応じて決まる。

eq5 = distance(0, 14, 10, y, -r1, r1) - r1^2
eq6 = distance(0, 14, 10, y, r2, r2) - r2^2
res = solve([eq5, eq6])

   1-element Vector{Dict{Any, Any}}:
    Dict(r1 => 8, y => 133/12)

甲円の半径は 8,直径は 16 寸である。

共通接線は (0, 14) と (10, 133/12) を通る。傾きは -7//24

(133//12 - 14)//10 

   -7//24

他のパラメータは以下の通り。

   甲円の直径 = 16;  r1 = 8;  r2 = 6;  r3 = 1.5;  y3 = 12;  y 切片 = 14;  傾き = -0.291667

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3) = (12, 3) .// 2
   y3 = 12
   r1 = 8
   intercept = 14
   slope = -7//24
   @printf("甲円の直径 = %g;  r1 = %g;  r2 = %g;  r3 = %g;  y3 = %g;  y 切片 = %g;  傾き = %g\n", 2r1, r1, r2, r3, y3, intercept, slope)
   plot()
   circle(r2, r2, r2)
   circle(r3, y3, r3, :blue)
   circle(-r1, r1, r1, :magenta)
   abline(0, intercept, slope, -18, 15)
   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(-r1, r1, "甲円:r1,(-r1,r1)", :magenta, :center, :top, delta=-delta)
       point(r2, r2, "乙円:r2,(r2,r2)", :red, :center, :top, delta=-delta)
       point(r3, y3, "丙円:r3,(r3,y3)", :blue, :left, :bottom, delta=7delta)
       point(0, 14, "14 ", :black, :right, :top)
       point(10, 133/12, "(10,133/12)", :black, :left, :bottom, delta=2delta)
   end
end;

 

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

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

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