裏 RjpWiki

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

算額(その827)

2024年03月30日 | Julia

算額(その827)

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子,谷垣美保,萬伸介:瀬峰泉谷熊野神社奉納算額をめぐる諸問題,仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)

https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

全円の中に水平な弦を引き,その上に正三角形を置き,その内部に内接する円を入れる。弦の下に三角形の内部の円と同じ大きさの円(等円)が弦に接し互いにも接し,さらに両端の円は全円とも接している。
等円の直径が 6 寸のとき,全円の直径はいかほどか。

全円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (0, a + r)
弦と y 軸の交点座標を (0, a)
とおき,以下の連立方程式を解く。

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

using SymPy

@syms R::positive, r::positive, a::positive
eq1 = R - a - 3r  # 正三角形の高さは内接する円の半径の 3 倍である
eq2 = (2r)^2 + (a - r)^2 - (R -r)^2  # 両端の等円が全円に内接する
res = solve([eq1, eq2], (R, a))[1]

   (19*r/6, r/6)

全円の半径 R は,等円の半径 r の 19/6 倍である。
等円の直径が 6 寸であれば,全円の直径は 19 寸である。

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 6/2
   (R, a) = (19*r/6, r/6)
   @printf("全円の直径 = %g\n", 2R)
   plot(√3r .* [1, 0, -1, 1], a .+ [0, 3r, 0, 0], color=:blue, lw=0.5)
   circle(0, 0, R, :black)
   circle(0, a + r, r)
   circle(0, a - r, r, :green)
   circle2(2r, a - r, r, :green)
   x = sqrt(R^2 - a^2)
   segment(-x, a, x, a, :magenta)
   if more == true
       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(0, a, " a", :magenta, :center, :bottom, delta=delta)
       point(0, a + r, " 等円:r,(0,a+r)", :red, :center, :bottom, delta=delta)
       point(0, a - r, " 等円:r,(0,a-r)", :red, :center, :bottom, delta=delta)
       point(2r, a - r, " 等円:r,(2r,a-r)", :red, :center, :bottom, delta=delta)
       point(0, R, " R", :black, :left, :bottom, delta=delta/2)
   end
end;

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

算額(その826)

2024年03月30日 | Julia

算額(その826)

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子,谷垣美保,萬伸介:瀬峰泉谷熊野神社奉納算額をめぐる諸問題,仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)

https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

正方形の中に四分円と正方形と小円が入っている。小円の直径が 5 寸のとき,内部の正方形の一辺の長さはいかほどか。

四分円の半径と中心座標を r1, (0, 0), (r1, 0)
小円の半径と中心座標を r2, (r1/2, r1 - r2)
内部の正方形の一辺の長さを a
とおき,以下の連立方程式を解く。

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

using SymPy

@syms r1::positive, r2::positive, a::positive
eq1 = (r1/2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (r1/2 + a/2)^2 + a^2 - r1^2;
res = solve([eq1, eq2], (r1, a))[1]

   (16*r2, 48*r2/5)

正方形の一辺の長さ a は小円の半径 r2 の 48/5 倍である。
小円の直径が 5 寸のとき,正方形の一辺の長さは 24 寸である。
ちなみに,外の正方形の一辺の長さは 40 寸である。

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 5/2
   (r1, a) = (16*r2, 48*r2/5)
   @printf("正方形の一辺の長さ = %g\n", a)
   plot(r1 .* [0, 1, 1, 0, 0], r1 .* [0,0,1,1,0], color=:green, lw=0.5)
   circle(r1/2, r1 - r2, r2)
   plot!((r1 - a)/2 .+ [0, a, a, 0, 0], [0, 0, a, a, 0], color=:magenta, lw=0.5)
   circle(0, 0, r1, :blue, beginangle=0, endangle=90)
   circle(r1, 0, r1, :blue, beginangle=90, endangle=180)
   if more == true
       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 - a)/2, 0, " (r1-a)/2", :magenta, :left, :bottom, delta=delta/2)
       point((r1 + a)/2, 0, " (r1+a)/2", :magenta, :left, :bottom, delta=delta/2)
       point(r1/2, 0, "r1/2", :magenta, :center, :bottom, delta=delta/2)
       point(r1/2, a, "(r1/2,a)", :magenta, :center, :bottom, delta=delta/2)
       point(r1, 0, " r1", :blue, :left, :bottom, delta=delta/2)
       point(0, r1, " r1", :blue, :left, :bottom, delta=delta/2)
       point(r1/2, r1 - r2, " 小円:r1,(r1/2,r1-r2", :black, :center, :bottom, delta=delta)
   end
end;

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

算額(その825)

2024年03月30日 | Julia

算額(その825)

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子,谷垣美保,萬伸介:瀬峰泉谷熊野神社奉納算額をめぐる諸問題,仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)

https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

直線上に甲円と乙円が載っており互いに接している。乙円の上に正方形が載り,正方形の二辺は甲円に接している。乙円の直径が 1 寸のとき,甲円の直径はいかほどか。

甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (0, r2)
正方形の大きさはこの問題では無関係であるが,右端の頂点座標を例えば (x1, 2r2 + x1)
とおき,以下の連立方程式を解く。

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

using SymPy

@syms r1::positive, x1::positive, r2::positive
eq1 = x1^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = dist(0, 2r2, x1, 2r2 + x1, x1, r1) - r1^2;
res = solve([eq1, eq2], (r1, x1))[1]

   (2*r2, 2*sqrt(2)*r2)

甲円の半径 r1 は乙円の半径 r2 の 2 倍である。
乙円の直径が 1 寸のとき,甲円の直径は 2 寸である。

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (r1, x1) = (2*r2, 2*sqrt(2)*r2)
   @printf("甲円の直径 = %g\n", 2r1)
   plot()
   circle(0, r2, r2, :green)
   circle2(x1, r1, r1)
   plot!([0, x1, 0, -x1, 0], 2r2 .+ [0, x1, 2x1, x1, 0], color=:blue, lw=0.5)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0, r1], color=:gray80, lw=0.5)
       vline!([0, x1], color=:gray80, lw=0.5)
       point(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=-delta/2)
       point(0, r2, "乙円:r2\n(0,r2)", :green, :center, delta=-delta/2)
       point(x1, 2r2 + x1, " (x1,2r2+x1)", :blue, :left, :vcenter)
   end
end;

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

算額(その824)

2024年03月30日 | Julia

算額(その824)

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子,谷垣美保,萬伸介:瀬峰泉谷熊野神社奉納算額をめぐる諸問題,仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)

https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

全円の中に 2 本の斜線と,圭(二等辺三角形),甲円 3 個,乙円 2 個を入れる。乙円の直径が 3 寸のとき,甲円の直径はいかほどか。

全円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1), (x1, R - 3r1)
乙円の半径と中心座標を r2, (x2, R - 2r1 +r2)
斜線と円の交点座標を (x0, y0)
とおき,以下の連立方程式を解く。

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

using SymPy

@syms R::positive, r1::positive, x1::positive,
     r2::positive, x2::negative, x0, y0
eq1 = x2^2 + (R - 2r1 +r2)^2 - (R - r2)^2
eq2 = dist(x0, y0, 0, -R, x1, R - 3r1) - r1^2
eq3 = dist(x0, y0, 0, -R, 0, R - r1) - r1^2
eq4 = dist(x0, y0, 0, -R, x2, R - 2r1 + r2) - r2^2
eq5 = dist(sqrt(R^2 - (R - 2r1)^2), R - 2r1, 0, -R, x1, R - 3r1) - r1^2
eq6 = x0^2 + y0^2 - R^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)
   (R, r1, x1, x2, x0, y0) = u
   return [
       x2^2 - (R - r2)^2 + (R - 2*r1 + r2)^2,  # eq1
       -r1^2 + (-x0 + x0*(-x0*(-x0 + x1) + (-R - y0)*(R - 3*r1 - y0))/(x0^2 + (-R - y0)^2) + x1)^2 + (R - 3*r1 - y0 - (-R - y0)*(-x0*(-x0 + x1) + (-R - y0)*(R - 3*r1 - y0))/(x0^2 + (-R - y0)^2))^2,  # eq2
       -r1^2 + (x0*(x0^2 + (-R - y0)*(R - r1 - y0))/(x0^2 + (-R - y0)^2) - x0)^2 + (R - r1 - y0 - (-R - y0)*(x0^2 + (-R - y0)*(R - r1 - y0))/(x0^2 + (-R - y0)^2))^2,  # eq3
       -r2^2 + (-x0 + x0*(-x0*(-x0 + x2) + (-R - y0)*(R - 2*r1 + r2 - y0))/(x0^2 + (-R - y0)^2) + x2)^2 + (R - 2*r1 + r2 - y0 - (-R - y0)*(-x0*(-x0 + x2) + (-R - y0)*(R - 2*r1 + r2 - y0))/(x0^2 + (-R - y0)^2))^2,  # eq4
       -r1^2 + (-r1 - (-2*R + 2*r1)*(-r1*(-2*R + 2*r1) - sqrt(R^2 - (R - 2*r1)^2)*(x1 - sqrt(R^2 - (R - 2*r1)^2)))/(R^2 + (-2*R + 2*r1)^2 - (R - 2*r1)^2))^2 + (x1 + sqrt(R^2 - (R - 2*r1)^2)*(-r1*(-2*R + 2*r1) - sqrt(R^2 - (R - 2*r1)^2)*(x1 - sqrt(R^2 - (R - 2*r1)^2)))/(R^2 + (-2*R + 2*r1)^2 - (R - 2*r1)^2) - sqrt(R^2 - (R - 2*r1)^2))^2,  # eq5
       -R^2 + x0^2 + y0^2,  # eq6
   ]
end;

r2 = 3/2
iniv = BigFloat[7.9, 2.1, 3.5, 3.5, 2.3, 7.7]
res = nls(H, ini=iniv)

   ([8.0, 2.0, 3.4641016151377544, 3.4641016151377544, 2.262270442538942, 7.673469387755102], true)

乙円の直径が 3 寸のとき,甲円の直径は 4 寸である。

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

R = 8;  r1 = 2;  x1 = 3.4641;  x2 = 3.4641;  x0 = 2.26227;  y0 = 7.67347

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 3//2
   (R, r1, x1, x2, x0, y0) = res[1]
   @printf("甲円の直径 = %g\n", 2r1)
   @printf("R = %g;  r1 = %g;  x1 = %g;  x2 = %g;  x0 = %g;  y0 = %g\n", R, r1, x1, x2, x0, y0)
   plot()
   circle(0, 0, R, :green)
   circle(0, R - r1, r1)
   circle2(x2, R - 2r1 + r2, r2, :blue)
   circle2(x1, R - 3r1, r1)
   plot!([-x0, 0, x0], [y0, -R, y0], color=:magenta, lw=0.5)
   y00 = R - 2r1
   x00 = sqrt(R^2 - y00^2)
   plot!([-x00, 0, x00, -x00], [y00, -R, y00, y00], color=:magenta, lw=0.5)
   if more == true
       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(0, R - r1, " 甲円:r1\n(0,R-r1)", :red, :center, :bottom, delta=delta/2)

       point(x1, R - 3r1, " 甲円:r1\n(x1,R-3r1)", :red, :center, delta=-delta/2)
       point(x2, R - 2r1 + r2, " 乙円:r2\n(x2,R-2r1+r2)", :blue, :center, :bottom, delta=delta/2)
       point(x0, y0, "(x0,y0)", :magenta, :left, :bottom, delta=delta/2)
       point(sqrt(R^2 - (R - 2r1)^2), R - 2r1, "sqrt(R^2-(R-2r1)^2),R-2r1 ", :black, :right, delta=-delta/2)
       point(0, R, " R", :green, :left, :bottom, delta=delta/2)
   end
end;

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

算額(その823)

2024年03月30日 | Julia

算額(その823)

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子,谷垣美保,萬伸介:瀬峰泉谷熊野神社奉納算額をめぐる諸問題,仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)

https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

全円の中に 2 本の斜線と,大円 1 個,小円 3 個を入れる。 大円と小円は,水平な弦に接している。小円の直径が 1 寸のとき,大円の直径はいかほどか。

全円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, r1 - R)
小円の半径と中心座標を r2, (0, r2 + 2r1 - R), (x2, 2r1 - R - r2)
とおき,以下の連立方程式を解く。

「答」,「術」は水平な弦が全円の中心を通る(つまり直径)と考えているようであるが,実際は水平な弦と y 軸の交点座標は (0, 2r1 - R) である。そして,r1 は R/2 ではない。

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

using SymPy

@syms R::positive, r1::positive, r2::positive,
     x2::positive, x0::positive, y0::negative
sinθ = x0/sqrt(x0^2 + (R - y0)^2)
eq1 = (2R - r1)*sinθ - r1
eq2 = (2R - 2r1 - r2)*sinθ - r2
eq3 = dist(0, R, x0, y0, x2, 2r1 - R - r2) - r2^2
eq4 = x2^2 + (2r1 - R - r2)^2 - (R - r2)^2
eq5 = x0^2 + y0^2 - R^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)
   (R, r1, x2, x0, y0) = u
   return [
       -r1 + x0*(2*R - r1)/sqrt(x0^2 + (R - y0)^2),  # eq1
       -r2 + x0*(2*R - 2*r1 - r2)/sqrt(x0^2 + (R - y0)^2),  # eq2
       -r2^2 + (-x0*(x0*x2 + (-R + y0)*(-2*R + 2*r1 - r2))/(x0^2 + (-R + y0)^2) + x2)^2 + (-2*R + 2*r1 - r2 - (-R + y0)*(x0*x2 + (-R + y0)*(-2*R + 2*r1 - r2))/(x0^2 + (-R + y0)^2))^2,  # eq3
       x2^2 - (R - r2)^2 + (-R + 2*r1 - r2)^2,  # eq4
       -R^2 + x0^2 + y0^2,  # eq5
   ]
end;

r2 = 1/2
iniv = BigFloat[2, 1, 1.4, 1.27, -1.6]
res = nls(H, ini=iniv)

   ([2.00016219313451, 1.0090868249613552, 1.4206243873461804, 1.2703923742568435, -1.5449116525791091], true)

大円の直径は 2.0181736499227103 である。2 ではない。

問には「乃小円者請至多数」とある。小円は3個で,直径が 1 寸と言っているので,この意味がわからなかった。
なお,小円,大円の直径がそれぞれ 1 寸,2 寸のときには,図のイメージはかなり異なるものになる。

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

R = 2.00016;  r1 = 1.00909;  r2 = 0.5;  x2 = 1.42062;  x0 = 1.27039;  y0 = -1.54491

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1//2
   (R, r1, x2, x0, y0) = res[1]
   @printf("大円の直径 = %g\n", 2r1)
   @printf("R = %g;  r1 = %g;  r2 = %g;  x2 = %g;  x0 = %g;  y0 = %g\n", R, r1, r2, x2, x0, y0)
   plot()
   circle(0, 0, R, :green)
   circle(0, r1 - R, r1)
   circle(0, r2 + 2r1 - R, r2, :blue)
   circle2(x2, 2r1 - R - r2, r2, :blue)
   plot!([-x0, 0, x0], [y0, R, y0], color=:black, lw=0.5)
   y = 2r1 - R
   x = sqrt(R^2 - y^2)
   segment(-x, y, x, y, :green)
   if more == true
       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(0, 2r1 - R, " 2r1-R", :green, :center, delta=-delta/2)
       point(x2, 2r1 - R - r2, "小円:r2\n(x2,2r1-R-r2)", :blue, :center, delta=-delta/2)
       point(0, r2 + 2r1 - R, "小円:r2\n(0,r2+2r1-R)", :blue, :center, delta=-delta/2)
       point(0, r1 - R, "大円:r1,(0,r1-R)", :red, :center, delta=-delta/2)
       point(0, -R, " -R", :red, :center, :bottom, delta=delta/2)
       point(R, 0, " R", :red, :left, :bottom, delta=delta/2)
   end
end;

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

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

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