算額（その622）

算額（その622）

和算図形問題あれこれ
令和 2 年 12 月の問題 - No.2?
飯塚文庫『関流叚書見題部円類五十問本書』より
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html

正方形内に斜線を引き，甲円，乙円，丙円を入れる。正方形の一辺の長さが 1 寸のとき，甲円の直径はいかほどか。

正方形の一辺の長さを a
甲円の半径と中心座標を r1, (r1, r1)
乙円の半径と中心座標を r2, (a - r2, a - r2)
丙円の半径と中心座標を r3, (r3, a - r2)
斜線の y 切片の座標を (0, b)
とおき，以下の連立方程式を解く。

乙円と丙円の半径と中心座標の間の関係を表す 2 つの式は同等の条件を表現しているので，どちらを使っても方程式は解ける。しかし，今回の場合は方程式1 を使うと SymPy は 'NotImplementedError' を起こしてしまい，解が求まらない。方程式2 を使えば解ける。

方程式1: (a - r2 - r3)^2 + (r2 - r3)^2 = (r2 + r3)^2
方程式2: a - r2 - r3 = sqrt((r2 + r3)^2 - (r2 - r3)^2)

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

using SymPy

@syms r1::positive, r2::positive, r3::positive, a::positive, b::positive
@syms r1, r2, r3, a, b

a = 1
eq1 = (a - r2 - r3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq2 = a - r2 - r3  - sqrt((r2 + r3)^2 - (r2 - r3)^2)
eq3 = a + b - sqrt(a^2 + b^2) - 2r1
eq4 = distance(a, 0, 0, b, a - r2, a - r2) - r2^2
eq5 = distance(a, 0, 0, b, r3, a - r3) - r3^2;res = solve([eq2, eq3, eq4, eq5], (r1, r2, r3, b));

2 組の解が得られるが，2 番目のものが適解である。
それぞれの式はかなり長い。甲円の半径は以下のような式になる。

res[2][1] |> println

   sqrt(7)*(-14^(3/4)*(3*sqrt(183) + 62)^(1/3)*sqrt(3*(-sqrt(14)*(3*sqrt(183) + 62)^(1/6)*sqrt(-44*sqrt(14)*(3*sqrt(183) + 62)^(1/3)*sqrt(-22*(3*sqrt(183) + 62)^(1/3) + 91 + 7*(3*sqrt(183) + 62)^(2/3)) - 7*sqrt(14)*(3*sqrt(183) + 62)^(2/3)*sqrt(-22*(3*sqrt(183) + 62)^(1/3) + 91 + 7*(3*sqrt(183) + 62)^(2/3)) - 91*sqrt(14)*sqrt(-22*(3*sqrt(183) + 62)^(1/3) + 91 + 7*(3*sqrt(183) + 62)^(2/3)) + 444*sqrt(21)*sqrt(3*sqrt(183) + 62)) + 2*2^(1/4)*sqrt(3)*7^(3/4)*(3*sqrt(183) + 62)^(1/3)*(-22*(3*sqrt(183) + 62)^(1/3) + 91 + 7*(3*sqrt(183) + 62)^(2/3))^(1/4) + 14^(3/4)*(3*sqrt(183) + 62)^(1/6)*(-22*(3*sqrt(183) + 62)^(1/3) + 91 + 7*(3*sqrt(183) + 62)^(2/3))^(3/4))^2 + 3087*sqrt(14)*(3*sqrt(183) + 62)^(2/3)*sqrt(-22*(3*sqrt(183) + 62)^(1/3) + 91 + 7*(3*sqrt(183) + 62)^(2/3)))*(-22*(3*sqrt(183) + 62)^(1/3) + 91 + 7*(3*sqrt(183) + 62)^(2/3))^(5/4)/4116 - 14^(1/4)*sqrt(3)*sqrt(3*sqrt(183) + 62)*(-22*(3*sqrt(183) + 62)^(1/3) + 91 + 7*(3*sqrt(183) + 62)^(2/3))^(5/4)*sqrt(-44*sqrt(14)*(3*sqrt(183) + 62)^(1/3)*sqrt(-22*(3*sqrt(183) + 62)^(1/3) + 91 + 7*(3*sqrt(183) + 62)^(2/3)) - 7*sqrt(14)*(3*sqrt(183) + 62)^(2/3)*sqrt(-22*(3*sqrt(183) + 62)^(1/3) + 91 + 7*(3*sqrt(183) + 62)^(2/3)) - 91*sqrt(14)*sqrt(-22*(3*sqrt(183) + 62)^(1/3) + 91 + 7*(3*sqrt(183) + 62)^(2/3)) + 444*sqrt(21)*sqrt(3*sqrt(183) + 62))/294 - sqrt(42)*(3*sqrt(183) + 62)^(2/3)*(-7*sqrt(3*sqrt(183) + 62) + 22*(3*sqrt(183) + 62)^(1/6))*(-22*(3*sqrt(183) + 62)^(1/3) + 91 + 7*(3*sqrt(183) + 62)^(2/3))/294 + sqrt(126*sqrt(183) + 2604)*(-28028*(3*sqrt(183) + 62)^(1/3) + 115934 + 8918*(3*sqrt(183) + 62)^(2/3))/4116 + 9*sqrt(7)*(3*sqrt(183) + 62)^(2/3)*(-22*(3*sqrt(183) + 62)^(1/3) + 91 + 7*(3*sqrt(183) + 62)^(2/3))^(3/2)/98)/((3*sqrt(183) + 62)^(2/3)*(-22*(3*sqrt(183) + 62)^(1/3) + 91 + 7*(3*sqrt(183) + 62)^(2/3))^(3/2))

それぞれの式を評価した結果を以下に示す。
順に，甲円の半径，乙円の半径，丙円の半径，斜線の y 切片 b

res[2][1].evalf(), res[2][2].evalf(), res[2][3].evalf(), res[2][4].evalf()

   (0.246773416321210, 0.336189121793897, 0.176552761905338, 0.734031532621976)

甲円の直径は 0.493546832642421 寸である。

2res[2][1].evalf() |> println

   0.493546832642421

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1
   (r1, r2, r3, b) = (0.246773416321210, 0.336189121793897, 0.176552761905338, 0.734031532621976)
   @printf("甲円の直径 = %g;  r1 = %g;  r2 = %g;  r3 = %g;  b = %g\n", 2r1, r1, r2, r3, b)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5, xlims=(-0.05, 1.05))
   circle(r1, r1, r1)
   circle(a - r2, a - r2, r2, :blue)
   circle(r3, a - r3, r3, :magenta)
   segment(0, b, a, 0, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       point(r1, r1, "甲円:r1,(r1,r1)", :red, :center, delta=-delta/2)
       point(a - r2, a - r2, "乙円:r2,(a-r2,a-r2)", :blue, :center, delta=-delta/2)
       point(r3, a - r3, "丙円:r3,(r3,a-r3)", :magenta, :center, delta=-delta/2)
       point(0, a, "a ", :green, :right, :vcenter)
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
       point(0, b, "b ", :green, :right, :vcenter)
   end
end;