裏 RjpWiki

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

算額(その740)

2024年03月01日 | Julia

算額(その740)

浅草観音堂の算額
会田安明:『他流諸国之表題集』
山口正義:やまぶき4,第84号

https://yamabukiwasan.sakura.ne.jp/ymbk84.pdf

直角三角形の中に,大円,中円,小円を入れる。
中円,小円の直径をそれぞれ 3.65767079 寸,1.95349271 寸としたとき,鈎,股,大円の直径を求めよ。
後述するが,1.92349271 の誤記である。

図形的には,算額(その453)と同じである。しかし,条件の付け方により,この問題は SymPy では解けず,NLsolve で数値解を求めるしかない。

直角三角形の直角を挟む二辺を「鈎」,「股」とする
大円の半径と中心座標を r1, (r1, r1)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (r3, y3)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1, r2, x2, r3, y3, 鈎, 股, 弦, d
eq1 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (r1 - r3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq3 = r1*(股 - x2) - r2*(股 - r1)
eq4 = r1*(鈎 - y3) - r3*(鈎 - r1)
eq5 = 鈎 + 股 - sqrt(鈎^2 + 股^2) - 2r1;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r1, 鈎, 股, x2, y3))

「問」で与えられた条件,中円直径 = 3.65767079,小円直径 = 1.95349271 では,「答」のように「きれいな数」は出ない。
r1, 鈎, 股 は 3.021927068553583, 8.099020725821061, 14.930767011037842 になる。

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)
   (r1, 鈎, 股, x2, y3) = u
   return [
       (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2,
       (r1 - r3)^2 + (y3 - r1)^2 - (r1 + r3)^2,
       r1*(股 - x2) - r2*(股 - r1),
       r1*(鈎 - y3) - r3*(鈎 - r1),
       鈎 + 股 - sqrt(鈎^2 + 股^2) - 2r1
   ]
end;

r2 = 3.65767079/2
r3 = 1.95349271/2  # 誤記のまま。正しくは 1.92349271 である
iniv = BigFloat[6/2, 8, 15, 7.7, 6.4]
res = nls(H, ini=iniv)

   ([3.021927068553583, 8.099020725821061, 14.930767011037842, 7.723674481624338, 6.458004046320843], true)

「答」がきれいな数になるように,中円,小円の直径をあれほど細かく書いたのだろうに。

「答」の方から逆算すると,小円直径 = 1.92349271 の誤記であることがわかる。
訂正して,計算し直すと r1, 鈎, 股 は 3.000046427350913, 8.000208871784782, 14.999849358214375 となり,3, 8, 15 というきれいな数になる。

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

r1 = 3.00005;  鈎 = 8.00021;  股 = 14.9998;  x2 = 7.68474;  y3=6.39727

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 3.65767079/2
   r3 = 1.92349271/2
   (r1, 鈎, 股, x2, y3) = res[1]
   @printf("r1 = %g;  鈎 = %g;  股 = %g;  x2 = %g;  y3=%g\n", r1, 鈎, 股, x2, y3)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:magenta, lw=0.5)
   circle(r1, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(r3, y3, r3, :green)
   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(r1, r1, "大円:r1,(r1,r1)", :red, :center, delta=-delta)
       point(x2, r2, "中円:r2,(x2,r2)", :blue, :center, delta=-delta)
       point(r3, y3, " 小円:r3,(r3,y3)", :green, :left, :vcenter)
       point(股, 0, "股", :magenta, :left, :bottom, delta=delta)
       point(0, 鈎, " 鈎", :magenta, :left, :bottom)
   end
end;

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

算額(その739)

2024年03月01日 | Julia

算額(その739)

千葉県富津市寺尾 六所神社 明治4年(1871)
山口正義:やまぶき4,第57号 

https://yamabukiwasan.sakura.ne.jp/ymbk57.pdf
山口正義:やまぶき4,第59号
https://yamabukiwasan.sakura.ne.jp/ymbk59.pdf

正五角形の内部に 5 本の矢を入れ,大きさの同じ二等辺三角形と小さな正五角形に区分する。二等辺三角形に内接する楕円の短径が 1 寸のとき,長径はいかほどか。

二等辺三角形と楕円の関係は以下の図のようになる。

以下では x 軸に平行な軸を長径 2a,y 軸上に短径 2b とする(実際は a 斜辺の長さを l
楕円と斜辺の交点座標を (x0, y0)
として以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms l, x, a, b, x0, y0
x = l*sind(Sym(18))
y0 = b*sqrt(1 - x0^2/a^2) + b
eq2 = -b^2*x0/(a^2*(y0 - b)) - 1/tand(Sym(162))
eq3 = (x - x0)/y0 - tand(Sym(18))
res = solve([eq2, eq3], (b, x0))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (sqrt(5)*sqrt((-16*a^2*l^2*(15 - 7*sqrt(5))^2 + (-40*a^2 + 16*sqrt(5)*a^2 - 25*l^2 + 11*sqrt(5)*l^2)^2)/(-40*a^2 + 16*sqrt(5)*a^2 - 25*l^2 + 11*sqrt(5)*l^2)^2)*(-40*a^2 + 16*sqrt(5)*a^2 - 25*l^2 + 11*sqrt(5)*l^2)/(4*l*sqrt(5 - 2*sqrt(5))*(15 - 7*sqrt(5))), 4*a^2*l*(15 - 7*sqrt(5))/(-40*a^2 + 16*sqrt(5)*a^2 - 25*l^2 + 11*sqrt(5)*l^2))

function draw(more=false)
   pyplot(size=(300, 300), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   l = 15
   a = 1//2
   x = l*sind(18)
   (b, x0) = (sqrt(5)*sqrt((-16*a^2*l^2*(15 - 7*sqrt(5))^2 + (-40*a^2 + 16*sqrt(5)*a^2 - 25*l^2 + 11*sqrt(5)*l^2)^2)/(-40*a^2 + 16*sqrt(5)*a^2 - 25*l^2 + 11*sqrt(5)*l^2)^2)*(-40*a^2 + 16*sqrt(5)*a^2 - 25*l^2 + 11*sqrt(5)*l^2)/(4*l*sqrt(5 - 2*sqrt(5))*(15 - 7*sqrt(5))), 4*a^2*l*(15 - 7*sqrt(5))/(-40*a^2 + 16*sqrt(5)*a^2 - 25*l^2 + 11*sqrt(5)*l^2))
   y0 = b*sqrt(1 - x0^2/a^2) + b
   @printf("楕円の長径 = %.10g;  a = %.10g;  b = %.10g\n", 2b, a, b)
   @printf("b = %.10g;  x0 = %.10g;  y0 = %.10g\n", b, x0, y0)
   plot([x, 0, -x, x], [0, l*cosd(18), 0, 0], color=:blue, lw=0.5)
   ellipse(0, b, a, b)
   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(x0, y0)
   end
end;

   楕円の長径 = 14.09985453;  a = 0.5;  b = 7.049927263

   b = 7.049927263;  x0 = 0.1066282377;  y0 = 13.93767977

三角形に内接する楕円の長径,短径を求めるには,「算法助術」の公式97 を使うことができる。
求まる数値 7.04992726271756 は上と一致する。

@syms x, h, a, b
x = 15sind(Sym(18))
h = sqrt(15^2 - x^2)
a = 1//2
b = solve(h*(2a)^2 - (h - 2b)*(2x)^2, b)[1]
b |> factor |> println
b.evalf() |> println

   -(-447*sqrt(2) + sqrt(10))*sqrt(sqrt(5) + 5)/240
   7.04992726271756

@printf("楕円の長径 = %.10g;  a = %.10g;  b = %.10g\n", 2b, a, b)

   楕円の長径 = 14.09985453;  a = 0.5;  b = 7.049927263

正五角形の一辺の長さが 15 寸のとき,楕円の長径は 14.09985453 である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   innerx = Vector{Float64}(undef, 6)
   innery = Vector{Float64}(undef, 6)
   outerx = Vector{Float64}(undef, 6)
   outery = Vector{Float64}(undef, 6)
   outerlength = 15/(2sind(36))
   innerlength = outerlength*sind(22.5)
   for i in 1:6
       θ = (72*i) - 18
       innerx[i] = innerlength*cosd(θ)
       innery[i] = innerlength*sind(θ)
       θ = (72*i) - 54
       outerx[i] = outerlength*cosd(θ)
       outery[i] = outerlength*sind(θ)
   end
   @printf("一辺の長さ = %g\n", outerx[5] - outerx[4])
   a = 1//2
   b = sqrt(12815212505 - 113904000*sqrt(5))/(60*sqrt(5 - 2*sqrt(5))*(112 + 113*sqrt(5)))
   @printf("楕円の長径 = %.8g;  a = %.8g;  b = %.8g\n", 2b, a, b)
   plot()
   for i in 1:5
       segment(outerx[i], outery[i], innerx[i+1], innery[i+1], :blue)
       #segment(innerx[i], innery[i], innerx[i+1], innery[i+1], :green)
       segment(outerx[i], outery[i], outerx[i+1], outery[i+1], :red)
   end
   ox = (outerx[1] + innerx[1])/2
   oy = -3.15
   l = sqrt(ox^2 + oy^2)
   θ = atand(oy/ox)
   for i in 0:4
       ellipse(l*cosd(72i + θ), l*sind(72i + θ), a, b, φ=72i, color=:green)
   end
end;

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

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

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