裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

算額(その56)

2022年12月09日 | Julia

算額(その56)

宮城県角田市横倉 愛宕神社 明治15年(1882)1月
http://www.wasan.jp/miyagi/yokokuraatago.html

外円の中に大小 8 個の円が入っている。それぞれの径を求めよ。

外円の半径を 1 として,図のように記号を定め,方程式を解く。


using SymPy
@syms or1::positive, ir1::positive
@syms or2::positive, ir2::positive
@syms or3::positive, oy3::positive, ir3::positive, ix3::positive, iy3::positive
eq1 = or3^2 + (oy3 - ir1)^2 - (ir1 + or3)^2;
eq2 = ix3^2 + (iy3 - or1)^2 - (or1 + ir3)^2;
eq3 = or3^2 + (2 - ir2 - oy3)^2 - (ir2 + or3)^2;
eq4 = ix3^2 + (2 - or2 - iy3)^2 - (or2 + ir3)^2;
eq5 = or3^2 + (oy3 - 1)^2 -  (1 - or3)^2;
eq6 = ix3^2 + (iy3 - 1)^2 - (1 - ir3)^2;
eq7 = or1 + or2 - 1;
eq8 = (oy3 - 1)/or3 - (iy3 - 1)/ix3;
eq9 = (ix3 - or3)^2 + (iy3 - oy3)^2 - (or3 - ir3)^2;

方程式が複雑になるとどうも solve() では解けないこともでてくるようだ。
そこで,nlsolve() を使うことにする。

# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9])

eq1 |> expand |> println
eq2 |> expand |> println
eq3 |> expand |> println
eq4 |> expand |> println
eq5 |> expand |> println
eq6 |> expand |> println
eq7 |> expand |> println
eq8 |> expand |> println
eq9 |> expand |> println


   -2*ir1*or3 - 2*ir1*oy3 + oy3^2
   -ir3^2 - 2*ir3*or1 + ix3^2 + iy3^2 - 2*iy3*or1
   -2*ir2*or3 + 2*ir2*oy3 - 4*ir2 + oy3^2 - 4*oy3 + 4
   -ir3^2 - 2*ir3*or2 + ix3^2 + iy3^2 + 2*iy3*or2 - 4*iy3 - 4*or2 + 4
   2*or3 + oy3^2 - 2*oy3
   -ir3^2 + 2*ir3 + ix3^2 + iy3^2 - 2*iy3
   or1 + or2 - 1
   oy3/or3 - 1/or3 - iy3/ix3 + 1/ix3
   -ir3^2 + 2*ir3*or3 + ix3^2 - 2*ix3*or3 + iy3^2 - 2*iy3*oy3 + oy3^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)
  or1, ir1, or2, ir2, or3, oy3, ir3, ix3, iy3 = u
  return [
       -2*ir1*or3 - 2*ir1*oy3 + oy3^2,
       -ir3^2 - 2*ir3*or1 + ix3^2 + iy3^2 - 2*iy3*or1,
       -2*ir2*or3 + 2*ir2*oy3 - 4*ir2 + oy3^2 - 4*oy3 + 4,
       -ir3^2 - 2*ir3*or2 + ix3^2 + iy3^2 + 2*iy3*or2 - 4*iy3 - 4*or2 + 4,
       2*or3 + oy3^2 - 2*oy3,
       -ir3^2 + 2*ir3 + ix3^2 + iy3^2 - 2*iy3,
       or1 + or2 - 1,
       oy3/or3 - 1/or3 - iy3/ix3 + 1/ix3,
       -ir3^2 + 2*ir3*or3 + ix3^2 - 2*ix3*or3 + iy3^2 - 2*iy3*oy3 + oy3^2
       ];
end;

iniv = [89, 8, 49, 27, 58, 185, 40, 75, 200] ./ 135;

res = nls(h, ini=iniv)
names =["or1", "ir1", "or2", "ir2", "or3", "oy3", "ir3", "ix3", "iy3"]
println(res[2])
for (name, value) in zip(names, res[1])
   println("$name = $value")
end

   true
   or1 = 0.7023666314613604
   ir1 = 0.5412673937726713
   or2 = 0.2976333685386396
   ir2 = 0.17483505787654022
   or3 = 0.41809549294196335
   oy3 = 1.4047332629227232
   ir3 = 0.2642988189720262
   ix3 = 0.5285976379440508
   iy3 = 1.5117037863118916

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=false; fill=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   or1 = 0.7023666314613604
   ir1 = 0.5412673937726713
   or2 = 0.2976333685386396
   ir2 = 0.17483505787654022
   or3 = 0.41809549294196335
   oy3 = 1.4047332629227232
   ir3 = 0.2642988189720262
   ix3 = 0.5285976379440508
   iy3 = 1.5117037863118916
   circle(0, 1, 1, :black)
   circle(0, ir1, ir1, :blue)
   circle(0, or1, or1, :magenta)
   circle(0, 2-ir2, ir2, :brown)
   circle(0, 2-or2, or2, :red)
   circle(or3, oy3, or3, :cyan)
   circle(ix3, iy3, ir3, :green)
   circle(-or3, oy3, or3, :cyan)
   circle(-ix3, iy3, ir3, :green)
   if more
       point(0, 1, " 1", :black)
       point(0, ir1, " ir1", :blue)
       point(0, or1, " or1", :magenta)
       point(0, 2-ir2, " 2-ir2", :brown)
       point(0, 2-or2, " 2-or2", :red)
       point(or3, oy3, "(or3,oy3)", :cyan)
       point(ix3, iy3, "(ix3,iy3),ir3", :green, :center)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   end
end;

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村