算額（その688）

算額（その688）

和算問題あれこれ 2 令和5年10月の問題-No.3（『五明算法』34問）
https://gunmawasan.web.fc2.com/k-n-mondai.html

扇面に 6 この円が入っている，甲円，乙円の直径が 1.75 寸，2.25 寸のとき，丁円の直径はいかほどか。

扇の長さ（半径）と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (x3, y3); r3 = r1 + r4
丁円の半径と中心座標を r4, (0, R - 2r1 - r4)
とおき，以下の連立方程式を解く。

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

using SymPy

@syms R, r1, r2, x2::positive, y2, r3, x3, y3, r4, x0, y0
(r1, r2) = (175, 225) .// 200
r3 = r1 + r4
eq1 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
eq2 = (x3 - x2)^2 + (y3 - y2)^2 - (r3 + r2)^2
eq3 = x3^2 + (R -2r1 - r4 - y3)^2 - (r3 + r4)^2
eq4 = x3^2 + y3^2 - (R - r3)^2
eq5 = x2^2 + (y2 -R + 2r1 + r4)^2 - (r2 + r4)^2
eq6 = x2^2 + y2^2 - (R - r2)^2
res1 = solve([eq1, eq2, eq3, eq4, eq5, eq6], (R, x2, y2, x3, y3, r4))

   1-element Vector{NTuple{6, Sym}}:
    (63/8, 9*sqrt(47)/32, 207/32, 189*sqrt(47)/320, 5607/1600, 329/200)

丁円の直径は 329/100 = 3.29 寸である。

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

   r1 = 0.875;  r2 = 1.125;  R = 7.875;  x2 = 1.92815;  y2 = 6.46875;  x3 = 4.04912;  y3 = 3.50438;  r4 = 1.645;  x0 = 7.67922;  y0 = 1.74504

問の答えはここまででよいが，図を描くために追加の計算をする。
扇の角の座標 (x0, y0) を求める。

@syms R, r1, r2, x2::positive, y2, r3, x3, y3, r4, x0::positive, y0::positive
(r1, r2) = (175, 225) .// 200
(R, x2, y2, x3, y3, r4) = res1[1]
r3 = r1 + r4
eq11 = x0^2 + y0^2 - R^2
eq12 = distance(0, 0, x0, y0, x3, y3) - r3^2
res2 = solve([eq11, eq12], (x0, y0))

   2-element Vector{Tuple{Sym, Sym}}:
    (-5607/2312 + 14175*sqrt(47)/18496, 945*sqrt(47)/2312 + 84105/18496)
    (5607/2312 + 14175*sqrt(47)/18496, 84105/18496 - 945*sqrt(47)/2312)

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

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (175, 225) .// 200
   (R, x2, y2, x3, y3, r4) = (63/8, 9*sqrt(47)/32, 207/32, 189*sqrt(47)/320, 5607/1600, 329/200)
   r3 = r1 + r4
   (x0, y0) = (5607/2312 + 14175*sqrt(47)/18496, 84105/18496 - 945*sqrt(47)/2312)
   θ = atand(y0/x0)
   println("丁円の直径 = $(2r4)")
   @printf("r1 = %g;  r2 = %g;  R = %g;  x2 = %g;  y2 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  x0 = %g;  y0 = %g\n", r1, r2, R, x2, y2, x3, y3, r4, x0, y0)
   plot()
   circle(0, 0, R, beginangle=θ, endangle=180-θ, n=1500)
   circle(0, 0, R - 2r3, :orange, beginangle=θ, endangle=180-θ, n=1500)
   segment(0, 0, x0, y0)
   segment(0, 0, -x0, y0)
   circle(0, R - r1, r1, :blue)
   circle(x2, y2, r2, :magenta)
   circle(-x2, y2, r2, :magenta)
   circle(x3, y3, r3, :green)
   circle(-x3, y3, r3, :green)
   circle(0, R - 2r1 - r4, r4, :tomato)
   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,(0,R-r1)", :black, :left, :vcenter)
       point(x2, y2, " 乙円:r2,(x2,y2)", :black, :left, :vcenter)
       point(x3, y3, " 丙円:r3\n (x3,y3)", :green, :left, :vcenter)
       point(0, R-2r1-r4, " 丁円:r4\n (0,R-2r1-r4)", :black, :left, :vcenter)
       point(0, R, " R", :black, :left, :bottom, delta=delta)
       point(0, R - 2r3, " R-2r3", :black, :left, delta=-delta)
       point(x0, y0, "(x0,y0) ", :red, :right, :bottom, delta=delta/2)
   end
end;

算額（その687）

算額（その687）

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

台形内に左斜と右斜の2本の線分を設ける。左斜，右斜，闊（高さ）の長さがそれぞれ 15 寸，5 寸，12 寸である。甲，乙を変化させたとき，台形の面積が最大になるときの甲，乙の長さを求めよ。

横倒しになっている台形であるが，下底（左側），下底（右側）の長さはそれぞれ固定されているので，甲によって増減する。したがって，上底，下底，闊により決まる面積も甲の関数になる。
面積(甲) = 闊*(sqrt(右斜^2 - 甲^2) + sqrt(左斜^2 - (-甲 + 闊)^2))/2
左斜，右斜，闊 がそれぞれ 15, 5, 12 ならば
面積(甲) = 6*sqrt(25 - 甲^2) + 6*sqrt(225 - (12 - 甲)^2)
である。甲を横軸，面積を縦軸に取り図を描いてみる。

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

using SymPy

@syms 甲, 左斜, 右斜, 闊, 下底, 上底, 面積
(左斜, 右斜, 闊) = ( 15, 5, 12)
下底 = sqrt(左斜^2 -  (闊 - 甲)^2)
上底 = sqrt(右斜^2 -  甲^2)
面積 = (下底 + 上底)*闊/2
面積 |> println

   6*sqrt(25 - 甲^2) + 6*sqrt(225 - (12 - 甲)^2)

pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(面積, xlims=(0, 5), xlabel="甲", ylabel="面積")

図より，甲が3前後のときに面積が最大になることがわかる。

面積が最大になるときの甲は，面積の導関数を求めそれが 0 になるときの値である。方程式を解くことにより，甲 = 3 である。また，このとき面積を表す関数の甲にその値を代入すれば，最大値が 96 であることがわかる。

solve(diff(面積))[1] |> println

   3

面積(甲 => 3) |> println

   96

 

 

算額（その686）

算額（その686）

三七 児玉郡上里村勅使河原 丹生神社 弘化3年(1846)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

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

長方形内に斜線を2本入れ，区分された領域に内接円を 4 個入れる。亨円，貞円の直径がそれぞれ 44 寸，33 寸のとき，利円の直径はいかほどか。

長方形の左下頂点を原点とし，長辺の長さを a とする（短辺の長さは元円の直径 2r1）
斜線と長辺の交点座標を (c, 0 ), (b, 0), (d, 2r1)
元円の半径と中心座標を r1, (a - r1, r1)
亨円の半径と中心座標を r2, (x2, 2r1 - r2)
利円の半径と中心座標を r3, (r3, r3)
貞円の半径と中心座標を r4, (x4, r4)
とおき，以下の連立方程式を解く。
eq1〜eq7 は円の中心から斜線への距離に関する方程式，eq8 は元円と亨円の半径の相似比についての方程式である。

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

using SymPy

@syms a, b, c, d, r1, r2, x2, r3, r4, x4

#(r2, r4) = (44, 33).//2

eq1 = dist(c, 0, d, 2r1, a - r1, r1) - r1^2
eq2 = dist(c, 0, d, 2r1, x2, 2r1 - r2) - r2^2
eq3 = dist(c, 0, d, 2r1, r3, r3) - r3^2
eq4 = dist(c, 0, d, 2r1, x4, r4) - r4^2
eq5 = dist(b, 0, 0, 2r1, x2, 2r1 - r2) - r2^2
eq6 = dist(b, 0, 0, 2r1, r3, r3) - r3^2
eq7 = dist(b, 0, 0, 2r1, x4, r4) - r4^2
eq8 = r2/x2 - r1/(a - r1);

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, b, c, d, r1, x2, r3, x4) = u
   return [
       -r1^2 + (-2*r1*(2*r1^2 + (-c + d)*(a - c - r1))/(4*r1^2 + (-c + d)^2) + r1)^2 + (a - c - r1 - (-c + d)*(2*r1^2 + (-c + d)*(a - c - r1))/(4*r1^2 + (-c + d)^2))^2,  # eq1
       -r2^2 + (-c + x2 - (-c + d)*(2*r1*(2*r1 - r2) + (-c + d)*(-c + x2))/(4*r1^2 + (-c + d)^2))^2 + (2*r1 - 2*r1*(2*r1*(2*r1 - r2) + (-c + d)*(-c + x2))/(4*r1^2 + (-c + d)^2) - r2)^2,  # eq2
       -r3^2 + (-2*r1*(2*r1*r3 + (-c + d)*(-c + r3))/(4*r1^2 + (-c + d)^2) + r3)^2 + (-c + r3 - (-c + d)*(2*r1*r3 + (-c + d)*(-c + r3))/(4*r1^2 + (-c + d)^2))^2,  # eq3
       -r4^2 + (-2*r1*(2*r1*r4 + (-c + d)*(-c + x4))/(4*r1^2 + (-c + d)^2) + r4)^2 + (-c + x4 - (-c + d)*(2*r1*r4 + (-c + d)*(-c + x4))/(4*r1^2 + (-c + d)^2))^2,  # eq4
       -r2^2 + (-b + b*(-b*(-b + x2) + 2*r1*(2*r1 - r2))/(b^2 + 4*r1^2) + x2)^2 + (2*r1 - 2*r1*(-b*(-b + x2) + 2*r1*(2*r1 - r2))/(b^2 + 4*r1^2) - r2)^2,  # eq5
       -r3^2 + (-2*r1*(-b*(-b + r3) + 2*r1*r3)/(b^2 + 4*r1^2) + r3)^2 + (-b + b*(-b*(-b + r3) + 2*r1*r3)/(b^2 + 4*r1^2) + r3)^2,  # eq6
       -r4^2 + (-2*r1*(-b*(-b + x4) + 2*r1*r4)/(b^2 + 4*r1^2) + r4)^2 + (-b + b*(-b*(-b + x4) + 2*r1*r4)/(b^2 + 4*r1^2) + x4)^2,  # eq7
       -r1/(a - r1) + r2/x2,  # eq8
   ]
end;

(r2, r4) = (44, 33)./2
iniv = BigFloat[205, 155, 40, 150, 42, 90, 30, 90]
res = nls(H, ini=iniv)

   (BigFloat[209.9999999999999999999999999999999999999999999999999999999999999998827654490999, 157.4999999999999999999999999999999999999999999999999999999999999999096530095862, 41.99999999999999999999999999999999999999999999999999999999999999999764231208424, 153.9999999999999999999999999999999999999999999999999999999999999998816882373083, 42.0000000000000000000000000000000000000000000000000000000000000000018357614881, 87.99999999999999999999999999999999999999999999999999999999999999994964199517249, 31.49999999999999999999999999999999999999999999999999999999999999999749742450348, 91.49999999999999999999999999999999999999999999999999999999999999994744891428587], true)

元円，利円の半径はそれぞれ 42, 31.5 である。
したがって，亨円，貞円の直径がそれぞれ 44 寸，33 寸のとき，利円の直径は 63 寸である。

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

   a = 210;  b = 157.5;  c = 42;  d = 154;  r1 = 42;  r2 = 22;  x2 = 88;  r3 = 31.5;  r4 = 16.5;  x4 = 91.5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, c, d, r1, x2, r3, x4) = [189, 113, 52, 100, 54, 62, 32, 76]
   (a, b, c, d, r1, x2, r3, x4) = res[1]
   (r2, r4) = (44, 33)./2
   println("利円の直径 = $(round(Int, Float64(2r3)))")
   @printf("a = %g;  b = %g;  c = %g;  d = %g;  r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  r4 = %g;  x4 = %g\n", a, b, c, d, r1, r2, x2, r3, r4, x4)
   plot([0, a, a, 0, 0], [0, 0, 2r1, 2r1, 0], color=:black, lw=0.5)
   circle(a - r1, r1, r1)
   circle(x2, 2r1 - r2, r2, :blue)
   circle(r3, r3, r3, :green)
   circle(x4, r4, r4, :magenta)
   segment(b, 0, 0, 2r1)
   segment(c, 0, d, 2r1)
   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 - r1, r1, "元円:r1,(a-r1,r1)", :red, :center, delta=-2delta)
       point(x2, 2r1 - r2, "亨円:r2\n(x2,2r1-r2", :blue, :center, :bottom, delta=2delta)
       point(r3, r3, "利円:r3,(r3,r3)", :green, :center, delta=-2delta)
       point(x4, r4, "貞円:r4,(x4,r4)", :black, :center, delta=-2delta)
       point(a, 0, "a", :black, :center, :top, delta=-delta)
       point(b, 0, "b", :black, :center, :top, delta=-delta)
       point(c, 0, "c", :black, :center, :top, delta=-delta)
       point(d, 2r1, "(d,2r1)", :black, :center, :bottom, delta=delta)
       point(0, 2r1, " 2r1", :black, :left, :bottom, delta=delta)
       plot!(ylims=(-10, 90))
   end
end;