算額（その174）

算額（その174）

岐阜県大垣市外野釜笛 釜笛八幡神社 慶応元年(1865)
http://ryugen3.sakura.ne.jp/toukou3/wasankibousya.PDF

第7問: 5 個の正方形，赤，黃，青，黒，白の一辺の長さの差は等しい。赤と黃の面積の和と青と黒と白の面積の和を知って，白の面積を求めよ。
白の一辺の長さを w，それぞれの一辺の差を d とすると，
赤の正方形の面積 = (w + 4d)^2
黃の正方形の面積 = (w + 3d)^2
青の正方形の面積 = (w + 2d)^2
黒の正方形の面積 = (w + d)^2
白の正方形の面積 = w^2
である。
赤と黃の面積の和を RY
青と黒と白の面積の和を BKW
とし，方程式を立てる。

以下は w = d = 1 の場合

using SymPy

@syms w::positive, d::positive, RY::positive, BKW::positive;

eq1 = (w + 4d)^2 + (w + 3d)^2 - RY |> expand
eq2 = (w +2d)^2 + (w + d)^2 + w^2 - BKW|> expand;

連立方程式を BKW, RY について解く。

res = solve([eq1, eq2], (w, d))

   2-element Vector{Tuple{Sym, Sym}}:
    ((-42*BKW + 18*RY - 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2))*sqrt(16*BKW/365 + 21*RY/365 - 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2)/365)/(6*(2*BKW - 3*RY)), sqrt(16*BKW/365 + 21*RY/365 - 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2)/365))
    ((-42*BKW + 18*RY + 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2))*sqrt(16*BKW/365 + 21*RY/365 + 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2)/365)/(6*(2*BKW - 3*RY)), sqrt(16*BKW/365 + 21*RY/365 + 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2)/365))

2 組の解が得られるが，BKW = 20, RY = 50 のときを計算してみると

res[1][1](BKW => 20, RY => 50).evalf(), res[1][2](BKW => 20, RY => 50).evalf()

   (1.43525126802703, 1.01117744530322)

res[2][1](BKW => 20, RY => 50).evalf(), res[2][2](BKW => 20, RY => 50).evalf()

   (-4.07737480549386, 2.54644251637035)

2 番目の解では w が負の値になるので不適切解である。

f(RY, BKW) = ((-42*BKW + 18*RY - 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2))*sqrt(16*BKW/365 + 21*RY/365 - 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2)/365)/(6*(2*BKW - 3*RY)), sqrt(16*BKW/365 + 21*RY/365 - 6*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2)/365));

(w, d) = f(41, 14)
println("w = $w,  d = $d")
println("RY  = $((w + 4d)^2 + (w + 3d)^2)")
println("BKW = $((w +2d)^2 + (w + d)^2 + w^2)")

   w = 1.0000000000000002,  d = 1.0000000000000002
   RY  = 41.000000000000014
   BKW = 14.000000000000007

白の正方形の面積は res[1][1]^2 で，簡約化して以下を得る。

res[1][1]^2 |> simplify |> factor |> println

   (213*BKW - 17*RY - 16*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2))/365

RY = 41
BKW = 14
g(RY, BKW) = (213*BKW - 17*RY - 16*sqrt(-BKW^2 + 43*BKW*RY - 6*RY^2))/365  # 白の面積
g(41, 14) |> println
g(50, 20) |> println

   1.0
   2.059946202373195

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
   end
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function rect(x1, y1, x2, y2, color=:pink)
   plot!([x1, x2, x2, x1, x1], [y1, y1, y2,  y2, y1], color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   colors = [:indianred1, :yellow, :steelblue1, :gray30, :white]
   plot()
   for i = 5:-1:1
       width = w + (i-1)d
       rect(0, 0, width, width, colors[6-i])
   end
   if more == true
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;