算額（その602）

算額（その602）

長野県東御市祢津 長明寺・大日如来堂 嘉永6年(1853)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html

外円内に長方形を 5 個入れる。それぞれの長方形の長辺の 2 頂点は外円の円周上にあり，残りの 2 頂点は隣の長方形と共有している。
5 個の面積の和が最大になるのは，長辺がいかほどのときか。

外円の半径を R, 内側の 5 角形が内接する円の半径を r として，長方形の頂点座標 (x, y), (x, sqrt(R^2 - x^2)) を求める。
長方形の面積の和は 5(2x * (sqrt(R^2 - x^2) - y)) である。

include("julia-source.txt");

using SymPy
@syms R, r
R = 55//2
(x, y) = r .* (cosd(Sym(54)), sind(Sym(54)))
eq = 5(2x * (sqrt(R^2 - x^2) - y))
eq |> println

   10*r*sqrt(5/8 - sqrt(5)/8)*(-r*(1/4 + sqrt(5)/4) + sqrt(-r^2*(5/8 - sqrt(5)/8) + 3025/4))

この式の値は r により変化するが，r = 14 近辺で最大値を取ることがわかる。

using Plots
pyplot(size=(400, 250), grid=false, aspectratio=:none, label="")
plot(eq, xlims=(0,55/2))
vline!([sqrt(-1/(-28 + 12*sqrt(5)))*sqrt(15125 - 6655*sqrt(5))])
hline!([1228.6026013806772])

eq が最大値となるのは，eq の導関数が 0 になるときである。

g = diff(eq);
plot(g, xlims=(0,55/2))
vline!([sqrt(-1/(-28 + 12*sqrt(5)))*sqrt(15125 - 6655*sqrt(5))])
hline!([0])

導関数が 0 になるのは r = 14.4576055832762 のときである。

res = solve(g)
res |> println

   Sym[-sqrt(-1/(-28 + 12*sqrt(5)))*sqrt(9075 - 3025*sqrt(5)), sqrt(-1/(-28 + 12*sqrt(5)))*sqrt(15125 - 6655*sqrt(5))]

2 番目のものが適解である。

res[2] |> println

   sqrt(-1/(-28 + 12*sqrt(5)))*sqrt(15125 - 6655*sqrt(5))

res[2].evalf() |> println

   14.4576055832762

最大値は 1228.60260138068 である。

eq(r => 14.4576055832762).evalf() |> println

   1228.60260138068

r = 14.4576055832762 のとき，長方形の長辺が 16.995934690622075，短辺が 14.457605583276203 である。

r = 14.4576055832762
2r*cosd(54) |> println
sqrt(R^2 - (r*cosd(54))^2) - r*sind(54) |> println

   16.99593469062214
   14.457605583276148

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 55/2
   r = sqrt(-1/(-28 + 12*sqrt(5)))*sqrt(15125 - 6655*sqrt(5))
   plot()
   circle(0, 0, R, :blue)
   circle(0, 0, r, :blue)
   (x0, y0) = R .* (cosd(54), sind(54))
   segment(0, 0, x0, y0)
   segment(0, 0, -x0, y0)
   (x, y) = r .* (cosd(54), sind(54))
   rect(-x, y, x, sqrt(R^2 - x^2), :red)
   long = 2x
   short = sqrt(R^2 - x^2) - y
   println("長方形: 長辺 = $long; 短辺 = $short;  5個の面積 = $(5long*short)")
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x, y, " (x,y)", :red, :left, :vcenter)
       point(x, sqrt(R^2 - x^2), " (x,sqrt(R^2-x^2))", :red, :left, :vcenter)
       point(R, 0, "R ", :blue, :right, :bottom, delta=delta/2)
       point(r, 0, "r ", :blue, :right, :bottom, delta=delta/2)
       circle(0, 0, 2, :green, beginangle=0, endangle=54)
       circle(0, 0, 2.2, :green, beginangle=0, endangle=54)
       point(2.4, 2, "θ=54°", mark=false)
   end
end;

算額（その601）

算額（その601）

長野県東御市祢津 長明寺・大日如来堂 嘉永6年(1853)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html

外円内に弦，矢，大円，小円を入れる。大円，小円の直径がそれぞれ 16 寸，12 寸，矢の長さが 15 寸のとき，外円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (x + r1, y + r1)
小円の半径と中心座標を r2, (x - r2, y + r2)
矢の長さを「矢」，弦と y 軸の交点座標を (x, y)
として以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R, 矢, r1, r2, x, y
eq1 = (x + r1)^2 + (y + r1)^2 - (R - r1)^2
eq2 = (x - r2)^2 + (y + r2)^2 - (R - r2)^2
eq3 = x^2 + (y + 矢)^2 - R^2;
res = solve([eq1, eq2, eq3], (R, x, y))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (-(r1^2 - 2*r1*矢 + r2^2 - 2*r2*矢 + 2*矢^2)/(4*(r1 + r2 - 矢)), -(r1 - r2)*(r1 + r2 - 2*矢)/(4*(r1 + r2 - 矢)), (-r1*r2 - r1*矢 - r2*矢 + 矢^2)/(2*r1 + 2*r2 - 2*矢))
    (-(2*r1^2*r2^2 - 2*r1^2*r2*矢 + r1^2*矢^2 - 2*r1*r2^2*矢 + r2^2*矢^2)/(4*r1*r2*矢), -(r1 - r2)*(2*r1*r2 - r1*矢 - r2*矢)/(4*r1*r2), (r1*r2 - r1*矢 - r2*矢 - 矢^2)/(2*矢))

2 組の解が得られるが，最初のものが適解である。

(r1, r2) = (16, 12) ./ 2
矢 = 15
(-(r1^2 - 2*r1*矢 + r2^2 - 2*r2*矢 + 2*矢^2)/(4*(r1 + r2 - 矢)), -(r1 - r2)*(r1 + r2 - 2*矢)/(4*(r1 + r2 - 矢)), (-r1*r2 - r1*矢 - r2*矢 + 矢^2)/(2*r1 + 2*r2 - 2*矢)),
(-(2*r1^2*r2^2 - 2*r1^2*r2*矢 + r1^2*矢^2 - 2*r1*r2^2*矢 + r2^2*矢^2)/(4*r1*r2*矢), -(r1 - r2)*(2*r1*r2 - r1*矢 - r2*矢)/(4*r1*r2), (r1*r2 - r1*矢 - r2*矢 - 矢^2)/(2*矢))

   ((32.5, -8.0, 16.5), (-2.4125, 1.1875, -12.9))

外円の直径は 65 寸である。

その他のパラメータは以下の通りである。
r1 = 8;  r2 = 6;  矢 = 15;  R = 32.5;  x = -8;  y = 16.5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (16, 12) .// 2
   矢 = 15
   (R, x, y) = (-(r1^2 - 2*r1*矢 + r2^2 - 2*r2*矢 + 2*矢^2)/(4*(r1 + r2 - 矢)), -(r1 - r2)*(r1 + r2 - 2*矢)/(4*(r1 + r2 - 矢)), (-r1*r2 - r1*矢 - r2*矢 + 矢^2)/(2*r1 + 2*r2 - 2*矢))
   @printf("外円の直径 = %g;  r1 = %g;  r2 = %g;  矢 = %g;  R = %g;  x = %g;  y = %g\n", 2R, r1, r2, 矢, R, x, y)
   plot()
   circle(0, 0, R, :blue)
   circle(x + r1, y + r1, r1, :green)
   circle(x - r2, y + r2, r2, :red)
   segment(-sqrt(R^2 - y^2), y, sqrt(R^2 - y^2), y)
   segment(x, y, x, y + 矢, :brown, lw=2)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x + r1, y + r1, " 大円:r1\n (x+r1,y+r1)", :black, :left, :vcenter)
       point(x - r2, y + r2, "小円:r2\n(x-r2,y+r2)", :black, :center, delta=-delta/2)
       point(x, y + 矢/2, " 矢", :brown, mark=false)
       point(0, R, " R", :blue, :left, :bottom, delta=delta/2)
       point(0, y, " y", :black, :left, :top, delta=-delta/2)
       point(x, y, "(x,y)", :black, :center, :top, delta=-delta/2)
   end
end;