裏 RjpWiki

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

算額(その2059)

2024年08月31日 | Julia

算額(その2059)

百二十六 群馬県倉渕村水沼 蓮華院 明治11年(1878)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:円3個,正方形,長方形

長方形の中に大円 1 個,小円 2 個,正方形 1 個を容れる。小円の直径が 3.432 寸のとき,大円の直径はいかほどか。

大円の半径と中心座標を r1, (-r1, 0)
小円の半径と中心座標を r2, (x2, r1 - r2)
とおき,以下の連立方程式を解く。
長方形の長辺と短辺は 4r1, 2r1 である。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, x2::positive;
eq1 = dist2(0, 0, r1, r1, x2, r1 - r2, r2)
eq2 = (x2 + r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
res = solve([eq1, eq2], (r1, x2))[3]  # 3 of 3

   (r2*(sqrt(2) + 3/2), r2/2)

大円の半径 r1 は,小円の半径 r2 の (√2 + 3/2) 倍である。
小円の直径が 3.432余りのとき,大円の直径は 3.432*(√2 + 3/2) = 10.001580946064461 で,ほぼ 10 寸と言いたいのだろう。

整数解にこだわるなら,小円の直径が 1562 寸のとき,大円の直径が 4552.001584426775 になるというのがある。

3.432*(√2 + 3/2)

   10.001580946064461

1562*(√2 + 3/2)

   4552.001584426775

function draw(r2, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, x2) = r2 .* (√2 + 3/2, 1/2)
   @printf("小円の直径が %g のとき,大円の直径は %g である。\n", 2r2, 2r1)
   @printf("r1 = %g;  r2 = %g; x2 = %g\n", r1, r2, x2)
   plot(r1 .* [2, 2, -2, -2, 2], r1 .* [-1, 1, 1, -1, -1], color=:green, lw=0.5)
   plot!(r1 .* [0, 1, 2, 1, 0], r1 .* [0, -1, 0, 1, 0], color=:blue, lw=0.5)
   circle(-r1, 0, r1)
   circle22(x2, r1 - r2, r2, :magenta)
   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, 0, "大円:r1,(-r1,0)", :red, :center, delta=-2delta)
       point(r1, 0, "(r1,0)", :blue, :center, delta=-2delta)
       point(x2, r1 - r2, "小円:r2\n(x2,r1-r2)", :magenta, :center, :bottom, delta=delta/2)
       point(2r1, r1, "(2r1,r1)", :green, :right, :bottom, delta=delta/2)
   end
end;

draw(3.432/2, true)

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

算額(その2058)

2024年08月31日 | Julia

算額(その2058)

百二十二 群馬県藤岡市東平井 諏訪神社 明治7年(1874)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:円5個,斜線2本

斜線 2 本で挟まれた上円,中円,下円とそれらが交差する領域に 2 個の等円を容れる。上円,中円,下円の直径がそれぞれ 3 寸,4 寸,5.4 寸のとき,等円の直径はいかほどか。

斜線(を延長したとき)の交点座標を (a, 0)
等円の半径と中心座標を r4, (x12, 0), (x23, 0); x12 = 2r3 + 2r2 - 3r4; x23 = 2r3 - r4
上円の半径と中心座標を r1, (x1, 0); x1 = 2r3 + 2r2 + r1 - 4r4
中円の半径と中心座標を r2, (x2, 0); x2 = 2r3 + r2 - 2r4
下円の半径と中心座標を r3, (x3, 0); x3 = r3
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive,
     r3::positive, r4::positive,
     x1::positive, x2::positive, x3::positive;
x1 = 2r3 + 2r2 + r1 - 4r4
x2 = 2r3 + r2 - 2r4
x3 = r3
eq1 = r1/(a - x1) - r2/(a - x2)
eq2 = r1/(a - x1) - r3/(a - x3)
res = solve([eq1, eq2], (r4, a));

等円の半径は (r1*r3 - r2^2)/(r1 - 2*r2 + r3) である。
上円,中円,下円の直径がそれぞれ 3 寸,4 寸,5.4 寸のとき,等円の直径は 0.5 寸である。

res[r4] |> println
res[r4](r1 => 3//2, r2 => 4//2, r3 => 54//20) |> println

   (r1*r3 - r2^2)/(r1 - 2*r2 + r3)
   1/4

斜線の交点座標は (18.9, 0) である。

res[a] |> println
res[a](r1 => 3//2, r2 => 4//2, r3 => 54//20) |> println

   2*r3*(-r2 + r3)/(r1 - 2*r2 + r3)
   189/10

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

   r1 = 1.5;  r2 = 2;  r3 = 2.7; , r4 = 0.25;  a = 18.9
   x1 = 9.9;  x2 = 6.9;  x3 = 2.7;  x12 = 8.65;  x23 = 5.15

function draw(r1, r2, r3, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r4, a) = (1/4, 189/10)
   x1 = 2r3 + 2r2 + r1 - 4r4
   x2 = 2r3 + r2 - 2r4
   x3 = r3
   x12 = 2r3 + 2r2 - 3r4
   x23 = 2r3 - r4
   θ = asind(r1/(a - x1))
   slope = tand(θ)
   @printf("等円の直径は %g である。\n", 2r4)
   @printf("r1 = %g;  r2 = %g;  r3 = %g; , r4 = %g;  a = %g\n", r1, r2, r3, r4, a)
   @printf("x1 = %g;  x2 = %g;  x3 = %g;  x12 = %g;  x23 = %g\n", x1, x2, x3, x12, x23)
   plot()
   circle(x1, 0, r1)
   circle(x2, 0, r2, :blue)
   circle(x3, 0, r3, :green)
   circle(x12, 0, r4, :orange)
   circle(x23, 0, r4, :orange)
   abline(a, 0, slope, 0, 1.05a)
   abline(a, 0, -slope, 0, 1.05a)
   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(x1, 0, " 上円:r1,(x1,0)", :red, :left, :vcenter)
       point(x2, 0, "中円:r2\n(x2,0)", :blue, :center, delta=-3delta)
       point(x3, 0, "下円:r3\n(x3,0)", :green, :center, delta=-3delta)
       point(a, 0, "a", :black, :center, delta=-3delta)
   end
end;

draw(3/2, 4/2, 5.4/2, true)

 

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

算額(その2057)

2024年08月31日 | Julia

算額(その2057)

百二十二 群馬県藤岡市東平井 諏訪神社 明治7年(1874)
群馬県和算研究会:群馬の算額,上武印刷株式会社,高崎市,1987年3月31日.
キーワード:円5個,正方形

正方形の中に,中円,東円,西円,南円,北円の 5 個を容れる。東円,西円,南円,北円の直径が 7 寸,6 寸,3 寸,9 寸 のとき,正方形の一辺はいかほどか。

正方形の一辺の長さを a
中円の半径と中心座標を r1, (x1, y1)
東円の半径と中心座標を r2, (a - r2, a - r2)
西円の半径と中心座標を r3, (r3, r3)
南円の半径と中心座標を r4, (a - r4, r4)
北円の半径と中心座標を r5, (r5, a - r5)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, x1::positive, y1::positive,
     r2::positive, r3::positive, r4::positive, r5::positive;
eq1 = (a - r2 - x1)^2 + (a - r2 - y1)^2 - (r1 + r2)^2
eq2 = (x1 - r3)^2 + (y1 - r3)^2 - (r1 + r3)^2
eq3 = (a - r4 - x1)^2 + (y1 - r4)^2 - (r1 + r4)^2
eq4 = (x1 - r5)^2 + (a - r5 - y1)^2 - (r1 + r5)^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, r1, x1, y1) = u
   return [
-(r1 + r2)^2 + (a - r2 - x1)^2 + (a - r2 - y1)^2,  # eq1
-(r1 + r3)^2 + (-r3 + x1)^2 + (-r3 + y1)^2,  # eq2
-(r1 + r4)^2 + (-r4 + y1)^2 + (a - r4 - x1)^2,  # eq3
-(r1 + r5)^2 + (-r5 + x1)^2 + (a - r5 - y1)^2,  # eq4
   ]
end;

(r2, r3, r4, r5) = [7, 6, 3, 9] ./ 2
iniv = BigFloat[20, 8, 13, 8]

res = nls(H, ini=iniv)

   ([21.0, 7.625, 12.625, 7.5], true)

東円,西円,南円,北円の直径が 7 寸,6 寸,3 寸,9 寸 のとき,正方形の一辺は 21 寸である。

しかし,残念ながら,この算額の問題は欠陥がある。
正方形の一辺の長さは 21 寸になるが,中円が正方形からはみ出す。

北円の直径を 8.9 寸程度にすると,正方形の一辺の長さは19.136 寸程度になり,「群馬の算額」の 140 ページにあるような図になる。

このような齟齬は,得られた解に基づいて正確な図を描いて検証するということがされていなかったことから生まれるのだろう。また,「群馬の算額」も含めて,第三者的立場でのチェックができていないということも根底にある。

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

  r2 = 3.5;  r3 = 3; , r4 = 1.5;  r5 = 4.45;  a = 19.1362;  r1 = 6.37033;  x1 = 11.6481;  y1 = 6.60738

function draw(r2, r3, r4, r5, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r1, x1, y1) = res[1]
   @printf("東円,西円,南円,北円の直径が %g,%g,%g,%g のとき,正方形の一辺は %g である。\n",
       2r2, 2r3, 2r4, 2r5, a)
   @printf("r2 = %g;  r3 = %g; , r4 = %g;  r5 = %g;  a = %g;  r1 = %g;  x1 = %g;  y1 = %g\n",
       r2, r3, r4, r5, a, r1, x1, y1)
   plot([0, a, a, 0,  0], [0, 0, a, a, 0], color=:green, lw=0.5)
   circle(x1, y1, r1)
   circle(a - r2, a - r2, r2, :blue)
   circle(r3, r3, r3, :magenta)
   circle(a - r4, r4, r4, :orange)
   circle(r5, a - r5, r5, :brown)
   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(a, a, "(a,a)", :green, :right, :bottom, delta=delta/2)
       point(x1, y1, "中円:r1,(x1,y1)", :red, :center, delta=-delta/2)
       point(a - r2, a - r2, "東円:r2,(a-r2,a-r2)", :blue, :center, delta=-delta/2)
       point(r3, r3, "西円:r3,(r3,r3)", :magenta, :center, delta=-delta/2)
       point(a - r4, r4, "南円:r4,(a-r4,r4)", :black, :right,:bottom, delta=delta/2, deltax=5delta)
       point(r5, a - r5, "北円:r5,(r5,a-r5)", :brown, :center, delta=-delta/2)
   end
end;

draw(7/2, 6/2, 3/2, 9/2, true)

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

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

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