算額（その1600）
宮城県角田市小田字斗蔵 斗蔵寺 明治42年(1909)
https://tajin.shiriagari.com/framepage5000.htm
（リンク先の左のインデックスから選択）
キーワード：二等辺三角形，正方形2個
#Julia, #SymPy, #算額, #和算, #数学
二等辺三角形内に大小の正方形を容れる。二等辺三角形の高さが 25 寸のとき，小正方形が最大になるときの二等辺三角形の底辺の長さと小正方形の一辺の長さはいかほどか。
二等辺三角形の高さを h，底辺の長さを 2a
大正方形の一辺の長さを 2b
小正方形の一辺の長さを 2c
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a, b, c, h
eq1 = 2b/(a - b) - h/a
eq2 = 2c/(b - c) - h/a
eq3 = (h - 2b - 2c)/c - h/a
res = solve([eq1, eq2], (c, b));
# c
res[c] |> println
a*h^2/(4*a^2 + 4*a*h + h^2)
h = 25 のとき，c は a の関数で下図のような関係にある。
pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(res[c](h => 25), xlims=(0, 30), xlabel="a", ylabel="c")
a が 10 近辺で，c が最大になるようだ。
最大値を取るときの a の値は，関数を微分し，0 になるときの値を求めればよい。
diff(res[c], a) |> simplify |> println
h^2*(-2*a + h)/(8*a^3 + 12*a^2*h + 6*a*h^2 + h^3)
ans_a = solve(diff(res[c], a), a)[1]
ans_a |> println
h/2
res[c](a => h/2)(h => 25) |> println
res[c](a => h/2)(h => 25).evalf() |> println
25/8
3.12500000000000
a = h/2 = 12.5 のとき，すなわち高さと底辺の長さが等しいときに小正方形の一辺の長さが最大 2*3.125 = 6.25 寸となる。
答では，小正方形の一辺の長さは 5 寸。
術では，高さを 5 で割る。
となっているが，それは上に示した最大値とは違う。小正方形の一辺の長さが 5 になるときの a を求めると 32.7254248593737，4.77457514062631 である。この間に c が最大になるときの a がある。
eq = a*h^2/(4*a^2 + 4*a*h + h^2) - c
ans_a = solve(eq, a)
ans_a |> println
Sym{PyCall.PyObject}[(-h*(4*c - h) + sqrt(-h^3*(8*c - h)))/(8*c), (-4*c*h + h^2 - sqrt(h^3*(-8*c + h)))/(8*c)]
ans_a[1](h => 25, c => 5//2) |> println
ans_a[1](h => 25, c => 5//2).evalf() |> println
25*sqrt(5)/4 + 75/4
32.7254248593737
ans_a[2](h => 25, c => 5//2) |> println
ans_a[2](h => 25, c => 5//2).evalf() |> println
75/4 - 25*sqrt(5)/4
4.77457514062631
function draw(h, more=false; a = 0)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
a == 0 && (a = h/2)
b = a*h/(2a + h)
c = a*h^2/(4*a^2 + 4*a*h + h^2)
str = @sprintf("h = %g\na = %g\nb = %g\nc = %g", h, a, b, c)
p = plot([a, 0, -a, a], [0, h, 0, 0], color=:green, lw=0.5)
rect(-b, 0, b, 2b, :blue)
rect(-c, 2b, c, 2b + 2c, :blue)
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
point(-b/2, 1.8b, str, :red, :left, mark=false)
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)
end
return p
end;
draw(25, true; a = 4.77457514062631)
