算額について中間報告

まあ，SymPy 最強ということかな

解けない算額は少ない（解けないのは，小生の力不足）

算額（その93）

算額（その93）

千葉県君津市 神野寺 明治20(1887)年4月
http://www.wasan.jp/chiba/jinyaji2.html

甲円の下に乙円，丙円，丁円がある。各円の径を求めよ。

甲円の半径を 1 として，以下のように記号を定め，方程式を解く。
なお，solve() では解が求まらないので，nlsolve() を使用し，数値的に解を求める。

using SymPy
@syms r1::positive, r2::positive, r3::positive, x1::positive, x2::positive, x3::positive, x::positive, y::positive;
r = 1
x1 = r1 
eq1 = x*x3 + y*(r3 - y) |> expand
eq2 = r3*x - y*(x - x3) |> expand
eq3 = x1^2 + (y - r1)^2 - (r + r1)^2 |> expand
eq4 = x2^2 + (y - r2)^2 - (r + r2)^2 |> expand
eq5 = x3^2 + (y - r3)^2 - (r + r3)^2 |> expand
eq6 = (x2 - x1)^2 + (r2 - r1)^2 - (r2 + r1)^2 |> expand
eq7 = (x3 - x2)^2 + (r3 - r2)^2 - (r3 + r2)^2 |> expand;

# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r1, r2, r3, x2, x3, x, y))

eq1 |> println
eq2 |> println
eq3 |> println
eq4 |> println
eq5 |> println
eq6 |> println
eq7 |> println

   r3*y + x*x3 - y^2
   r3*x - x*y + x3*y
   r1^2 - 2*r1*y - 2*r1 + y^2 - 1
   -2*r2*y - 2*r2 + x2^2 + y^2 - 1
   -2*r3*y - 2*r3 + x3^2 + y^2 - 1
   r1^2 - 4*r1*r2 - 2*r1*x2 + x2^2
   -4*r2*r3 + x2^2 - 2*x2*x3 + x3^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, r3, x2, x3, x, y) = u
return [
       r3*y + x*x3 - y^2,
       r3*x - x*y + x3*y,
       r1^2 - 2*r1*y - 2*r1 + y^2 - 1,
       -2*r2*y - 2*r2 + x2^2 + y^2 - 1,
       -2*r3*y - 2*r3 + x3^2 + y^2 - 1,
       r1^2 - 4*r1*r2 - 2*r1*x2 + x2^2,
       -4*r2*r3 + x2^2 - 2*x2*x3 + x3^2
      ];
end;

iniv = [0.1, 0.2, 0.3, 0.4, 0.7, 1.0, 1.1];
res = nls(H, ini=iniv)

   ([0.13839876084731936, 0.17954306798152128, 0.3287706665765487, 0.45366714256276236, 0.9395827489780465, 1.2683534155545952, 1.268353415554595], 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; fontsize=10, mark=true)
  mark && scatter!([x], [y], color=color, markerstrokewidth=0)
  annotate!(x, y, text(string, fontsize, vertical, position, 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=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, x2, x3, x, y) = [0.13839876084731936, 0.17954306798152128, 0.3287706665765487,
       0.45366714256276236, 0.9395827489780465, 1.2683534155545952, 1.268353415554595]
   r = 1
   plot()
   circle(0, y, 1)
   circle(r1, r1, r1, :green)
   circle(-r1, r1, r1, :green)
   circle(x2, r2, r2, :blue)
   circle(-x2, r2, r2, :blue)
   circle(x3, r3, r3, :magenta)
   circle(-x3, r3, r3, :magenta)
   hline!([0], color=:black, lw=0.5)
   if more
       point(0, y, "y ", :red, :right)
       point(0, y-r, "y-r ", :red, :right)
       point(r1, r1, "(r1,r1)", :green, :top)
       point(x2, r2, "(x2,r2)", :green, :top)
       point(x3, r3, "(x3,r3)", :magenta, :top)
       point(x, 0, "   x", :magenta, :bottom)
       segment(0, y, x, 0, :black)
       vline!([0], color=:black, lw=0.5)
   end
end;

算額（その92）

算額（その92）

千葉県君津市 鹿野山神野寺 明治20(1887)年4月
http://www.wasan.jp/chiba/jinyaji4.html

山口正義(2018): やまぶき4（第58号）
https://yamabukiwasan.sakura.ne.jp/ymbk58.pdf

長方形内に大円，中円，小円入っている。長方形の長辺の長さが25寸のとき，小円の径を求めよ。

長方形の長辺の長さを 2a，短辺の長さを b とおく。
大円の半径と中心座標を r1, (r1, b - r1)
中円の半径と中心座標を r2, (0, r2), (2r2, r2)
小円の半径と中心座標を r3, (0, 2r2 + r3)
とおき，連立方程式を解く。

include("julia-source.txt")

using SymPy
@syms r1::positive, r2::positive, r3::positive, a::positive, b::positive;
(r1, r2) = a .* (1//2, 1//3)
eq1 =(r1 - r2)^2 + (b - r1 - r2)^2 - (r1 + r2)^2
eq2 = r1^2 + (b - r1 - (2r2 + r3))^2 - (r1 + r3)^2
res = solve([eq1, eq2], (b, r3))

   1-element Vector{Tuple{Sym, Sym}}:
    (23*a*(-31/69 + 16*sqrt(6)/69)/16 + 71*a/48, a*(-31/69 + 16*sqrt(6)/69))

   r3 = 1.48403;  b = 20.6229
   小円の直径 = 2.96806

小円の半径は a*(16√6 - 31)/69 である。

2a = 25 寸のとき，小円の半径は 1.4840282399512403寸（直径は 2.968056479902481寸）である。

(25/2)*(16√6 - 31)/69

   1.4840282399512403

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 25//2
   (r1, r2) = a .* (1//2, 1//3)
   (b, r3) = (23*a*(-31/69 + 16*sqrt(6)/69)/16 + 71*a/48, a*(-31/69 + 16*sqrt(6)/69))
   @printf("r3 = %g;  b = %g\n", r3, b)
   @printf("小円の直径 = %g\n", 2r3)
   plot([a, a, -a, -a, a], [0, b, b, 0, 0], color=:black, lw=0.5)
   circle(r1, b - r1, r1)
   circle(-r1, b - r1, r1)
   circle(0, r2, r2, :green)
   circle(a - r2, r2, r2, :green)
   circle(r2 - a, r2, r2, :green)
   circle(0, 2r2 + r3, r3, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r1, b - r1, "大円:r1,(r1,b-r1)", :red, :center, delta=-delta)
       point(a - r2, r2, "中円:r2,(a-r2,r2)", :green, :center, delta=-delta)
       point(0, 2r2 + r3, " 小円:r3,(0,2r2+r3)", :blue, :left, :vcenter)
       hline!([0], color=:gray, lw=0.5)
       vline!([0], color=:gray, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

 

算額（その91）

算額（その91）

千葉県君津市 神野寺 明治20(1887)年4月
http://www.wasan.jp/chiba/jinyaji4.html

円内に正三角形，甲円，乙円，丙円がある。それぞれの径を求めよ。

外円の半径を 1 とする。甲円の半径は 1/2 である。次いで，順に乙円，丙円の中心座標と径を求める。SymPy を使わあなくても求めることができる。

乙円の半径は 1/4 中心座標は cosd(i)*3/4, sind(i)*3/4, i = 30, 150, 270
丙円の半径は 1/6 中心座標は cosd(i+60)*2/3, sind(i+60)*2/3, i = 30, 150, 270

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; fontsize=10, mark=true)
  mark && scatter!([x], [y], color=color, markerstrokewidth=0)
  annotate!(x, y, text(string, fontsize, vertical, position, 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=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   circle(0, 0, 1, :black)
   plot!([0, cosd(-30), cosd(210), 0], [1, sind(-30), sind(210), 1], color=:magenta, lw=0.5)
   circle(0, 0, 1/2)
   for i = 30:120:270
       circle(cosd(i)*3/4, sind(i)*3/4, 1/4, :green)
       circle(cosd(i+60)*2/3, sind(i+60)*2/3, 1/6)
   end
   hline!([0], color=:black, lw=0.5)
   if more
       point(0, 0, "0 ", :red, :right)
       point(0, 1/2, "1/2 ", :red, :right)
       point(0, 2/3, "2/3 ", :blue, :right)
       point(cosd(30)*3/4, sind(30)*3/4, "(x,y)\nx=cosd(30)*3/4\ny=sind(30)*3/4", :green, :top)
       vline!([0], color=:black, lw=0.5)
   end
end;