算額（その691）

算額（その691）

神壁算法 東都本郷春木町一丁目 佐野富次郎安行 天明七年
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

楕円の中に大円 2 個，中円 2 個，小円 1 個が入っている。
楕円の長径，短径がそれぞれ 2168 寸，1084 寸のとき，大円の直径はいかほどか。

楕円の長半径と短半径を a, b
大円の半径と中心座標を r1, (x1, 0); x1 = r1 + r3
中円の半径と中心座標を r2, (0, b - r2)
小円の半径と中心座標を r3, (0, 0)
楕円と大円の共通接点の座標を (x0, y0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms a::positive, b::positive, x0::positive, y0::positive,
     r1::positive, x1::positive,r2::positive, r3::positive
(a, b) = (2168, 1084) .// 2
x1 = r1 + r3
eq1 = (x0 - x1)^2 + y0^2 - r1^2
eq2 = -b^2*x0/(a^2*y0) - (x1 - x0)/y0
eq3 = x0^2/a^2 + y0^2/b^2 - 1
eq4 = r3 + 2r2 - b
eq5 = x1^2 + (b - r2)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (x0, y0, r1, r2, r3))

   1-element Vector{NTuple{5, Sym}}:
    (-1626/7 + 2710*sqrt(57)/21, sqrt(6315926/147 + 734410*sqrt(57)/49), 271*sqrt(57)/14 + 4065/14, 3523/7 - 271*sqrt(57)/7, -3252/7 + 542*sqrt(57)/7)

大円の直径は (271√57 + 4065)/7 = 873.000733136910 寸である。
ちなみに，中円，小円の直径はそれぞれ 421.998533726179 寸，240.002932547642 寸である。

res[1][3].evalf()*2 |> println
res[1][4].evalf()*2 |> println
res[1][5].evalf()*2 |> println

   873.000733136910
   421.998533726179
   240.002932547642

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (2168, 1084) .// 2
   (x0, y0, r1, r2, r3) = (-1626/7 + 2710*sqrt(57)/21, sqrt(6315926/147 + 734410*sqrt(57)/49), 271*sqrt(57)/14 + 4065/14, 3523/7 - 271*sqrt(57)/7, -3252/7 + 542*sqrt(57)/7)
   x1 = r1 + r3
   plot()
   ellipse(0, 0, a, b, color=:orange)
   circle(x1, 0, r1)
   circle(-x1, 0, r1)
   circle(0, 0, r3, :green)
   circle(0, b - r2, r2, :blue)
   circle(0, r2 - b, r2, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(a, 0, " a", :black, :left, :bottom, delta=delta)
       point(0, b, " b", :black, :left, :bottom, delta=delta)
       point(x0, y0, "(x0,y0)", :red, :left, :bottom, delta=delta)
       point(x1, 0, "大円:r1,(x1,y1)", :red, :center, delta=-delta)
       point(0, b - r2, " 中円:r2,(0,b-r2)", :black, :left, :vcenter)
       point(0, r3, " r3", :black, :left, :bottom, delta=delta)
       point(0, r3, " 小円", :black, :left, :top, delta=-delta)
   end
end;

 

算額（その690）

算額（その690）

神壁算法 一柳土佐守家士 高瀬恆右衛門信之 寛政八年
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

外円の中に大円 5 個，小円 10 個が入っている。外円の直径が 33 寸 99 分のとき，小円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, R - r1)
小円の半径と中心座標を r2, (x1, y1), (x2, y2)
とおき，以下の連立方程式を解く。

単なる二元連立方程式であるが，SymPy の性能では，解析解を求めることができない。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms R::positive, r1::positive, r2::positive,
     x1::positive, y1::positive, x2::positive, y2::positive

R = 33.99/2
x1 = (R - r2)*cosd(54)
y1 = (R - r2)*sind(54)
x2 = (R - r1)*cosd(36)*cosd(54)
y2 = (R - r1)*cosd(36)*sind(54)
eq1 = x1^2 + (R - r1 - y1)^2 - (r1 + r2)^2
eq2 = x2^2 + (R - r1 - y2)^2 - (r1 - r2)^2
res = solve([eq1, eq2], (r1, r2))

   2-element Vector{Tuple{Sym, Sym}}:
    (0.972560942096724, 10.3903143260872)
    (7.45655928660198, 1.85000450540052)

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

小円の直径は 2\*1.85000450540052 = 3.70000901080104
大円の直径は 2\*7.45655928660198 = 14.91311857320396

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

   R = 16.995;  r1 = 7.45656;  r2 = 1.85;  x1 = 8.902;  y1 = 12.2526;  x2 = 4.5358;  y2 = 6.24299

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 33.99/2
   (r1, r2) = (7.45655928660198, 1.85000450540052)
   x1 = (R - r2)*cosd(54)
   y1 = (R - r2)*sind(54)
   x2 = (R - r1)*cosd(36)*cosd(54)
   y2 = (R - r1)*cosd(36)*sind(54)
   @printf("小円の直径 = %g\n", 2r2)
   @printf("R = %g;  r1 = %g;  r2 = %g;  x1 = %g;  y1 = %g;  x2 = %g;  y2 = %g\n",
       R, r1, r2, x1, y1, x2, y2)
   plot()
   circle(0, 0, R, :orange)
   rotate(0, R - r1, r1, angle=72)
   rotate(x1, y1, r2, :blue, angle=72)
   rotate(x2, y2, r2, :blue, angle=72)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(0, R - r1, " 大円:r1\n (0,R-r1)", :red, :left, :bottom, delta=delta)
       point(x1, y1, " 小円:r2,(x1,y1)", :black, :left, :vcenter)
       point(x2, y2, " 小円:r2,(x2,y2)", :black, :left, :vcenter)
   end
end;

 

算額（その689）

算額（その689）

和算問題あれこれ 2 令和5年4月の問題-No.2（『神壁算法』15問）
https://gunmawasan.web.fc2.com/k-n-mondai.html

大円内に等脚台形と小円をいれる。台形の上底（上頭），下底（下頭），対角線（内斜）がそれぞれ 36 寸，300 寸，280 寸のとき，小円の直径は何寸か。
注：「神壁算法」の図の寸法を10倍している。

 

大円の半径と中心座標を R, (0, 0)
小円の半径と中心座標を r0, (x0, y0)
上底，下底と円の交点座標を (x1, y1), (x2, y2)
とおき，以下の連立方程式を解く。

なお，4 元連立方程式として一度に解くことはできるが，R, r0 が負の値になる（絶対値を取ればもんだいはない）ので，eq1, eq2 の 2 元連立方程式として R, r0 を求め，ついで eq3, eq4 で x0，y0 を求める。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms R::positive, x1::positive, y1::positive, x2::positive, y2::negative, r0::positive, x0::positive, y0::positive

(x1, x2) = (300, 36) .// 2
y1 = -sqrt(R^2 - x1^2)
y2 = sqrt(R^2 - x2^2)
eq1 = (x1 + x2)^2 + (y2 - y1)^2 - 280^2
eq2 = 2sqrt(R^2 - (R - 2r0)^2) - sqrt((x1 - x2)^2 + (y2 - y1)^2);
solve([eq1, eq2], (R, r0))

   2-element Vector{Tuple{Sym, Sym}}:
    (325/2, 65/2)
    (325/2, 130)

最初のものが適解である。

小円の直径は 65 寸である。ちなみに，外円の直径は 325 寸である。

図を描くために，小円の中心座標 (x0, y0) を求める。

(R, r0) = (325//2, 65//2)
eq3 = x0^2 + y0^2 - (R - r0)^2
eq4 = y0/x0 * (sqrt(R^2 - x2^2) + sqrt(R^2 - x1^2))/(x2 - x1) + 1;
solve([eq3, eq4], (x0, y0))

   1-element Vector{Tuple{Sym, Sym}}:
    (112.000000000000, 66.0000000000000)

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x1, x2) = (300, 36) .// 2
   (R, r0, x0, y0) = [325/2, 65/2, 112, 66]
   y1 = -sqrt(R^2 - x1^2)
   y2 = sqrt(R^2 - x2^2)
   @printf("小円の直径 = %g\n", 2r0)
   @printf("R = %g;  x0 = %g;  y0 = %g;  r0 = %g\n", R, x0, y0, r0)
   plot([x1, x2, -x2, -x1, x1], [y1, y2, y2, y1, y1], color=:black, lw=0.5)
   circle(0, 0, R)
   circle(x0, y0, r0, :blue)
   segment(-x2, y2, x1, y1, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(x0, y0, "小円:r0\n(x0,y0)", :blue, :center, delta=-delta/2)
       point(R, 0, " R", :red, :left, :bottom, delta=delta)
       point(x1, y1, "(x1,y1)")
       point(x2, y2, "(x2,y2)", :green, :left, :bottom, delta=delta)
   end
end;