川田保知：『算法極数小補解義』 文化15年(1818)
山口正義：やまぶき，第26号
https://yamabukiwasan.sakura.ne.jp/ymbk26.pdf
キーワード：円5個
#Julia #SymPy #算額 #和算 #数学
2 個の甲円が交差してできる区画に，乙円と丙円を容れる。丙円の直径が最大になるときの乙円の直径を求めよ。
甲円の半径と中心座標を \(r_1,\ (r_1,\ 0),\ (r_1 + 2r_2,\ 0)\)
乙円の半径と中心座標を \(r_2,\ (r_2,\ 0)\)
丙円の半径と中心座標を \(r_3,\ (x_3,\ y_3)\)
とおき，\(r_2\) を変数のまま \(r_3,\ x_3,\ y_3\) を求める。
\(r_3\) は \(r_2\) の関数になるので，\(r_3\) を \(r_2\) で微分し，接線の傾きが 0 になるときの \(r_2\) を求める。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms r1::positive, r2::positive, r3::positive, x3::positive, y3::positive
r1 = 35//2
eq1 = (r1 - x3)^2 + y3^2 - (r1 - r3)^2
eq2 = (x3 - r2)^2 + y3^2 - (r2 + r3)^2
eq3 = (r1 + 2r2 - x3)^2 + y3^2 - (r1 + r3)^2
res1 = solve([eq1, eq2, eq3], (r3, x3, y3));
丙円の半径は，乙円の半径の関数である。
r3 = res1[1][1] |> simplify
r3 |> println
r2*(1225 - 4*r2^2)/(4*r2^2 + 1225)
res1[1][2] |> simplify |> println
res1[1][3] |> simplify |> println
r2*(2*r2 + 35)^2/(4*r2^2 + 1225)
70*r2*sqrt(1225 - 4*r2^2)/(4*r2^2 + 1225)
乙円の半径 \(r_2\) が 7.5〜10.0 の範囲内で丙円の半径 \(r_3\) が最大になる。
using Plots
pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(r3, xlims=(0, 17.5), xlabel="乙円の半径 r2", ylabel="丙円の半径 r3")
\(r_3\) を \(r_2\) で微分して，接線の傾きが 0 になるときの \(r_2\) を求めると 8.50269475574130 である。つまり，乙円の直径が 17.0053895114826 のとき丙円の直径が最も大きくなる。
f = diff(r3, r2)
res2 = solve(f, r2)[1]
res2 |> println
res2.evalf() |> println
35*sqrt(-1/2 + sqrt(5)/4)
8.50269475574130
もっとも大きいときの丙円の半径は 5.25495435501361（直径は 10.5099087100272 である）。
2r3(r2 => res2[1]).evalf() |> println
10.5099087100272
その他のパラメータは以下のとおりである。
\(r_1 = 17.5;\ r_2 = 8.50269;\ r_3 = 5.25495\)
\(x_3 = 15.1871;\ y_3 = 12.0246\)
描画関数プログラムのソースを見る
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 35//2
r2 = 35sqrt(√5 - 2)/2
t = r2/(4r2^2 + 1225)
(r3, x3, y3) = t .* (
1225 - 4r2^2,
(2*r2 + 35)^2,
70sqrt(1225 - 4r2^2))
@printf("乙円の直径が %g のとき，丙円の直径は最大値 %g になる。\n", 2r2, 2r3)
@printf("r1 = %g; r2 = %g; r3 = %g; x3 = %g; y3 = %g\n", r1, r2, r3, x3, y3)
plot()
circle(r1, 0, r1)
circle(r1 + 2r2, 0, r1)
circle(r2, 0, r2, :blue)
circle22(x3, y3, r3, :green)
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(r1, 0, " 甲円:r1,(r1,0)", :red, :left, delta=-delta)
point(r1 + 2r2, 0, " 甲円:r1,(r1+2r2,0)", :red, :left, delta=-delta)
point(r2, 0, "乙円:r2,(r2,0)", :blue, :center, delta=-delta)
point(x3, y3, "丙円:r3\n(x3,y3)", :green, :center, :bottom, delta=delta)
end
end;
以下のアイコンをクリックして応援してください
