算額（その1602）
群馬県みどり市大間々町 鎮守社 文政4年(1821)
深川英俊，トニー・ロスマン：聖なる数学：算額，森北出版株式会社，2010年４月22日.
キーワード：菱形2個，面積
#Julia, #SymPy, #算額, #和算, #数学
菱形の中に正方形を容れる。菱形の面積から正方形の面積を差し引いた面積（黒積と呼ぼう）を考える。菱形の一辺の長さが与えられたとき，黒積が最大となるときの正方形の一辺の長さはいかほどか。
菱形の一辺の長さを a
正方形の対角線の長さを 2b
正方形の一辺の長さを x
とおき，以下の方程式を解く。
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, x::positive
c = sqrt(a^2 - b^2)
x = √Sym(2)*b
黒積 = 2b*c - x^2;
黒積は a, b の関数として表される。
黒積 |> println
-2*b^2 + 2*b*sqrt(a^2 - b^2)
たとえば a = 10 のとき，下図のように b が 4 前後で黒積が最大になる。
pyplot(size=(300, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(黒積(a => 10), xlims=(0, 7), xlabel="b", ylabel="黒積")
黒積が最大値を取るときの b の値は，黒積を表す関数を微分し，微分係数が 0 になるときの値を求めればよい。
diff(黒積, b) |> simplify |> println
2*(a^2 - 2*b^2 - 2*b*sqrt(a^2 - b^2))/sqrt(a^2 - b^2)
ans_b = solve(diff(黒積, b), b)[1] |> simplify
ans_b |> println
a*sqrt(2 - sqrt(2))/2
b = a*sqrt(2 - √2)/2 のとき，導関数は 0 になり，黒積は最大値 a^2*(√2 - 1) を取る。
黒積 = 黒積(b => ans_b) |> simplify
黒積 |> println
a^2*(-1 + sqrt(2))
たとえば，a = 10 のとき，b = 3.82683432365090 のときに黒積は最大値 41.4213562373095 を取る。
ans_b(a => 10).evalf() |> println
黒積(a => 10).evalf() |> println
3.82683432365090
41.4213562373095
正方形の一辺の長さは x = √2b = √2*ans_b = a*sqrt(4 - 2√2)/2 である。
x = √Sym(2)*ans_b
x |> simplify |> println
a*sqrt(4 - 2*sqrt(2))/2
function draw(a, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
b = a*sqrt(2 - √2)/2
c = sqrt(a^2 - b^2)
x = a*sqrt(4 - 2√2)/2
@printf("a = %g; b = %g; c = %g; x = %g\n", a, b, c, x)
plot([c, 0, -c, 0, c], [0, b, 0, -b, 0], seriestype=:shape, color=:gray80, fillcolor=:gray80, lw=0.5)
plot!([b, 0, -b, 0, b], [0, b, 0, -b, 0], seriestype=:shape, color=:white, fillcolor=:white, lw=0.5)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray40, lw=0.5)
vline!([0], color=:gray40, lw=0.5)
point(0, b, "b ", :blue, :right, :bottom, delta=delta)
point(c, 0, "c", :blue, :center, delta=-delta)
point(b, 0, "b", :blue, :center, :bottom, delta=delta)
dimension_line(2delta, b + 4delta, c + 2delta, 4delta, "a", :blue, :left, :bottom, delta=delta)
dimension_line(- 3delta, b-3delta, b - 3delta, -3delta, "x", :blue, :right, :vcenter, deltax = -6delta)
end
end;
draw(10, true)
