裏 RjpWiki

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

算額(その645)

2024年01月19日 | Julia

算額(その645)

長野県小諸市八幡町 八幡神社 寛政11年(1799)
中村信弥「改訂増補 長野県の算額」

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

鈎股弦(直角三角形)内に内接円があり,中鈎(直角の頂点から斜辺(弦)への垂線)を引く。短弦と円の直径の積が 3.6 平方寸,中鈎と円の直径の積が 4.8 平方寸のとき,円の直径を求めよ。


直角三角形の 3 辺を鈎,股,弦,弦 = 短弦 + 長弦,中鉤と弦の交点座標を (x, y),円の直径を円径として以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
     中鈎::positive, 短弦::positive, 長弦::positive,
     円径::positive, x::positive, y ::positive

eq1 = 短弦^2 + 中鈎^2 - 鈎^2
eq2 = 短弦 * 円径 - 3.6
eq3 = 中鈎 * 円径 - 4.8
eq4 = 鈎 + 股 - 弦 - 円径
eq5 = 鈎^2 + 股^2 - 弦^2
eq6 = 長弦^2 + 中鈎^2 - 股^2
eq7 = 長弦 + 短弦 - 弦
eq8 = (x^2 + y^2) + (x^2 + (鈎 - y)^2) - 鈎^2
eq9 = (x^2 + y^2) + ((股 - x)^2 + y^2) - 股^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (円径, 鈎, 股, 弦, 中鈎, 短弦, 長弦, x, y))

   1-element Vector{NTuple{9, Sym}}:
    (2.00000000000000, 3.00000000000000, 4.00000000000000, 5.00000000000000, 2.40000000000000, 1.80000000000000, 3.20000000000000, 1.44000000000000, 1.92000000000000)

円の直径は 2 寸である。

その他のパラメータは以下の通りである。

鈎 = 3;  股 = 4;  弦 = 5;  中鈎 = 2.4;  短弦 = 1.8;  長弦 = 3.2;  x = 1.44;  y = 1.92

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (円径, 鈎, 股, 弦, 中鈎, 短弦, 長弦, x, y) = res[1]
   @printf("円径 = %g;  鈎 = %g;  股 = %g;  弦 = %g;  中鈎 = %g;  短弦 = %g;  長弦 = %g;  x = %g;  y = %g\n",
       円径, 鈎, 股, 弦, 中鈎, 短弦, 長弦, x, y)
   r = 円径/2
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(r, r, r, :blue)
   segment(0, 0, x, y, :red)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(x, y, " (x,y)", :red, :left, :vcenter)
       point(0, 鈎, " 鈎", :black, :left, :vcenter)
       point(股, 0, "股", :black, :left, :bottom, delta=delta/2)
       point(r, r, " (r,r)", :blue, :left, :vcenter)
       point(0.6, 2.7, "短弦", mark=false)
       point(2.8, 1.05, "長弦", mark=false)
       point(0.75, 1.5, "中鈎", mark=false)
   end
end;

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

算額(その644)

2024年01月19日 | Julia

算額(その644)

長野県軽井沢町峠 熊野神社 安政4年(1857)
中村信弥「改訂増補 長野県の算額」

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

二等辺三角形の中に大円,中円,荘園が入っている。中円は大円と等辺の接点を結ぶ直線に接している。また,小円は沖苑に内接し,大円に外接している。小円の直径が 1 寸のとき,中円の直径を求めよ。

二等辺三角形の底辺の長さを 2a,高さを h とする。また,大円と等辺の接点座標を (c, h) とする。
大円の半径と中心座標を r1, (0, r1)
中円の半径と中心座標を r2, (0, b + r2)
小円の半径と中心座標を r3, (0, 2r1 + r3)
とおき,以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r1::positive, r2::positive, r3::positive, a::positive, b::positive, c::positive, h::positive

r3 = 1//2
eq1 = r1/(h - r1) - a/sqrt(a^2 + h^2)
eq2 = r2/(h - b - r2) - a/sqrt(a^2 + h^2)
eq3 = c/(h - b) - a/h
eq4 = r2/r1 - a/h
eq5 = b^2 + (a - c)^2 - a^2
eq6 = b + 2r2 - (2r1 + 2r3)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 v, r.f_converged
end;

function H(u)
   (r1, r2, a, b, c, h) = u
   return [
       -a/sqrt(a^2 + h^2) + r1/(h - r1),  # eq1
       -a/sqrt(a^2 + h^2) + r2/(-b + h - r2),  # eq2
       -a/h + c/(-b + h),  # eq3
       -a/h + r2/r1,  # eq4
       -a^2 + b^2 + (a - c)^2,  # eq5
       b - 2*r1 + 2*r2 - 1,  # eq6
   ]
end;

iniv = BigFloat[1.67, 1.1, 2.42, 2.14, 2.2, 6.43]
res = nls(H, ini=iniv)

   (BigFloat[1.883203505913525864168947465362055090560951328672239179557958847642330007152124, 0.9999999999999999999999999999999999999999999999999999999915012754481286025568294, 3.132241882311900195655072338157300657084786186095832984118051754533183394146556, 2.766407011827051728337894930724110181121902657344478359132915144388402809190606, 1.663251938771469380287813597004943645093506732800541877736976628888590678029407, 5.898648894138951923992967268881410838206688843440311343250654955537926982258505], true)

中円の半径は 1寸(直径は 2 寸)である。
その他のパラメータは以下の通り。
r1 = 1.8832;  r2 = 1;  a = 3.13224;  b = 2.76641;  c = 1.66325;  h = 5.89865

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 1//2
   (r1, r2, a, b, c, h) = res[1]
   @printf("r1 = %g;  r2 = %g;  a = %g;  b = %g;  c = %g;  h = %g\n", r1, r2, a, b, c, h)
   plot([a, 0, -a, a], [0, h, 0, 0], color=:black, lw=0.5)
   circle(0, r1, r1, :blue)
   circle(0, b + r2, r2)
   circle(0, 2r1 + r3, r3, :green)
   segment(c, b, -c, b)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(c, b, " (c,b)", :black, :left, :vcenter)
       point(0, r1, " 大円:r1,(0, r1)", :blue, :left, :vcenter)
       point(0, b + r2, " 中円:r2\n (0,b+r2)", :red, delta=-2delta)
       point(0, 2r1 + r3, " 小円:r3\n (0,2r1+r3)", :black, :left, :vcenter)
       point(0, h, " h", :black, :left, :vcenter)
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
   end
end;

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

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

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