算額（その790）
寛政八年丙辰十一月 丸山因平良玄門人 参州苅屋 林政右衛門盛保
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf
十字線を隔てて乾円，坤円，巽円，艮円の4円と，中央に容円を置く。艮円，坤円，巽円の直径がそれぞれ 15 寸，10 寸，6 寸のとき，乾円の直径を求めよ。
容円の半径と中心座標を r0, (x0, y0)
坤円の半径と中心座標を r1, (r1, r1)
乾円の半径と中心座標を r2, (r2, -r2)
巽円の半径と中心座標を r3, (-r3, r3)
艮円の半径と中心座標を r4, (-r4, -r4)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r0::positive, x0::negative, y0::positive,
r1::positive, r2::positive, r3::positive, r4::positive
eq1 = (r1 - x0)^2 + (r1 - y0)^2 - (r0 + r1)^2
eq2 = (r2 - x0)^2 + (-r2 - y0)^2 - (r0 + r2)^2
eq3 = (-r3 - x0)^2 + (r3 - y0)^2 - (r0 + r3)^2
eq4 = (-r4 - x0)^2 + (-r4 - y0)^2 - (r0 + r4)^2
res = solve([eq1, eq2, eq3, eq4], (r2, r0, x0, y0));
4 番目の組が適解であるが，SymPy では簡約化できない長い式になる。
res[4][1] |> println
sqrt(8*r1^2*r3*r4^2*(r1*r3 - r1*r4 + r3^2 + r3*r4)*(r1^2*r3 + r1^2*r4 + r1*r3^2 + 6*r1*r3*r4 - r1*r3*sqrt(r1^2 + 6*r1*r4 + r4^2) + r1*r4^2 - r1*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) + r3^2*r4 - r3^2*sqrt(r1^2 + 6*r1*r4 + r4^2) + r3*r4^2 - r3*r4*sqrt(r1^2 + 6*r1*r4 + r4^2)) + (r1^3*r3^2 - 2*r1^3*r3*r4 + r1^3*r4^2 + r1^2*r3^3 + 6*r1^2*r3^2*r4 - r1^2*r3^2*sqrt(r1^2 + 6*r1*r4 + r4^2) - 10*r1^2*r3*r4^2 + r1^2*r4^3 + r1^2*r4^2*sqrt(r1^2 + 6*r1*r4 + r4^2) + 4*r1*r3^3*r4 - r1*r3^3*sqrt(r1^2 + 6*r1*r4 + r4^2) + 6*r1*r3^2*r4^2 - r1*r3^2*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) - 2*r1*r3*r4^3 + r3^3*r4^2 - r3^3*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) + r3^2*r4^3 - r3^2*r4^2*sqrt(r1^2 + 6*r1*r4 + r4^2))^2)/(4*r1*r4*(r1*r3 - r1*r4 + r3^2 + r3*r4)) + (-r1^3*r3^2 + 2*r1^3*r3*r4 - r1^3*r4^2 - r1^2*r3^3 - 6*r1^2*r3^2*r4 + r1^2*r3^2*sqrt(r1^2 + 6*r1*r4 + r4^2) + 10*r1^2*r3*r4^2 - r1^2*r4^3 - r1^2*r4^2*sqrt(r1^2 + 6*r1*r4 + r4^2) - 4*r1*r3^3*r4 + r1*r3^3*sqrt(r1^2 + 6*r1*r4 + r4^2) - 6*r1*r3^2*r4^2 + r1*r3^2*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) + 2*r1*r3*r4^3 - r3^3*r4^2 + r3^3*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) - r3^2*r4^3 + r3^2*r4^2*sqrt(r1^2 + 6*r1*r4 + r4^2))/(4*r1*r4*(r1*r3 - r1*r4 + r3^2 + r3*r4))
res[4][2] |> println
(-r1^3*r3^2 + r1^3*r3*r4 - 2*r1^3*r4^2 - r1^2*r3^3 - 5*r1^2*r3^2*r4 + r1^2*r3^2*sqrt(r1^2 + 6*r1*r4 + r4^2) + 2*r1^2*r3*r4^2 + r1^2*r3*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) - 2*r1^2*r4^3 - 2*r1*r3^3*r4 + r1*r3^3*sqrt(r1^2 + 6*r1*r4 + r4^2) - 5*r1*r3^2*r4^2 + 2*r1*r3^2*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) + r1*r3*r4^3 + r1*r3*r4^2*sqrt(r1^2 + 6*r1*r4 + r4^2) - r3^3*r4^2 + r3^3*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) - r3^2*r4^3 + r3^2*r4^2*sqrt(r1^2 + 6*r1*r4 + r4^2))/(8*r1*r4*(r1*r3 - r1*r4 + r3^2 + r3*r4))
res[4][3] |> println
(3*r1^3*r3 - r1^3*r4 + 9*r1^2*r3*r4 - r1^2*r3*sqrt(r1^2 + 6*r1*r4 + r4^2) + r1^2*r4^2 - r1^2*r4*sqrt(r1^2 + 6*r1*r4 + r4^2) - 3*r1*r3^3 - 7*r1*r3^2*r4 - r3^3*r4 + r3^3*sqrt(r1^2 + 6*r1*r4 + r4^2) - r3^2*r4^2 + r3^2*r4*sqrt(r1^2 + 6*r1*r4 + r4^2))/(8*r1*(r1*r3 - r1*r4 + r3^2 + r3*r4))
res[4][4] |> println
res[4][4](r1 => 5, r3 => 3, r4 => 7.5).evalf() |> println
-(r1 + r3)*(r3 - r4)*(r3 + r4)*sqrt(r1^2 + 6*r1*r4 + r4^2)/(8*r4*(r1*r3 - r1*r4 + r3^2 + r3*r4)) + (r3 - r4)*(r1^2*r3 + r1^2*r4 + r1*r3^2 + 8*r1*r3*r4 - r1*r4^2 + 3*r3^2*r4 + 3*r3*r4^2)/(8*r4*(r1*r3 - r1*r4 + r3^2 + r3*r4))
res[4][1](r1 => 5, r3 => 3, r4 => 7.5).evalf() |> println
res[4][2](r1 => 5, r3 => 3, r4 => 7.5).evalf() |> println
res[4][3](r1 => 5, r3 => 3, r4 => 7.5).evalf() |> println
res[4][4](r1 => 5, r3 => 3, r4 => 7.5).evalf() |> println
10.0000000000000
2.70833333333333
0.666666666666667
-1.37500000000000
まえもって r1, r3, r4 に定数を代入しておいてから解くと数値解が求まる。
@syms r0::positive, x0::negative, y0::positive,
r1::positive, r2::positive, r3::positive, r4::positive
@syms r0, x0, y0, r1, r2, r3, r4
(r1, r3, r4) = (10, 6, 15) .// 2
eq1 = (r1 - x0)^2 + (r1 - y0)^2 - (r0 + r1)^2
eq2 = (r2 - x0)^2 + (-r2 - y0)^2 - (r0 + r2)^2
eq3 = (-r3 - x0)^2 + (r3 - y0)^2 - (r0 + r3)^2
eq4 = (-r4 - x0)^2 + (-r4 - y0)^2 - (r0 + r4)^2
res = solve([eq1, eq2, eq3, eq4], (r2, r0, x0, y0))
4-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
(-18, -305/8, 17, -207/8)
(55/2, -305/8, 17, -207/8)
(-1/2, 65/24, 2/3, -11/8)
(10, 65/24, 2/3, -11/8)
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r3, r4) = (10, 6, 15) .// 2
(r2, r0, x0, y0) = [10, 65/24, 2/3, -11/8]
plot()
circle(r1, r1, r1)
circle(r2, -r2, r2, :blue)
circle(-r3, r3, r3, :green)
circle(-r4, -r4, r4, :magenta)
circle(x0, y0, r0, :orange)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:black, lw=0.5)
vline!([0], color=:black, lw=0.5)
point(x0, y0, "容円:r0,(x0,y0)", :black, :left, delta=-delta/2)
point(r1, r1, "坤円:r1,(r1,r1)", :red, :center, delta=-delta/2)
point(r2, -r2, "乾円:r2,(r2,-r2)", :blue, :center, delta=-delta/2)
point(-r3, r3, "巽円:r3\n(-r3,r3)", :green, :center, :bottom, delta=delta/2)
point(-r4, -r4, "艮円:r4\n(-r4,-r4)", :magenta, :center, :bottom, delta=delta/2)
end
end;
