算額（その47）

算額（その47）

第2版を公開しました。
https://blog.goo.ne.jp/r-de-r/e/e79c4d483dbfe55f6973570fa73a0c04

岩手県東磐井郡藤沢町 藤勢寺 元治 2 年
http://www.wasan.jp/iwate/fujise.html

1 辺が 1 寸の正方形に2種類の円が収まっている。それぞれの径を求めよ。

下図のように記号を定め解を求める。

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;

using SymPy
function intersection(x1, y1, x2, y2, x3, y3, x4, y4)
   p1, p2, p3, p4 = sympy.Point(x1, y1), sympy.Point(x2, y2), sympy.Point(x3, y3), sympy.Point(x4, y4)
   l1 = sympy.Line(p1, p2)
   l2 = sympy.Line(p3, p4)
   l1.intersection(l2)
end;

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

x1, x2 は半径 r1, r2 で表すことができる。

x1 = (1 + sqrt(Sym(2))) * r1
x2 = (1 - sqrt(Sym(2)) * r2)/2
(x1, x2) |> println

   (r1*(1 + sqrt(2)), -sqrt(2)*r2/2 + 1/2)

以下の連立方程式を解くが，solve() では解けないので，nlsolve() で解く。

eq1 = distance(r1, x1, 1, 0, x1, r1) - r1;
eq2 = distance(r1, x1, 1, 0, x2, x2) - r2;

# solve([eq1, eq2], (r1, r2))

eq1 |> println

   -r1 + sqrt((r1 - r1*(14*r1^2 + 10*sqrt(2)*r1^2 - 10*r1 - 7*sqrt(2)*r1 + 2*sqrt(2) + 3)/(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2)))^2 + (r1*(1 + sqrt(2)) - r1*(4*sqrt(2)*r1^2 + 6*r1^2 - 2*r1 - sqrt(2)*r1 + 2*sqrt(2) + 3)/(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2)))^2)

eq2 |> println

   -r2 + sqrt((-sqrt(2)*r2/2 + 1/2 - (-6*r1^2*r2 - 4*sqrt(2)*r1^2*r2 + 18*r1^2 + 13*sqrt(2)*r1^2 + 5*sqrt(2)*r1*r2 + 8*r1*r2 - 4*sqrt(2)*r1 - 5*r1 - 2*r2 - sqrt(2)*r2 + 1 + sqrt(2))/(2*(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2))))^2 + (-r1*(-10*sqrt(2)*r1*r2 - 14*r1*r2 + 4*r1 + 3*sqrt(2)*r1 + 4*r2 + 3*sqrt(2)*r2 + 2*sqrt(2) + 3)/(2*(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2))) - sqrt(2)*r2/2 + 1/2)^2)

using NLsolve

function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
   r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params...), [ini], ftol=1e-14)
   v = r.zero[1]
else
   r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
   v = r.zero
end
return v, r.f_converged
end;

function h(u)
   r1, r2 = u
   return [
       -r1 + sqrt((r1 - r1*(14*r1^2 + 10*sqrt(2)*r1^2 - 10*r1 - 7*sqrt(2)*r1 + 2*sqrt(2) + 3)/(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2)))^2 + (r1*(1 + sqrt(2)) - r1*(4*sqrt(2)*r1^2 + 6*r1^2 - 2*r1 - sqrt(2)*r1 + 2*sqrt(2) + 3)/(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2)))^2),
       -r2 + sqrt((-sqrt(2)*r2/2 + 1/2 - (-6*r1^2*r2 - 4*sqrt(2)*r1^2*r2 + 18*r1^2 + 13*sqrt(2)*r1^2 + 5*sqrt(2)*r1*r2 + 8*r1*r2 - 4*sqrt(2)*r1 - 5*r1 - 2*r2 - sqrt(2)*r2 + 1 + sqrt(2))/(2*(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2))))^2 + (-r1*(-10*sqrt(2)*r1*r2 - 14*r1*r2 + 4*r1 + 3*sqrt(2)*r1 + 4*r2 + 3*sqrt(2)*r2 + 2*sqrt(2) + 3)/(2*(8*r1^2 + 6*sqrt(2)*r1^2 - 2*sqrt(2)*r1 - 2*r1 + 1 + sqrt(2))) - sqrt(2)*r2/2 + 1/2)^2)
   ];
end;

iniv = [0.1, 0.2];

res = nls(h, ini=iniv)

   ([0.09990042250462969, 0.18946869098150596], true)

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
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 draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1, r2 = [0.09990042250462969, 0.18946869098150596]
   x1, x2 = r1*(1 + sqrt(2)), -sqrt(2)*r2/2 + 1/2
   res = intersection(0, 1, x1, r1, r1, x1, 1, 0)
   cros = res[1].coordinates[1]  # 3 直線の交点
   p = plot([0, 1, 1, 0, 0], [0, 0, 1, 1, 0], color=:black, lw=0.5)
   circle(r1, x1, r1)
   circle(x1, r1, r1)
   circle(1-r1, 1-x1, r1)
   circle(1-x1, 1-r1, r1)
   circle(x2, x2, r2, :red)
   circle(sqrt(2)*r2 + x2, sqrt(2)*r2 + x2, r2, :red)
   plot!([0, 1], [1, 0], color=:black, lw=0.5)
   plot!([cros, 0], [cros, 1], color=:black, lw=0.5)
   plot!([cros, 1], [cros, 0], color=:black, lw=0.5)
   plot!([cros, 0], [cros, 0], color=:black, lw=0.5)
   plot!([1-cros, 0], [1-cros, 1], color=:black, lw=0.5)
   plot!([1-cros, 1], [1-cros, 0], color=:black, lw=0.5)
   plot!([1-cros, 1], [1-cros, 1], color=:black, lw=0.5)
   if more
       point(x1, r1, "(x1,r1)", :green, :center)
       point(r1, x1, "(r1,x1)", :green, :center)
       point(x2, x2, "(x2,x2)", :green, :center)
   end
   display(p)
   savefig("fig1.png")
end
draw(true)

 

算額（その46）

算額（その46）

大阪府茨木市 井於神社 弘化3年
http://www.wasan.jp/osaka/iyo.html

円の下半分に垂直方向に等間隔の弦を引く。弦の長さを求めよ。

下図のように記号を定め解を求める。

using SymPy
@syms y(), x, y;
@syms R::positive;
@syms x1::positive, x2::positive, x3::positive;

using Plots

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
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 answer(more)
   R = 8  # 直径
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   p = plot()
   circle(0, 0, R/2)
   y = -sqrt((R//2)^2 - x^2)  # SymPy での関数定義（円の下半分）
   for i in 0:4
       y1 = -i*(4//5)  # 弦の y 座標
       eq1 = y - y1    # 方程式
       x1, x2 = solve(eq1, x)  # 円との交点
       println("$(i+1)番: $eq1;  x1 = $x1;  x2 = $x2;  長さ = $(2*x2.evalf())")
       plot!([x1, x2], [y1, y1])
       if more && i == 3
           point(x1, y1, "x1 ", :green, :right)
           point(x2, y1, " x2", :green, :left)
       end
   end
   savefig("fig1.png")
end
answer(true)

   1番: -sqrt(16 - x^2);  x1 = -4;  x2 = 4;  長さ = 8.00000000000000
   2番: 4/5 - sqrt(16 - x^2);  x1 = -8*sqrt(6)/5;  x2 = 8*sqrt(6)/5;  長さ = 7.83836717690617
   3番: 8/5 - sqrt(16 - x^2);  x1 = -4*sqrt(21)/5;  x2 = 4*sqrt(21)/5;  長さ = 7.33212111192934
   4番: 12/5 - sqrt(16 - x^2);  x1 = -16/5;  x2 = 16/5;  長さ = 6.40000000000000
   5番: 16/5 - sqrt(16 - x^2);  x1 = -12/5;  x2 = 12/5;  長さ = 4.80000000000000

算額では 2 番が 7寸8分6厘1毛，3 番が 7寸，，2毛となっているが，それぞれ 7.838, 7.332 である。