算額（その818）

算額（その818）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

長方形と，長方形の長辺を底辺として共有する二等辺三角形の中に，甲円，乙円，丙円を入れる。
丙円の直径が 1 寸のとき，乙円の直径はいかほどか。

甲円の半径と中心座標を r1, (r1, r1)
乙円の半径と中心座標を r2, (0, 2r1 + r2)
丙円の半径と中心座標を r3, (0, r3), (x3, y3)
長方形の長辺と短辺を 2r1, r1
二等辺三角形の高さを h
とおき，以下の連立方程式を解く。

乙円についての方程式 eq5 は，他の方程式とは独立なので，まずは eq1〜eq4 の連立方程式を解いて
h, r1, x3, y3 を求める。

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

using SymPy

@syms h::positive, r1::positive, r2::positive,
     r3::positive, x3::positive, y3::positive, d
r3 = 1//2
eq1 = r1^2 + (r1 - r3)^2 - (r1 + r3)^2
eq2 = (x3 - r1)^2 + (y3 - r1)^2 - (r1 - r3)^2
eq3 = (x3 - r1)*2r1 - h*(y3 - r1)
eq4 = dist(0, h, 2r1, 0, x3, y3) - r3^2
eq4 = numerator(apart(eq4, d));
(h, r1, x3, y3) = solve([eq1, eq2, eq3, eq4], (h, r1, x3, y3))[2]

   (2*sqrt(2*sqrt(7) + 8)/3 + 2*sqrt(7)*sqrt(2*sqrt(7) + 8)/3, 2, 3*sqrt(2*sqrt(7) + 8)/8 + 2, 3*sqrt(7)/8 + 13/8)

h, r1, x3, y3 が求められた前提で，次の方程式を解く。

eq5 = r2/(h - 2r1 - r2) - 2r1/sqrt(h^2 + 4r1^2)
res2 = solve(eq5, r2)[1]
res2 |> println

   -144/(36 + 3*sqrt(144 + (2*sqrt(2*sqrt(7) + 8) + 2*sqrt(7)*sqrt(2*sqrt(7) + 8))^2)) + 24*sqrt(2*sqrt(7) + 8)/(36 + 3*sqrt(144 + (2*sqrt(2*sqrt(7) + 8) + 2*sqrt(7)*sqrt(2*sqrt(7) + 8))^2)) + 24*sqrt(7)*sqrt(2*sqrt(7) + 8)/(36 + 3*sqrt(144 + (2*sqrt(2*sqrt(7) + 8) + 2*sqrt(7)*sqrt(2*sqrt(7) + 8))^2))

長い式は簡約化できる。

res2 |> factor |> println

   4*(-6 + sqrt(2)*sqrt(sqrt(7) + 4) + sqrt(14)*sqrt(sqrt(7) + 4))/(6 + sqrt(32*sqrt(7) + 128))

乙円の直径は 2.83398951148328 である。

2res2.evalf() |> println

   2.83398951148328

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

   h = 8.861;  r1 = 2;  r2 = 1.41699;  x3 = 3.36716;  y3 = 2.61716

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 1/2
   (h, r1, x3, y3) = (2*sqrt(2*sqrt(7) + 8)/3 + 2*sqrt(7)*sqrt(2*sqrt(7) + 8)/3, 2, 3*sqrt(2*sqrt(7) + 8)/8 + 2, 3*sqrt(7)/8 + 13/8)
   r2 = 4*(-6 + sqrt(2)*sqrt(sqrt(7) + 4) + sqrt(14)*sqrt(sqrt(7) + 4))/(6 + sqrt(32*sqrt(7) + 128))
   @printf("乙円の直径 = %g\n", 2r2)
   @printf("h = %g;  r1 = %g;  r2 = %g;  x3 = %g;  y3 = %g\n", h, r1, r2, x3, y3)
   plot([2r1, 2r1, -2r1, -2r1, 2r1], [0, 2r1, 2r1, 0, 0], color=:black, lw=0.5)
   circle2(r1, r1, r1)
   circle(0, 2r1 + r2, r2, :blue)
   circle(0, r3, r3, :green)
   circle2(x3, y3, r3, :green)
   plot!([-2r1, 0, 2r1], [0, h, 0], color=:orange, 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(x3, y3, "丙円:r3,(x3,y3) ", :black, :right, :vcenter)
       point(0, r3, " 丙円:r3,(0,r3)", :black, :left, :vcenter)
       point(r1, r1, "甲円:r1,(r1,r1)", :red, :center, delta=-delta/2)
       point(0, 2r1 + r2, "乙円:r2,(0,2r1+r2)", :blue, :center, delta=-delta/2)
       point(0, h, " h", :orange, :left, :vcenter)
       point(2r1, 2r1, "(2r1,2r1)", :black, :right, :bottom, delta=delta/2)
   end
end;

算額（その817）

算額（その817）

宮城県栗原市瀬峰泉谷 瀬峰泉谷熊野神社奉納算額
徳竹亜紀子，谷垣美保，萬伸介：瀬峰泉谷熊野神社奉納算額をめぐる諸問題，仙台高等専門学校名取キャンパス 研究紀要 第60号(2024)
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2024/03/kiyo2024-1.pdf

楕円内に等円 2 個が入っている。楕円の短径が 1 のとき，原点を点対称とする円と楕円の 2 接点間の距離を求めよ。

条件が「楕円の短径」のみなので，解けるかどうか不安であったが，「楕円の長径」は無関係であることがわかった。

楕円の長半径と短半径，中心座標を a, b, (0, 0)
円の半径と中心座標を r, (r, 0)
円と，楕円の接点を (x0，y0)
とおき，以下の連立方程式を解く。

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

using SymPy

@syms a::positive, b::positive, x0::positive, y0::positive, r::positive
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*x0/(a^2*y0) + (x0 - r)/y0
eq3 = (x0 - r)^2 + y0^2 - r^2
res = solve([eq1, eq2, eq3], (x0, y0, r))

   2-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (a*b/sqrt(a^2 - b^2), -b*sqrt((a^2 - 2*b^2)/(a - b))/sqrt(a + b), b*sqrt(a^2 - b^2)/a)
    (a*b/sqrt(a^2 - b^2), b*sqrt((a^2 - 2*b^2)/(a - b))/sqrt(a + b), b*sqrt(a^2 - b^2)/a)

2 組の解が得られるが，符号の違いのみなのでどちらでもよい。

接点間の距離 diag は以下のようになり，点対称な接点間の距離は，楕円の短径のみで決まる。常に，楕円の短径の √2 倍である。
すなわち，問のように短径が 1 寸ならば，接点間の距離は 1.41421356... である。

@syms d
diag = 2sqrt(res[2][1]^2 + res[2][2]^2)
diag |> println

   2*sqrt(a^2*b^2/(a^2 - b^2) + b^2*((a^2 - 2*b^2)/(a - b))/(a + b))

長い式は，以下のように簡約化される。

apart(diag, d) |> simplify |> println

   2*sqrt(2)*b

b に特定の値を代入して数値解を求める。

apart(diag, d)(b => 1/2) |> simplify |> println

   1.41421356237310

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

   b = 1.234; x0 = 1.35383;  y0 = 1.1012;  r = 1.12477

function draw(more)
    pyplot(size=(500, 500), grid=false, showaxis=true, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (3, 1.234)
   (x0, y0, r) = (a*b/sqrt(a^2 - b^2), b*sqrt((a^2 - 2*b^2)/(a - b))/sqrt(a + b), b*sqrt(a^2 - b^2)/a)
   println("点対称な接点間の距離 = $(2sqrt(x0^2 + y0^2)) = $(2√2b)")
   @printf("b = %g; x0 = %g;  y0 = %g;  r = %g\n", b, x0, y0, r)
   plot()
   ellipse(0, 0, a, b, color=:red)
   circle2(r, 0, r, :blue)
   segment(x0, y0, -x0, -y0, :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, b, " b", :red, :left, :bottom, delta=delta/2)
       point(x0, y0, "(x0,y0)", :red, :left, :bottom, delta=delta)
       point(-x0, -y0, "(-x0,-y0)", :red, :left, :bottom, delta=delta)
       #plot!(xlims=(-0.1,0.1), ylims=(-0.1,0.1))
   end
end;

 

算額（その816）

算額（その816）

文化三年丙寅正月 藤田貞資門人 石州津和野 二介泉尹
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf

直角三角形の田がある。直角を挟む二辺の短い方（鈎），長い方（股）が 6 間，8 間のとき，鈎，股に平行な直線で区切ってできる長方形（長辺（長），短辺（平）で，頂点の一つは斜辺上にある）の田の面積が最大になるとき，長，平はいかほどか。

鈎，股を A, B，鈎，股上にある点を a（平）, b（長） とおくと，長方形の面積 S は S = a*b である。
以下の連立方程式を解く。

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

using SymPy

@syms A, B, a, b, S

eq1 = (A - a)/b - A/B
eq2 = a/(B - b) - A/B
eq3 = a*b - S
res = solve([eq1, eq2, eq3], (a, b, S))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (A*(B - b)/B, b, A*b*(B - b)/B)

S = res[1][3]
S |> println

   A*b*(B - b)/B

問のように，A = 6, B = 8 のとき，b の取り方で，面積 S は b の関数で，S(b) は以下のように変化する。

pyplot(size=(300, 150), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(S(A => 6, B => 8), xlims=(0, 8), xlabel="b", ylabel="S(b)")

面積が最大になるの上図の曲線の接線の傾きが 0 になるときである。数値的に解くには，S(b) を b で微分し，S'(b) = 0 になるときの b の値を求める。

diff(S, b) |> println

   -A*b/B + A*(B - b)/B

solve(diff(S, b), b)[1] |> println

   B/2

b = B/2，a = A/2 が解である。
すなわち，長方形の短辺が 3 間，長辺が 6 間である。

function draw(more)
    pyplot(size=(300, 300), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (A, B) = (6, 8)
   (a, b) = (3, 4)
   plot([0, B, 0, 0], [0, 0, A, 0])
   rect(0, 0, b, a, :red)
   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)
       plot!(showaxis=false)
       point(0, A, " A（鈎）", :blue, :left, :bottom, delta=delta/2)
       point(B, 0, " B（股）", :blue, :left, :bottom, delta=delta/2)
       point(0, a, " a（平）", :red, :left, :bottom, delta=delta/2)
       point(b, 0, " b（長）", :red, :left, :bottom, delta=delta/2)
   end
end;