算額（その1213）

算額（その1213）

(20) 兵庫県青垣町遠坂字後岶 熊野神社 明治18年(1885)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：長方形，直角三角形

長方形の中に曲尺（かねじゃく）を置く。長方形の長辺，短辺が 40 寸，38 寸，短曲（曲尺の短い方）が 10 寸のとき，曲（曲尺の長い方）と長方形の長辺を二辺とする直角三角形の第三の辺（鈎）の長さはいかほどか。

注：曲尺の短曲と曲とのなす角は直角である。

算額（その1198）と基本的に同じ問題である。
https://blog.goo.ne.jp/r-de-r/e/0e430beba5823ce19bc2449c0ff9b448

パラメータはそのまま方程式の変数名とする。
左上にある直角三角形の直角を挟む 2 辺の長さを x，38 -「鈎」，とおき，以下の連立方程式を立てる。
たかが三元連立方程式であるが，「長」，「平」，「短曲」を変数のままにして一度に「鈎」，および「曲」と x を求めることは(SymPy で自動的には）できない。

include("julia-source.txt");

using SymPy

@syms 長::positive, 平::positive, 曲::positive, 短曲::positive, 鈎::positive, x::positive
eq1 = (鈎*長 + x*(平 - 鈎) + 短曲*曲 + 平*(長 - x))/2 - 平*長
eq2 = (平 - 鈎)^2 + x^2 - 短曲^2
eq3 = x/鈎 - (平 - 鈎)/長;
#res = solve([eq1, eq2, eq3], (鈎, 曲, x))

eq3 を x について解き，解を ans_x とおく。

ans_x = solve(eq3, x)[1]
ans_x |> println

   鈎*(平 - 鈎)/長

eq2 の x に ans_x を代入し（式を簡約化して）eq12 とおく。

eq12 = eq2(x => ans_x) |> factor |> numerator
eq12 |> println

   平^2*鈎^2 + 平^2*長^2 - 2*平*鈎^3 - 2*平*鈎*長^2 - 短曲^2*長^2 + 鈎^4 + 鈎^2*長^2

eq12 を「鈎」について解くと，非常に複雑な長い式が得られる。SymPy では簡約化できないが，「長」，「平」，「短曲」を特定すれば，その場合の「鈎」が得られる。

ans_鈎 = solve(eq12, 鈎)[3];  # 3 of 4

「長」，「平」，「短曲」がそれぞれ 40 寸, 38 寸, 10 寸のとき，「鈎」は 30 寸である。

ans_鈎(長 => 40, 平 => 38, 短曲 => 10).evalf() |> println

   30.0000000000000

「鈎」が決まれば，すべての図形要素が逐次的に決定される。

ans_x(鈎 => 30, 長 => 40, 平 => 38) |> println

   6

solve(eq1(長 => 40, 平 => 38, 短曲 => 10, x => 6, 鈎 => 30), 曲)[1] |> println

   50

function draw(more=false)
   pyplot(size=(300, 300), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (長, 平, 短曲) = (40, 38, 10)
   (鈎, x) = (30, 6)
   plot([0, 平, 平, 0, 0], [0, 0, 長, 長, 0], color=:blue, lw=0.5)
   segment(0, 長 -x, 平 - 鈎, 長, :red)
   segment(平 - 鈎, 長, 平, 0, :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(平, 長, " (平,長)", :red, :left, :vcenter)
       dimension_line(平 - 鈎, 長 + 0.5, 平, 長 + 0.5, "鈎", :blue, :center, :bottom, delta=delta)
       dimension_line(-0.5, 長, -0.5, 長 - x, "x ", :blue, :right, :vcenter)
       point((平 - 鈎)/2, 長 - x/2, " 短曲", :red, mark=false)
       point(平 - 鈎/2, 長/2, "  曲", :green, mark=false)
       xlims!(-6delta, 平 + 20delta)
   end
end;

draw(true)

算額（その1212）

算額（その1212）

(20) 兵庫県青垣町遠坂字後岶 熊野神社 明治18年(1885)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円8個

互いに外接し合う，木円 2 個，火円 2 個，土円 2 個，金円 1 個，水円 1 個がある。
木円，火円の直径が与えられたとき，金円の直径を求める術を述べよ。

木円の半径と中心座標を r1, (r1, y1)
火円の半径と中心座標を r2, (r2, y2)
土円の半径と中心座標を r3, (r3, 0)
金円の半径と中心座標を r4, (0, y4)
水円の半径と中心座標を r5, (0, y5)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, y1::negative, r2::positive, y2::positive, r3::positive, r4::positive, y4::negative, r5::positive, y5::positive

eq1 = (r1 - r2)^2 + (y2 - y1)^2 - (r1 + r2)^2
eq2 = (r1 - r3)^2 + y1^2 - (r1 + r3)^2
eq3 = (r2 - r3)^2 + y2^2 - (r2 + r3)^2;

res1 = solve([eq1, eq2, eq3], (y1, y2, r3))[1]

   (-2*sqrt(r1)*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2))/sqrt(r2), 2*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2)), (-2*r1^(3/2)*r2^(3/2)*(r1 - r2)^2 + r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2))/((r1 - r2)^2*(r1^2 - 2*r1*r2 + r2^2)))

eq4 = r1^2 + (y1 - y4)^2 - (r1 + r4)^2
eq5 = r3^2 + y4^2 - (r3 + r4)^2
eq6 = r2^2 + (y2 - y5)^2 - (r2 + r5)^2
eq7 = r3^2 + y5^2 - (r3 + r5)^2;
#res2 = solve([eq4, eq5, eq6, eq7], (r4, y4, r5, y5))

次いで，r4, y4, r5, y5 を求めるが，前段で求めた y1, y2, r3 は数式が複雑なので代入せず，変数のまま解く。

res2 = solve([eq4, eq5, eq6, eq7], (r4, y4, r5, y5))[4]

   (y1*(y1*(r1 - r3)*sqrt(4*r1*r3 + y1^2)/(-r1^2 + 2*r1*r3 - r3^2 + y1^2) + y1 - y1*(2*r1*r3 - 2*r3^2 + y1^2)/((-r1 + r3 + y1)*(r1 - r3 + y1)))/(2*(r1 - r3)), -y1*(r1 - r3)*sqrt(4*r1*r3 + y1^2)/(2*(-r1^2 + 2*r1*r3 - r3^2 + y1^2)) + y1*(2*r1*r3 - 2*r3^2 + y1^2)/(2*(-r1 + r3 + y1)*(r1 - r3 + y1)), y2*(y2*(r2 - r3)*sqrt(4*r2*r3 + y2^2)/(-r2^2 + 2*r2*r3 - r3^2 + y2^2) + y2 - y2*(2*r2*r3 - 2*r3^2 + y2^2)/((-r2 + r3 + y2)*(r2 - r3 + y2)))/(2*(r2 - r3)), -y2*(r2 - r3)*sqrt(4*r2*r3 + y2^2)/(2*(-r2^2 + 2*r2*r3 - r3^2 + y2^2)) + y2*(2*r2*r3 - 2*r3^2 + y2^2)/(2*(-r2 + r3 + y2)*(r2 - r3 + y2)))

res2[1] |> simplify |> println

   y1^2*(-r1 - r3 + sqrt(4*r1*r3 + y1^2))/(2*(-r1^2 + 2*r1*r3 - r3^2 + y1^2))

res2[2] |> simplify |> println

   y1*(2*r1*r3 - r1*sqrt(4*r1*r3 + y1^2) - 2*r3^2 + r3*sqrt(4*r1*r3 + y1^2) + y1^2)/(2*(-r1^2 + 2*r1*r3 - r3^2 + y1^2))

res2[3] |> simplify |> println

   y2^2*(-r2 - r3 + sqrt(4*r2*r3 + y2^2))/(2*(-r2^2 + 2*r2*r3 - r3^2 + y2^2))

res2[4] |> simplify |> println

   y2*(2*r2*r3 - r2*sqrt(4*r2*r3 + y2^2) - 2*r3^2 + r3*sqrt(4*r2*r3 + y2^2) + y2^2)/(2*(-r2^2 + 2*r2*r3 - r3^2 + y2^2))

金円の直径の式に y1, y2, r3 を代入すればよいが，SymPy では簡約化できないとてつもなく長い式になる。

y1 = -2*sqrt(r1)*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2))/sqrt(r2)
y2 = 2*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2))
r3 = (-2*r1^(3/2)*r2^(3/2)*(r1 - r2)^2 + r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2))/((r1 - r2)^2*(r1^2 - 2*r1*r2 + r2^2))
金径 = y1^2*(-r1 - r3 + sqrt(4*r1*r3 + y1^2))/(2*(-r1^2 + 2*r1*r3 - r3^2 + y1^2))
金径 |> println

   4*r1*(-r1 + sqrt(4*r1*(r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2) - 2*r1^1.5*r2^1.5*(r1 - r2)^2)/((r1 - r2)^2*(r1^2 - 2*r1*r2 + r2^2)) + 4*r1*((r1^2*r2^2 + r1*r2^3 - 2*r1^1.5*r2^2.5)/(r1^2 - 2*r1*r2 + r2^2))/r2) - (r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2) - 2*r1^1.5*r2^1.5*(r1 - r2)^2)/((r1 - r2)^2*(r1^2 - 2*r1*r2 + r2^2)))*(r1^2*r2^2 + r1*r2^3 - 2*r1^1.5*r2^2.5)/(r2*(r1^2 - 2*r1*r2 + r2^2)*(-2*r1^2 + 4*r1*(r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2) - 2*r1^1.5*r2^1.5*(r1 - r2)^2)/((r1 - r2)^2*(r1^2 - 2*r1*r2 + r2^2)) + 8*r1*((r1^2*r2^2 + r1*r2^3 - 2*r1^1.5*r2^2.5)/(r1^2 - 2*r1*r2 + r2^2))/r2 - 2*(r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2) - 2*r1^1.5*r2^1.5*(r1 - r2)^2)^2/((r1 - r2)^4*(r1^2 - 2*r1*r2 + r2^2)^2)))

r1，r2 に具体例を代入すれば，正しい金円の直径が求まるので，式は間違いない。

金径(r1 => 8/2, r2 => 10/2) |> println

   0.804254212447020

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

   r1 = 4;  y1 = 4.22291;  r2 = 5;  y2 = -4.72136;  r3 = 1.11456;  r4 = 0.804254;  y4 = -1.64362;  r5 = 0.871325;  y5 = -1.64362

function draw(r1, r2, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   y1 = -2*sqrt(r1)*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2))/sqrt(r2)
   y2 = 2*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2))
   r3 = (-2*r1^(3/2)*r2^(3/2)*(r1 - r2)^2 + r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2))/((r1 - r2)^2*(r1^2 - 2*r1*r2 + r2^2))
   (y1, y2, y3) =  (-2*sqrt(r1)*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2))/sqrt(r2), 2*sqrt((-2*r1^(3/2)*r2^(5/2) + r1^2*r2^2 + r1*r2^3)/(r1^2 - 2*r1*r2 + r2^2)), (-2*r1^(3/2)*r2^(3/2)*(r1 - r2)^2 + r1*r2*(r1 + r2)*(r1^2 - 2*r1*r2 + r2^2))/((r1 - r2)^2*(r1^2 - 2*r1*r2 + r2^2)))
   r4 = y1^2*(-r1 - r3 + sqrt(4*r1*r3 + y1^2))/(2*(-r1^2 + 2*r1*r3 - r3^2 + y1^2))
   y4 = y1*(2*r1*r3 - r1*sqrt(4*r1*r3 + y1^2) - 2*r3^2 + r3*sqrt(4*r1*r3 + y1^2) + y1^2)/(2*(-r1^2 + 2*r1*r3 - r3^2 + y1^2))
   r5 = y2^2*(-r2 - r3 + sqrt(4*r2*r3 + y2^2))/(2*(-r2^2 + 2*r2*r3 - r3^2 + y2^2))
   y5 = y2*(2*r2*r3 - r2*sqrt(4*r2*r3 + y2^2) - 2*r3^2 + r3*sqrt(4*r2*r3 + y2^2) + y2^2)/(2*(-r2^2 + 2*r2*r3 - r3^2 + y2^2))
   @printf("r1 = %g;  y1 = %g;  r2 = %g;  y2 = %g;  r3 = %g;  r4 = %g;  y4 = %g;  r5 = %g;  y5 = %g\n", r1, y1, r2, y2, r3, r4, y5, r5, y5)
   plot()
   circle2(r1, y1, r1)
   circle2(r2, y2, r2, :blue)
   circle2(r3, 0, r3, :green)
   circle(0, y4, r4, :magenta)
   circle(0, y5, r5, :orange)
   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, y1, "木円:r1,(r1,y1)", :red, :center, delta=-delta/2)
       point(r2, y2, "火円:r2,(r2,y2)", :blue, :center, delta=-delta/2)
       point(r3, 0, "土円:r3,(r3,0)", :green, :left, delta=-delta/2)
       point(0, y4, " 金円:r4,(0,r4)", :magenta, :left, :vcenter)
       point(0, y5, " 水円:r5,(0,r5)", :orange, :left, :vcenter)
   end
end;

draw(10/2, 8/2, true)