長野県佐久市広川原 禅昌寺 文政11年(1828)
中村信弥「改訂増補 長野県の算額」県内の算額(P.108)
http://www.wasan.jp/zoho/zoho.html
キーワード：3次元，球，減弧円錐
#Julia #SymPy #算額 #和算 #数学
減弧円錐の中に甲球 3 個，乙球 1 個，丙球 3 個を容れる。甲球は互いに外接し，底面と側面に接し，乙球とも外接する。丙球は乙球と外接し，側面とも接している。乙球は甲球，丙球と外接し，側面にも接している。
丙球の直径が甲球の直径の 1/5(2 分)のとき，乙球の直径の最小値を求めよ。
減弧円錐とは，左側の図において，赤で描いた 2 個の円を z 軸を中心として回転させたときにできる先の尖った円錐状の立体である。
大球の半径と中心座標を \(R, (R, 0, R)\)
甲球の半径と中心座標を\( r_1, (2r_1/\sqrt{3}, 0, r_1)\)
乙球の半径と中心座標を \(r_2, (0, 0, z_2)\)
丙球の半径と中心座標を \(r_3, (2r_3/\sqrt{3}, 0, z_3)\)
とおき，以下の連立方程式を解く。
\(r_1, r_2, z_2, r_3, z_3\) は \(R\) の関数として得られる。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms R, r1, r2, z2, r3, z3
eq1 = (R -2r1/√Sym(3))^2 + (R - r1)^2 - (R + r1)^2
eq2 = R^2 + (R - z2)^2 - (R + r2)^2
eq3 = (2r1/√Sym(3))^2 + (z2 - r1)^2 - (r1 + r2)^2
eq4 = (R - 2r3/√Sym(3))^2 + (R - z3)^2 - (R + r3)^2
eq5 = (2r3/√Sym(3))^2 + (z3 - z2)^2 - (r2 + r3)^2;
ans_r1 = solve(eq1, r1)[2] # 2 of 2
ans_r1 |> println
R*(-sqrt(9 + 6*sqrt(3)) + sqrt(3) + 3)/2
eq12 = eq2(r1 => ans_r1)
eq13 = eq3(r1 => ans_r1);
res = solve([eq12, eq13], (r2, z2))[1]
(R*(-1940*sqrt(3 + 2*sqrt(3)) - 1120*sqrt(9 + 6*sqrt(3)) + 2848*sqrt(3) + 4933)/(4*(-277*sqrt(3)*sqrt(3 + 2*sqrt(3)) - 478*sqrt(3 + 2*sqrt(3)) + 1218 + 704*sqrt(3))), R*(-427*sqrt(9 + 6*sqrt(3)) - 736*sqrt(3 + 2*sqrt(3)) + 1081*sqrt(3) + 1877)/(4*(-39*sqrt(3)*sqrt(3 + 2*sqrt(3)) - 67*sqrt(3 + 2*sqrt(3)) + 98*sqrt(3) + 171)))
# r2
@syms d
ans_r2 = apart(res[1]/R, d) |> simplify |> sympy.sqrtdenest |> simplify |> x -> x*R
ans_r2 |> println
R*(-37*sqrt(2)*3^(3/4)/552 - 21*sqrt(2)*3^(1/4)/184 + 13/46 + 4*sqrt(3)/23)
# z2
ans_z2 = apart(res[2]/R, d) |> simplify |> sympy.sqrtdenest |> simplify |> x -> x*R
ans_z2 |> println
R*(-67*sqrt(2)*3^(1/4)/184 - 83*sqrt(2)*3^(3/4)/552 + 4*sqrt(3)/23 + 59/46)
eq15 = eq5(r2 => ans_r2, z2 => ans_z2) |> simplify
eq14 = eq4 |> expand
res2 = solve([eq14, eq15], (r3, z3))[2] # 2 of 2
( (-192*sqrt(3)*R - 312*R + 120*sqrt(2)*3^(3/4)*R + 264*sqrt(2)*3^(1/4)*R + (3*R*sqrt(-133695296*sqrt(2)*3^(3/4) - 231566928*sqrt(2)*3^(1/4) + 497992069 + 287516218*sqrt(3))/(2*(2630*sqrt(2)*3^(3/4) + 4590*sqrt(2)*3^(1/4) + 17642*sqrt(3) + 30675)) + 3*R*(-5537*sqrt(2)*3^(3/4) - 9587*sqrt(2)*3^(1/4) + 16095*sqrt(3) + 27957)/(2*(2630*sqrt(2)*3^(3/4) + 4590*sqrt(2)*3^(1/4) + 17642*sqrt(3) + 30675)))*(-201*sqrt(2)*3^(1/4) - 83*sqrt(2)*3^(3/4) + 156 + 96*sqrt(3)))/(63*sqrt(2)*3^(1/4) + 37*sqrt(2)*3^(3/4) + 396 + 272*sqrt(3)), 3*R*sqrt(-133695296*sqrt(2)*3^(3/4) - 231566928*sqrt(2)*3^(1/4) + 497992069 + 287516218*sqrt(3))/(2*(2630*sqrt(2)*3^(3/4) + 4590*sqrt(2)*3^(1/4) + 17642*sqrt(3) + 30675)) + 3*R*(-5537*sqrt(2)*3^(3/4) - 9587*sqrt(2)*3^(1/4) + 16095*sqrt(3) + 27957)/(2*(2630*sqrt(2)*3^(3/4) + 4590*sqrt(2)*3^(1/4) + 17642*sqrt(3) + 30675)))
# r3
@syms d
ans_r3 = apart(res2[1]/R, d) |> simplify |> sympy.sqrtdenest |> simplify |> x -> x*R
ans_r3 |> println
R*(-1457609*sqrt(2)*3^(3/4) - 2524500*sqrt(2)*3^(1/4) + 3^(1/4)*(-201 - 83*sqrt(3))*sqrt(-267390592*sqrt(2)*3^(3/4) - 463133856*sqrt(2)*3^(1/4) + 995984138 + 575032436*sqrt(3))/2 + 78*sqrt(-133695296*sqrt(2)*3^(3/4) - 231566928*sqrt(2)*3^(1/4) + 497992069 + 287516218*sqrt(3)) + 48*sqrt(-401085888*sqrt(2)*3^(3/4) - 694700784*sqrt(2)*3^(1/4) + 1493976207 + 862548654*sqrt(3)) + 6137076 + 3544118*sqrt(3))/(2618169*sqrt(2)*3^(1/4) + 1512127*sqrt(2)*3^(3/4) + 9518764 + 5497344*sqrt(3))
# z3
ans_z3 = apart(res2[2]/R, d) |> simplify |> sympy.sqrtdenest |> simplify |> x -> x*R
ans_z3 |> println
3*R*(-5537*sqrt(2)*3^(3/4) - 9587*sqrt(2)*3^(1/4) + sqrt(-133695296*sqrt(2)*3^(3/4) - 231566928*sqrt(2)*3^(1/4) + 497992069 + 287516218*sqrt(3)) + 16095*sqrt(3) + 27957)/(2*(2630*sqrt(2)*3^(3/4) + 4590*sqrt(2)*3^(1/4) + 17642*sqrt(3) + 30675))
\(R = 1\) のとき，\(r_1 = 0.164191; r_2 = 0.155332; z_2 = 0.421387; r_3 = 0.0357628; z_3 = 0.607967\) となる。
\(r_3/r_1 = 0.2178121821537112\) である。
「答」は，乙球径 = 1.32678 である。なぜこのような答えになるのか？？
描画関数プログラムのソースを見る
function draw(R)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = R*(-sqrt(9 + 6*sqrt(3)) + sqrt(3) + 3)/2
r2 = R*(-37*sqrt(2)*3^(3/4)/552 - 21*sqrt(2)*3^(1/4)/184 + 13/46 + 4*sqrt(3)/23)
z2 = R*(-67*sqrt(2)*3^(1/4)/184 - 83*sqrt(2)*3^(3/4)/552 + 4*sqrt(3)/23 + 59/46)
r3 = R*(-1457609*sqrt(2)*3^(3/4) - 2524500*sqrt(2)*3^(1/4) + 3^(1/4)*(-201 - 83*sqrt(3))*sqrt(-267390592*sqrt(2)*3^(3/4) - 463133856*sqrt(2)*3^(1/4) + 995984138 + 575032436*sqrt(3))/2 + 78*sqrt(-133695296*sqrt(2)*3^(3/4) - 231566928*sqrt(2)*3^(1/4) + 497992069 + 287516218*sqrt(3)) + 48*sqrt(-401085888*sqrt(2)*3^(3/4) - 694700784*sqrt(2)*3^(1/4) + 1493976207 + 862548654*sqrt(3)) + 6137076 + 3544118*sqrt(3))/(2618169*sqrt(2)*3^(1/4) + 1512127*sqrt(2)*3^(3/4) + 9518764 + 5497344*sqrt(3))
z3 = 3*R*(-5537*sqrt(2)*3^(3/4) - 9587*sqrt(2)*3^(1/4) + sqrt(-133695296*sqrt(2)*3^(3/4) - 231566928*sqrt(2)*3^(1/4) + 497992069 + 287516218*sqrt(3)) + 16095*sqrt(3) + 27957)/(2*(2630*sqrt(2)*3^(3/4) + 4590*sqrt(2)*3^(1/4) + 17642*sqrt(3) + 30675))
@printf("R = %g; r1 = %g; r2 = %g; z2 = %g; r3 = %g; z3 = %g\n", R, r1, r2, z2, r3, z3)
p1 = plot(xlabel="x-axis", ylabel="z-axis")
circle2(R, R, R)
circle(2r1/√3, r1, r1, :blue)
circle(-r1/√3, r1, r1, :blue)
circle(0, z2, r2, :green)
circle(2r3/√3, z3, r3, :magenta)
circle(-r3/√3, z3, r3, :magenta)
point(0.8, 0.9, "甲円:r2,(2r1/√3,0,r1)", :blue, :right, :bottom, delta=0.1R, mark=false)
point(0.8, 1.1, "乙円:r2,(0,0,z2)", :green, :right, :bottom, delta=0.1R, mark=false)
point(0.8, 1.3, "丙円:r3,(2r3/√3,0,z3)", :magenta, :right, :bottom, delta=0.08R, mark=false)
point(R, R, "大円:R,(R,R)", :red, :right, delta=-0.1R)
plot!(xlims=(-1, 1), ylims = (0, 1))
p2 = plot(xlabel="x-axis", ylabel="y-axis")
circle2(R, 0, R)
rotate(2r1/√3, 0, r1, :blue)
circle(0, 0, r2, :green)
rotate(2r3/√3, 0, r3, :magenta)
plot!(xlims=(-0.5, 0.5), ylims = (-0.5, 0.5))
plot(p1, p2)
end;
draw(1)
以下のアイコンをクリックして応援してください
