算額（その799）

算額（その799）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

等脚台形の中に斜線を引き，区画された領域に甲円と乙円を入れる。乙円の直径が 15 寸，上頭（上底）の長さが 12 寸，斜線の長さが 20 寸のとき，甲円の直径はいかほどか。

等脚台形の高さを h，下頭（下底），上頭の長さをそれぞれ 2a, 2b
甲円の半径と中心座標を r1, (0, r1)
乙円の半径と中心座標を r2, (0, h - r2)
とおいて以下の連立方程式を解く。

なお，斜の長さについては算法助術の公式 39 を援用した。

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

using SymPy

@syms a::positive, b::positive, h::positive,
     r1::negative, r2::positive, 斜::positive, d
eq1 = h^2 + (a - b)^2 - (斜 + a + b)^2
eq2 = dist(a, 0, b, h, 0, r1) - r1^2
eq2 = numerator(apart(eq2, d))
eq3 = dist(a, 0, b, h, 0, h - r2) - r2^2
eq3 = numerator(apart(eq3, d))
res = solve([eq1, eq2, eq3], (r1, a, h));
res[2]

   (-b*斜/(2*r2) + r2 + r2*斜/(2*b), -斜/2 + r2^2/b + r2^2*斜/(2*b^2), r2*(2*b + 斜)/b)

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

b = 12//2
斜 = Sym(20)
r2 = 15//2
(-b*斜/(2*r2) + r2 + r2*斜/(2*b), -斜/2 + r2^2/b + r2^2*斜/(2*b^2), r2*(2*b + 斜)/b)

   (12, 15, 40)

甲円の半径の一般解は r2*(1 + 斜/2b) - b*斜/2r2 である

r2*(1 + 斜/2b) - b*斜/2r2 |> println

   12

図を描くために斜と等脚辺の交点座標 (x1, y1), (x2, y2) を求める。
一般解は求めにくいようなので，問の条件下で交点座標を求める。

@syms x1::positive, y1::positive, x2::negative, y2::positive, d
(a, h, r1) = (15, 40, 12)
eq4 = dist(x1, y1, x2, y2, 0, r1) - r1^2
eq4 = numerator(apart(eq4, d))
eq5 = dist(x1, y1, x2, y2, 0, h - r2) - r2^2
eq5 = numerator(apart(eq5, d))
eq6 = y1*(b - a) - h*(x1 - a)
eq7 = y2*(a - b) - h*(x2 + a)
res2 = solve([eq4, eq5, eq6, eq7], (x1, y1, x2, y2));
res2[2]

   (9*sqrt(10)/41 + 390/41, 1000/41 - 40*sqrt(10)/41, -390/41 + 9*sqrt(10)/41, 40*sqrt(10)/41 + 1000/41)

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

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = 12//2
   斜 = Sym(20)
   r2 = 15//2
   (a, h, r1) = (15, 40, 12)
   @printf("甲円の直径 = %g;  下頭 = %g;  高さ = %g\n", 2r1, 2a, h)
   plot([a, b, -b, -a, a], [0, h, h, 0, 0], color=:black, lw=0.5)
   circle(0, r1, r1)
   circle(0, h - r2, r2, :blue)
   if more == true
       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)
       (x1, y1, x2, y2) = (9*sqrt(10)/41 + 390/41, 1000/41 - 40*sqrt(10)/41, -390/41 + 9*sqrt(10)/41, 40*sqrt(10)/41 + 1000/41)
       segment(x1, y1, x2, y2, :green)
       point(b, h, "(b,h)", :black, :right, :bottom, delta=delta/2)
       point(a, 0, "(a,0) ", :black, :right, :bottom, delta=delta/2)
       point(x1, y1, "(x1,y1)", :black, :right, delta=-delta)
       point(x2, y2, "(x2,y2) ", :black, :left, delta=-2delta)
       point(0, r1, " 甲円:r1,(0,r1)", :red, :left, :vcenter)
       point(0, h - r2, " 乙円:r2\n (0,h-r2)", :blue, :left, :vcenter)
   end
end;

算額（その798）

算額（その798）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

半円の中に離矢を隔てて大円と小円が入っている。両円の直径の和が 4 寸，離矢が 2 寸 4 分のとき，外円（半円）の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (x1, r1)
小円の半径と中心座標を r2, (x2, r2)
2 円の直径の和を α，離矢を β
とおいて以下の連立方程式を解く。

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

using SymPy

@syms r1::positive, x1::positive, r2::positive,
     x2::negative, R::positive,
     α::positive, β::positive
eq1 = 2r1 + 2r2 - α
eq2 = x1^2 + r1^2 - (R - r1)^2
eq3 = x2^2 + r2^2 - (R - r2)^2
eq4 = (R^2 - (x1 - r1)^2) - β^2
eq5 = (x2 + r2) - (x1 - r1)
res = solve([eq1, eq2, eq3, eq4, eq5], (R, r1, x1, r2, x2));

2 通りの解が得られるが，大円と小円の位置関係だけが違うものである。
問の位置関係からいえば，2 番目のものが適解である。
この場合，外円の直径は α^2/(4(2β - α)) となる。
2 円の直径の和 α が 4 寸，離矢 β が 2.4 寸ならば，外円の直径は 5 寸である。

α = 4
β = 2.4
α^2/(4(2β - α))

   5.000000000000001

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

R = 2.5;  r1 = 1.2;  x1 = 0.5;  r2 = 0.8;  x2 = -1.5

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (α, β) = (4, 2.4)
   (R, r1, x1, r2, x2) = (5//2, 1.2, 0.5, 0.8, -1.5000000000000002)
   @printf("2円の直径の和が %g，離矢が %g のとき，外円の直径は %g 寸である。\n", α, β, 2R)
   @printf("R = %g;  r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g\n", R, r1, x1, r2, x2)
   plot()
   circle(0, 0, R, beginangle=0, endangle=180)
   circle(x1, r1, r1, :blue)
   circle(x2, r2, r2, :green)
   segment(x1 - r1, 0, x1 - r1, β, :magenta)
   if more == true
       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(x1, r1, "大円:r1,(x1,r1)", :blue, :center, delta=-delta)
       point(x2, r2, "小円:r2,(x2,r2)", :green, :center, delta=-delta)
       point(x1 - r1, β, "(x1-r1,β) ", :magenta, :right, :bottom, delta=delta)
       point(x1 - r1, 3β/4, "離矢 ", :magenta, :right, :vcenter, mark=false)
       point(R, 0, " R", :red, :left, :bottom, delta=delta)
   end
end;

算額（その797）

算額（その797）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

正方形の中に2本の斜線を引き，できた区画の中に大円，中円，小円を入れる。大円，中円，小円の直径および2本の斜線の長さと正方形の一辺の長さの和が 55 寸のとき，正方形の一辺の長さはいかほどか。

大円の半径と中心座標を r1, (r1, a - r1)
中円の半径と中心座標を r2, (a - r2, y2)
小円の半径と中心座標を r3, (r3, r3)
とおいて以下の連立方程式を解く。

SymPy の能力的に，解析解を求めることができないので，数値解を求める。

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

using SymPy

@syms r1::positive, r2::positive, y2::positive,
     r3::positive, a::positive, b::positive,
     c::positive, d::positive
eq1 = dist(b, 0, a, a, r1, a - r1) - r1^2
eq2 = dist(b, 0, a, a, a - r2, y2) - r2^2
eq3 = dist(b, 0, a, a, r3, r3) - r3^2
eq4 = dist(c,0, 0, d, r1, a - r1) - r1^2
eq5 = dist(c,0, 0, d, a - r2, y2) - r2^2
eq7 = c + d - sqrt(c^2 + d^2) - 2r3
eq8 = a + (a - b) - sqrt(a^2 + (a - b)^2) - 2r2
eq9 = 2r1 + 2r2 + 2r3 + sqrt(a^2 + (a - b)^2) + sqrt(c^2 + d^2) + a - 55;

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 Float64.(v), r.f_converged
end;

function H(u)
   (r1, r2, y2, r3, a, b, c, d) = u
   return [
       -r1^2 + (a - a*(a*(a - r1) + (a - b)*(-b + r1))/(a^2 + (a - b)^2) - r1)^2 + (-b + r1 - (a - b)*(a*(a - r1) + (a - b)*(-b + r1))/(a^2 + (a - b)^2))^2,  # eq1
       -r2^2 + (-a*(a*y2 + (a - b)*(a - b - r2))/(a^2 + (a - b)^2) + y2)^2 + (a - b - r2 - (a - b)*(a*y2 + (a - b)*(a - b - r2))/(a^2 + (a - b)^2))^2,  # eq2
       -r3^2 + (-a*(a*r3 + (a - b)*(-b + r3))/(a^2 + (a - b)^2) + r3)^2 + (-b + r3 - (a - b)*(a*r3 + (a - b)*(-b + r3))/(a^2 + (a - b)^2))^2,  # eq3
       -r1^2 + (a - d*(-c*(-c + r1) + d*(a - r1))/(c^2 + d^2) - r1)^2 + (-c + c*(-c*(-c + r1) + d*(a - r1))/(c^2 + d^2) + r1)^2,  # eq4
       -r2^2 + (-d*(-c*(a - c - r2) + d*y2)/(c^2 + d^2) + y2)^2 + (a - c + c*(-c*(a - c - r2) + d*y2)/(c^2 + d^2) - r2)^2,  # eq5
       c + d - 2*r3 - sqrt(c^2 + d^2),  # eq7
       2*a - b - 2*r2 - sqrt(a^2 + (a - b)^2),  # eq8
       a + 2*r1 + 2*r2 + 2*r3 + sqrt(a^2 + (a - b)^2) + sqrt(c^2 + d^2) - 55,  # eq9
   ]
end;

iniv = BigFloat[44, 30, 30, 21, 124, 35, 85, 62]./10
res = nls(H, ini=iniv)

   ([4.0, 3.0, 3.0, 2.0, 12.0, 3.0, 8.0, 6.0], true)

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

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

r1 = 4;  r2 = 3;  y2 = 3;  r3 = 2;  a = 12;  b = 3;  c = 8;  d = 6

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, y2, r3, a, b, c, d) = res[1]
   @printf("正方形の一辺の長さは %g 寸である。\n", a)
   @printf("r1 = %g;  r2 = %g;  y2 = %g;  r3 = %g;  a = %g;  b = %g;  c = %g;  d = %g\n", r1, r2, y2, r3, a, b, c, d)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
   circle(r1, a - r1, r1)
   circle(a - r2, y2, r2, :blue)
   circle(r3, r3, r3, :green)
   segment(c, 0, 0, d, :magenta)
   segment(b, 0, a, a, :magenta)
   if more == true
       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(r1, a - r1, "大円:r1,(r1,a-r1)", :red, :center, delta=-delta)
       point(a - r2, y2, "中円:r2,(a-r2,y2)", :blue, :center, delta=-delta)
       point(r3, r3, "小円:r3,(r3,r3)", :green, :center, delta=-delta)
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
       point(b, 0, "   b", :magenta, :left, :bottom, delta=delta/2)
       point(c, 0, "c   ", :magenta, :right, :bottom, delta=delta/2)
       point(0, d, "d ", :magenta, :right, :vcenter)
       plot!(xlims=(-0.5, 12.5))
   end
end;

算額（その796）

算額（その796）

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

直角三角形内に界斜（斜線）によって隔てられた領域に全円と小円が入っている。全円は三篇と界斜に接し，小円は二辺と界斜に接している。全円，小円の直径がそれぞれ 8 寸，2 寸のとき，界斜の長さはいかほどか。

全円の半径と中心座標を r1, (r1, r1)
小円の半径と中心座標を r2, (x2, r1)
斜線と直角三角形の辺との交点座標を (c, 0), (d, e)
とおいて以下の連立方程式を解く。

SymPy の能力的に，解析解を求めることができないので，数値解を求める。

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

using SymPy

@syms a::positive, b::positive, c::positive,
     d::positive, e::positive, r1::positive,
     r2::positive, x2::positive
#x2 = b - b*r2/r1 + r2
eq1 = (a + b - 2r1) - sqrt(a^2 + b^2)
eq4 = dist(0, a, b, 0, r1, r1) - r1^2
eq4 = numerator(apart(eq4, d))
eq5 = r2/(b - x2) - r1/(b - r1)
#eq5 = r2*(b - r1) - r1*(b - x2);
#eq5 = dist(0, a, b, 0, x2, r2) - r2^2
#eq5 = (sqrt(a^2 + b^2) + b + a)*r1 - a*b

eq2 = dist(c, 0, d, e, r1, r1) - r1^2
eq2 = numerator(apart(eq2, d))/e
eq3 = dist(c, 0, d, e, x2, r2) - r2^2
eq3 = numerator(apart(eq3, d))/e
eq6 = e*b - a*(b - d);

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 Float64.(v), r.f_converged
end;

function H(u)
   (a, b, c, d, e, x2) = u
   return [
       a + b - 2*r1 - sqrt(a^2 + b^2),  # eq1
       c^2*e - 2*c^2*r1 + 2*c*d*r1 - 2*c*e*r1 + 2*c*r1^2 - 2*d*r1^2,  # eq2
       c^2*e - 2*c^2*r2 + 2*c*d*r2 - 2*c*e*x2 + 2*c*r2*x2 - 2*d*r2*x2 - e*r2^2 + e*x2^2,  # eq3
       a^2*b^2 - 2*a^2*b*r1 - 2*a*b^2*r1 + 2*a*b*r1^2,  # eq4
       -r1/(b - r1) + r2/(b - x2),  # eq5
       -a*(b - d) + b*e,  # eq6
   ]
end;

(r1, r2) = (8, 2) .// 2
iniv = BigFloat[12, 16, 4.6, 13.2, 2.1, 13]
res = nls(H, ini=iniv)

   ([12.0, 16.0, 4.468871125850725, 13.22490309931942, 2.0813226755104353, 13.0], true)

c = 4.468871125850725, d = 13.22490309931942, e = 2.0813226755104353 となり，斜線の長さは 9 寸である。

c = 4.468871125850725
d = 13.22490309931942
e = 2.0813226755104353
sqrt((d - c)^2 + e^2)

   9.0

その他のパラメータは以下のとおりである。
a = 12;  b = 16;  c = 4.46887;  d = 13.2249;  e = 2.08132;  x2 = 13
なお，x2 = b - b*r2/r1 + r2 の関係がある。

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (8, 2) .// 2
   (a, b, c, d, e, x2) = res[1]
   # x2 = b - b*r2/r1 + r2
   l = sqrt((d - c)^2 + e^2)
   @printf("全円の直径が %g，小円の直径が %g のとき，界斜の長さは %g である。\n", 2r1, 2r2, l)
   @printf("a = %g;  b = %g;  c = %g;  d = %g;  e = %g;  x2 = %g\n", a, b, c, d, e, x2)
   plot([0, b, 0, 0], [0, 0, a, 0], color=:black, lw=0.5)
   circle(r1, r1, r1)
   circle(x2, r2, r2, :blue)
   segment(c, 0, d, e, :green)
   if more == true
       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(r1, r1, "全円:r1,(r1,r1)", :red, :center, delta=-delta)
       point(x2, r2, "全円:r2,(x2,r2) ", :black, :right, delta=-delta)
       point(0, a, " a", :black, :left, :bottom)
       point(b, 0, "b", :black, :left, :bottom, delta=delta)
       point(c, 0, "(c,0)", :green, :center, :bottom, delta=2delta)
       point(d, e, "(d,e)", :green, :center, :bottom, delta=2delta)
   end
end;

算額（その795）

算額（その795）

宮崎智也，吉岡倖佑，原本博史，安部利之：大西佐兵衛『雑題』の現代解について ~第8巻 第1問から第4問~，愛媛大学教育学部附属 科学教育研究センター紀要，Vol.2，pp.17-22，2023/03/31.
https://ed.ehime-u.ac.jp/CRESE/wp-content/uploads/2023/03/Vol2No3.pdf

円弧内に等円と小円を入れる。等円の直径と元の長さを用いて，分数の形で小円の直径を求めよ。

弦の長さを 2x0，弦と y 軸の交点座標を (0, y0)
円弧がその一部である円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r1, (r1, y0 + r1)
小円の半径と中心座標を r2, (x2, y0 + r2)
とおき，以下の連立方程式を解く。

たとえば，等円の直径が 18 寸，弦の長さが 54 寸のとき，小円の直径は 8 寸である。

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,
     x2::positive, x0::positive, y0::positive
eq1 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = r1^2 + (y0 + r1)^2 - (R - r1)^2
eq3 = x2^2 + (y0 + r2)^2 - (R - r2)^2
eq4 = x0^2 + y0^2 - R^2;
res = solve([eq1, eq2, eq3, eq4], (R, r2, x2, y0));

小円の半径は，r1, x0 を使って以下のように表すことができる。

r2 = res[1][2]
r2|> println

   r1*(-r1 + x0)^2*(r1 + x0)^2/(3*r1^2 + x0^2)^2

等円の直径が 18 寸，弦の長さが 54 寸のとき，小円の直径は 8 寸である。

r2(r1 => 18//2, x0 => 54//2) |> println

   4

その他のパラメータは以下のとおりである。
r1 = 9;  x0 = 27;  R = 28.125;  r2 = 4;  x2 = 21;  y0 = 7.875

SymPy の簡約化は中途半端なので，手動で簡約化する。

r22 = r1*(r1^2 - x0^2)^2 / (3r1^2 + x0^2)^2;
r2(r1 => 18//2, x0 => 54//2) |> println
r22(r1 => 18//2, x0 => 54//2) |> println

   4
   4

直径 d1 = 2r1, 弦の長さ x = 2x0 を使って等円の「直径」を求めるときは，次のようになる。

@syms d1, x
r23 = d1*(d1^2 - x^2)^2/(3*d1^2 + x^2)^2
r23(d1 => 18, x => 54) |> println

   8

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, x0) = (18, 54) .// 2
   (R, r2, x2, y0) = (225/8, 4, 21, 63/8)
   @printf("等円の直径が %g，弦の長さが %g のとき，小円の直径は %g\n", 2r1, 2x0, 2r2)
   @printf("r1 = %g;  x0 = %g;  R = %g;  r2 = %g;  x2 = %g;  y0 = %g\n", r1, x0, R, r2, x2, y0)
   θ = atand(y0, x0)
   plot()
   circle(0, 0, R, :black, beginangle=θ, endangle=180-θ, n=500)
   circle2(r1, y0 + r1, r1, :green)
   circle2(x2, y0 + r2, r2, :blue)
   segment(-x0, y0, x0, y0)
   segment(0, 0, x0, y0)
   if more == true
       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, R, " R", :black, :left, :bottom, delta=delta)
       point(0, y0, " y0", :black, :left, delta=-delta)
       point(x0, y0, "(x0,y0)", :black, :right, delta=-4delta)
       point(r1, y0 + r1, "等円:r1,(r1,y0+r1)", :green, :center, delta=-delta)
       point(x2, y0 + r2, "小円:r2,(x2,y0+r2)", :black, :right, delta=-delta)
   end
end;