裏 RjpWiki

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

算額(その1087)

2024年06月22日 | Julia

算額(その1087)

二十六 一関市萩荘 赤萩観音寺 弘化4年(1847)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html

キーワード:円3個,正三角形,二等辺三角形

正三角形と二等辺三角形が重なってできる領域に等円を 3 個容れる。等円の直径が 1 寸のとき正三角形の一辺の長さはいかほどか。

正三角形の一辺の長さを 2a
等円の半径と中心座標を r, (x, r), (0, -r)
二等辺三角形の底辺の頂点座標を (x0, -2r)
とおいて以下の連立方程式を解く。
しかし,SymPy では一度に解けない。さらに,eq1, eq3 から x,x0 を求めて eq2 に代入し,方程式を解くも,a = √3r/3 という不適切回しかでてこない。方程式のグラフを描くと a = √3 という解があるのを確認できるのに不思議だ。
解析解を諦めて数値解を求めると解は求まる。

しかし,ここでめげてはいけない。初心に帰って図に示すような補助線を描くと,答えが見えてくる。



できる正三角形は全て合同で,一辺の長さは内接円の直径の √3 倍である。
したがって,等円の直径が 1 寸のとき,元の正三角形の一辺の長さは 2√3 = 3.4641016151377544 寸である。

include("julia-source.txt")

using SymPy
@syms a::positive, r::positive, x::positive, x0::positive
eq1 = dist2(a, 0, 0,√Sym(3)a, x, r, r)/√Sym(3)
eq2 = dist(x0, -2r, 0, √Sym(3)a, x, r)- r^2 |> expand |> simplify |> numerator |> factor
eq3 = x0*(√Sym(3)a + r) - r*sqrt(x0^2 + (√Sym(3)a + 2r)^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=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (a, x, x0) = u
   return [
       sqrt(3)*(3*a^2/4 - sqrt(3)*a*r/2 - 3*a*x/2 - 3*r^2/4 + sqrt(3)*r*x/2 + 3*x^2/4)/3,  # eq1
       -27*a^6*r^2 + 27*a^6*x^2 - 54*a^6*x*x0 + 27*a^6*x0^2 - 108*sqrt(3)*a^5*r^3 + 108*sqrt(3)*a^5*r*x^2 - 162*sqrt(3)*a^5*r*x*x0 + 54*sqrt(3)*a^5*r*x0^2 - 540*a^4*r^4 + 540*a^4*r^2*x^2 - 540*a^4*r^2*x*x0 + 54*a^4*r^2*x0^2 + 18*a^4*x^2*x0^2 - 36*a^4*x*x0^3 + 18*a^4*x0^4 - 480*sqrt(3)*a^3*r^5 + 480*sqrt(3)*a^3*r^3*x^2 - 240*sqrt(3)*a^3*r^3*x*x0 - 96*sqrt(3)*a^3*r^3*x0^2 + 48*sqrt(3)*a^3*r*x^2*x0^2 - 60*sqrt(3)*a^3*r*x*x0^3 + 12*sqrt(3)*a^3*r*x0^4 - 720*a^2*r^6 + 720*a^2*r^4*x^2 - 288*a^2*r^4*x0^2 + 144*a^2*r^2*x^2*x0^2 - 72*a^2*r^2*x*x0^3 - 27*a^2*r^2*x0^4 + 3*a^2*x^2*x0^4 - 6*a^2*x*x0^5 + 3*a^2*x0^6 - 192*sqrt(3)*a*r^7 + 192*sqrt(3)*a*r^5*x^2 + 96*sqrt(3)*a*r^5*x*x0 - 96*sqrt(3)*a*r^5*x0^2 + 64*sqrt(3)*a*r^3*x^2*x0^2 + 16*sqrt(3)*a*r^3*x*x0^3 - 20*sqrt(3)*a*r^3*x0^4 + 4*sqrt(3)*a*r*x^2*x0^4 - 2*sqrt(3)*a*r*x*x0^5 - 2*sqrt(3)*a*r*x0^6 - 64*r^8 + 64*r^6*x^2 + 64*r^6*x*x0 - 32*r^6*x0^2 + 32*r^4*x^2*x0^2 + 32*r^4*x*x0^3 - 4*r^4*x0^4 + 4*r^2*x^2*x0^4 + 4*r^2*x*x0^5,  # eq2
       -r*sqrt(x0^2 + (sqrt(3)*a + 2*r)^2) + x0*(sqrt(3)*a + r),  # eq3
   ]
end;

r = 1/2
iniv = BigFloat[1.7, 0.87, 0.58]
res = nls(H, ini=iniv)

   ([1.7320508075688774, 0.8660254037844387, 0.5773502691896257], true)

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   (a, x, x0) = [1.7320508075688774, 0.8660254037844387, 0.5773502691896257]
   #@printf("小円の直径が %g のとき,中円の直径は %g である。\n", 2r4, 2r3)
   #@printf("r4 = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  y3 = %g\n", r4, r1, r2, r3, x3, y3)
   plot([a, 0, -a, a], [0, √3a, 0, 0], color=:blue, lw=0.5)
   plot!([x0, 0, -x0, x0], [-2r, √3a, -2r, -2r], color=:green, lw=0.5)
   plot!([0, x, -x, 0], [0, √3x, √3x, 0], color=:blue, lw=1)
   circle2(x, r, r)
   circle(0, -r, r)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(x, r, "等円:r,(x,r)", :red, :center, delta=-delta/2)
       point(0, -r, "等円:r,(0,-r)", :red, :center, delta=-delta/2)
       point(a, 0, "a", :blue, :left, :bottom, delta=delta/2)
       point(x0, -2r, " (x0,-2r)", :green, :left, :vcenter)
   end
end;

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

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

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