算額（その107）

算額（その107）

埼玉県行田市 河原神社 万延2年(1861)3月
http://www.wasan.jp/saitama/kawara.html

https://yamabukiwasan.sakura.ne.jp/ymbk27.pdf

正方形の一辺の長さが 31 寸，甲円の直径が 18 寸のとき，乙円の直径を求めなさい。

乙円の半径を r1 とし，図のように記号を定め甲円，乙円が外接するという方程式を解く。

using SymPy

@syms r1::positive, r2::positive;

eq = (31//2 - r1)^2 + r1^2 - (r1 + 9)^2
solve(eq)[1] |> println

   7/2

乙円の半径は 7/2 である。もとの単位での直径は 7 寸である。

using Plots

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, color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   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 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=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 7/2
   plot()
   circle(0, 0, 9)
   circle(31/2 - r1, r1, r1, :blue)
   circle(31/2 - r1, -r1, r1, :blue)
   circle(-(31/2 - r1), r1, r1, :blue)
   circle(-(31/2 - r1), -r1, r1, :blue)

   circle(r1, 31/2 - r1, r1, :blue)
   circle(-r1, 31/2 - r1, r1, :blue)
   circle(r1, -(31/2 - r1), r1, :blue)
   circle(-r1, -(31/2 - r1), r1, :blue)

   a = 31/2
   plot!([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:black, lw=0.5)
    if more
       point(9, 0, "9 ", :red, :bottom, :right)
       point(31/2 - r1, r1, "(31/2-r1,r1)", :blue, :center)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
    end
end;

 

算額（その106）

算額（その106）

埼玉県行田市 河原神社 万延2年(1861)3月
http://www.wasan.jp/saitama/kawara.html

https://yamabukiwasan.sakura.ne.jp/ymbk27.pdf

菱形の中に大小三つの円があります。菱形の長軸が20寸、短軸が15寸の時、小円の直径は何寸ですか？

大円，小円の半径を r1, r2 とする。下図のように，三角形 OAB, AaO, Abc が相似であることを使えば，SymPy など使わなくても解ける。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
l = sympy.Line(p1, p2)
l.distance(sympy.Point(x0, y0))
end;

@syms r1::positive, r2::positive;

原点から直線 AB までの距離が r1 である。

eq1 = distance(0, 15/2, 20/2, 0, 0, 0) - r1;
solve(eq1)[1] |> println

   6

同じように「点 c から直線 AB までの距離が r2 である」を解こうとすると，'NotImplementedError' になる。

eq2 = distance(0, 15/2, 20/2, 0, r1 + r2, 0) - r2;
eq2 |> println
1.　solve(eq2(r1 => 6))

   -r2 + sqrt((9*r1/25 + 9*r2/25 - 18/5)^2 + (12*r1/25 + 12*r2/25 - 24/5)^2)

原因は sqrt() の中が複雑なためのようだ。

二乗距離にすれば解ける。

eq3 = distance(0, 15/2, 20/2, 0, r1 + r2, 0)^2 - r2^2
eq3 |> println

   -r2^2 + (9*r1/25 + 9*r2/25 - 18/5)^2 + (12*r1/25 + 12*r2/25 - 24/5)^2

小円の半径は 3/2 である。もとの単位での直径は 3 寸である。

solve(eq3(r1 => 6))[1] |> println

   3/2

@syms r1::positive, r2::positive;
solve([eq1, eq3], (r1, r2))

   1-element Vector{Tuple{Sym, Sym}}:
    (6, 3/2)

using Plots

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, color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   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 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=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (6, 3/2)
   plot()
   circle(0, 0, r1)
   circle(15/2, 0, 3/2, :blue)
   circle(-15/2, 0, 3/2, :blue)
   plot!([20/2, 0, -20/2, 0, 20/2], [0, 15/2, 0, -15/2, 0], color=:black, lw=0.5)
    if more
       cosine = 4/5
       sine = 3/5
       (ax, ay) = (20/2*sine^2, 20/2*cosine*sine)
       (bx, by) = ((20/2 - r1 - r2)sine^2 + r1 + r2, (20/2 - r1 - r2)cosine*sine)
       point(ax, ay, " a", :red, :left, :bottom) 
       point(bx, by, " b", :blue, :left, :bottom) 
       point(r1 + r2, 0, "c", :blue)
       point(0, 0, " O", :black)
       point(20/2, 0, " A", :black)
       point(0, 15/2, " B", :black, :left, :bottom)
       segment(ax, ay, 0, 0)
       segment(bx, by, r1 + r2, 0)
       point(ax/2, ay/2, "r1", :red, mark=false)
       point((bx + r1 + r2)/2, by/2, "r2", :blue, mark=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
    end
end;