算額（その439）

算額（その439）

東都本郷真光寺天満宮 天明7年丁未11月(1787)
東都湯嶋御徒町 池田重次郎教政

藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

外円の一部分である円弧内に，大中小の 5 個の円が入っている。大円と小円は弦に接している。
弦の長さが 8 寸，小円の直径が 1 寸のとき，矢（弦と円の距離，図では r0 - a）の長さはいかほどか。

外円の半径と中心座標を r0, (0, 0)
大円の半径と中心座標を r1, (x1, a + r1)
中円の半径と中心座標を r2, (0, r0 - r2)
小円の半径と中心座標を r3, (r3, a + r3)
弦と y 軸の交点座標を (0, a)
として以下の連立方程式の解を求める。

include("julia-source.txt");

using SymPy
@syms r0::positive,
     r1::positive, x1::positive,
     r2::positive, r3::positive,
     弦::positive, a::positive;

eq1 = x1^2 + (a + r1)^2 - (r0 - r1)^2
eq2 = (x1 - r3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = x1^2 + (r0 - r2 - a - r1)^2 - (r1 + r2)^2
eq4 = r3^2 + (r0 - r2 - a - r3)^2 - (r2 + r3)^2 
eq5 = (r0^2 - (弦/2)^2) - a^2;

res = solve([eq1, eq2, eq3, eq4, eq5], (r0, r1, x1, r2, a))

   2-element Vector{NTuple{5, Sym}}:
    ((1536*r3^6 + 288*r3^4*弦^2 - 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 84*r3^2*弦^4 + 52*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 - 弦^5*sqrt(16*r3^2 + 5*弦^2))/(32*r3*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)), r3*(128*r3^4 - 8*r3^2*弦^2 - 16*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 3*弦^4 + 弦^3*sqrt(16*r3^2 + 5*弦^2))/(2*(256*r3^4 + 32*r3^2*弦^2 + 弦^4)), r3*弦*(2*弦 + sqrt(16*r3^2 + 5*弦^2))/(16*r3^2 + 弦^2), r3*(96*r3^4 + 44*r3^2*弦^2 + 20*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 7*弦^4 + 3*弦^3*sqrt(16*r3^2 + 5*弦^2))/(8*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)), (-1536*r3^6 + 1440*r3^4*弦^2 + 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 180*r3^2*弦^4 + 20*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 - 弦^5*sqrt(16*r3^2 + 5*弦^2))/(4608*r3^5 - 1280*r3^3*弦^2 + 32*r3*弦^4))
    ((1536*r3^6 + 288*r3^4*弦^2 + 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 84*r3^2*弦^4 - 52*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 + 弦^5*sqrt(16*r3^2 + 5*弦^2))/(32*r3*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)), r3*(128*r3^4 - 8*r3^2*弦^2 + 16*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 3*弦^4 - 弦^3*sqrt(16*r3^2 + 5*弦^2))/(2*(256*r3^4 + 32*r3^2*弦^2 + 弦^4)), r3*弦*(2*弦 - sqrt(16*r3^2 + 5*弦^2))/(16*r3^2 + 弦^2), r3*(96*r3^4 + 44*r3^2*弦^2 - 20*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 7*弦^4 - 3*弦^3*sqrt(16*r3^2 + 5*弦^2))/(8*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)), (-1536*r3^6 + 1440*r3^4*弦^2 - 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 180*r3^2*弦^4 - 20*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 + 弦^5*sqrt(16*r3^2 + 5*弦^2))/(4608*r3^5 - 1280*r3^3*弦^2 + 32*r3*弦^4))

最初の組が適解である。

矢の長さは 2*r3*(3*r3^2 - 4*sqrt(r3^2 + 20) - 24)/(9*r3^2 - 16)

res[1][1] - res[1][5] |> simplify |> println

   r3*(24*r3^2 - 3*弦^2 - 弦*sqrt(16*r3^2 + 5*弦^2))/(36*r3^2 - 弦^2)

弦の長さが 8 寸，小円の直径が 1 寸 のとき，矢の長さは 3 寸である。

(r3, 弦) = (1/2, 8)
r3*(24*r3^2 - 3*弦^2 - 弦*sqrt(16*r3^2 + 5*弦^2))/(36*r3^2 - 弦^2)

   3.0

   r0 = 4.16667;  r1 = 1.125;  x1 = 2;  r2 = 1.04167;  a = 1.16667
   矢 = 3

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r3, 弦) = (1/2, 8)
   (r0, r1, x1, r2, a) = ((1536*r3^6 + 288*r3^4*弦^2 - 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 84*r3^2*弦^4 + 52*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 - 弦^5*sqrt(16*r3^2 + 5*弦^2))/(32*r3*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)),
       r3*(128*r3^4 - 8*r3^2*弦^2 - 16*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 3*弦^4 + 弦^3*sqrt(16*r3^2 + 5*弦^2))/(2*(256*r3^4 + 32*r3^2*弦^2 + 弦^4)),
       r3*弦*(2*弦 + sqrt(16*r3^2 + 5*弦^2))/(16*r3^2 + 弦^2),
       r3*(96*r3^4 + 44*r3^2*弦^2 + 20*r3^2*弦*sqrt(16*r3^2 + 5*弦^2) + 7*弦^4 + 3*弦^3*sqrt(16*r3^2 + 5*弦^2))/(8*(144*r3^4 - 40*r3^2*弦^2 + 弦^4)),
       (-1536*r3^6 + 1440*r3^4*弦^2 + 64*r3^4*弦*sqrt(16*r3^2 + 5*弦^2) - 180*r3^2*弦^4 + 20*r3^2*弦^3*sqrt(16*r3^2 + 5*弦^2) + 3*弦^6 - 弦^5*sqrt(16*r3^2 + 5*弦^2))/(4608*r3^5 - 1280*r3^3*弦^2 + 32*r3*弦^4))
   @printf("r0 = %g;  r1 = %g;  x1 = %g;  r2 = %g;  a = %g\n",
           r0, r1, x1, r2, a)
   @printf("矢 = %g\n", r0 - a)
   plot()
   circle(0, 0, r0, :gray)
   circle(x1, a + r1, r1)
   circle(-x1, a + r1, r1)
   circle(0, r0 - r2, r2, :blue)
   circle(r3, a + r3, r3, :orange)
   circle(-r3, a + r3, r3, :orange)
   b = sqrt(r0^2 - a^2)
   segment(-b, a, b, a)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, a + r1, "大円:r1,(x1,a+r1)", :red, :center, :bottom, delta=delta/2)
       point(0, r0 - r2, "中円:r2,(0,r0-r2)", :blue, :center, :bottom, delta=delta/2)
       point(r3, a + r3, " 小円:r3,(r3,a+r3)", :black, :left, :vcenter)
       point(0, a, " a", :black, :left, :top, delta=-delta/2)
       point(弦/2, a, "(弦/2,a)", :black, :right, :top, delta=-delta/2)
       point(0, r0, " r0", :black, :left, :bottom, delta=delta/2)
       point(r0, 0, "r0 ", :black, :right, :bottom, delta=delta/2)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5, ylims=(-0.3, 4.3))
   else
       plot!(showaxis=false)
   end
end;