算額（その1601）
宮城県角田市小田字斗蔵 斗蔵寺 明治42年(1909)
https://tajin.shiriagari.com/framepage5000.htm
（リンク先の左のインデックスから選択）
キーワード：円5個，最大値
#Julia, #SymPy, #算額, #和算, #数学
交差する大円 2 個の隙間に，中円 1 個，小円 2 個を容れる。大円の直径が 10 寸のとき，小円の直径が最大になるのはどのようなときか，そしてそのときの小円の直径はいかほどか。
大円の半径と中心座標を r1, (r1 - r2, 0)
中円の半径と中心座標を r2, (0, 0)
小円の半径と中心座標を r3, (0, r2 + r3)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r1, r2, r3
eq = (r1 - r2)^2 + (r2 + r3)^2 - (r1 - r3)^2
res = solve(eq, r3)[1]
res |> println
r2*(r1 - r2)/(r1 + r2)
r3 は r2 の関数で，r1 が定数として固定されたとき，r3:r2 は下図のような関係にある。
たとえば，r1 が 10/2 のとき，r2 が 2 近辺で，r3 が最大になるようだ。
pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(res(r1 => 10/2), xlims=(0, 5), xlabel="r2", ylabel="r3")
最大値を取るときの r2 の値は，関数を微分し，0 になるときの値を求めればよい。
diff(res, r2) |> simplify |> println
(-r2*(r1 - r2) + (r1 - 2*r2)*(r1 + r2))/(r1 + r2)^2
ans_r2 = solve(diff(res, r2), r2)[1]
ans_r2 |> println
r1*(-1 + sqrt(2))
r1 = 10/2 のとき，r2 = 2.07106781186548 のとき，導関数は 0 になり，r3 は最大値を取る。
ans_r2(r1 => 10/2) |> println
ans_r2(r1 => 10/2).evalf() |> println
-5.0 + 5.0*sqrt(2)
2.07106781186548
r1 = 10/2，r2 = 2.07106781186548 のとき，r3 = 0.857864376269050 が最大値になる。
res(r2 => ans_r2)(r1 => 10/2) |> println
res(r2 => ans_r2)(r1 => 10/2).evalf() |> println
0.5*sqrt(2)*(-1 + sqrt(2))*(10.0 - 5.0*sqrt(2))
0.857864376269050
function draw(r1, more=false; r2 = 0)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 == 0 && (r2 = r1*(-1 + sqrt(2)))
r3 = r2*(r1 - r2)/(r1 + r2)
str = @sprintf("r1 = %g;; r2 = %.3g; r3 = %.3g\n", r1, r2, r3)
p = plot()
circle2(r1 - r2, 0, r1)
circle(0, 0, r2, :blue)
circle22(0, r2 + r3, 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(0, r1, str, :black, :center, :bottom, delta=delta/2, mark=false)
end
return p
end;
draw(10/2, true)
