長野県上水内牟礼村牟礼 渋薬師堂 嘉永2年(1849)
長野県上水内郡牟飯綱町 岩崎観音堂 文化15年(1818)
長野県下高井郡山ノ内町 渋温泉医王殿 慶應元年(1865)
長野県上水内郡飯綱町 牟礼神社 明治31年(1898)
中村信弥「改訂増補 長野県の算額」県内の算額(P.99, 179, 216, 261)
http://www.wasan.jp/zoho/zoho.html
キーワード：3次元，球，回転楕円体
#Julia #SymPy #算額 #和算 #数学
回転楕円体の中に，甲球 2 個，乙球 3 個を容れる。回転楕円体の長径と短径が与えられたとき，甲球の直径を得る術を述べよ。
回転楕円体の長径，短径を \(2a,\ 2b\)
甲球の半径と中心座標を \(r_1,\ (0,\ 0,\ z_1)\)
乙球の半径と中心座標を \(r_2,\ (b - r_2,\ 0,\ 0)\)
甲球と回転楕円体の \(x-z\) 平面上の交点座標を \( (x_0,\ 0,\ z_0)\)
とおき，以下の連立方程式を解く。
まず，eq1 を解き \(r_2\) を求める。
次いで， eq3, eq4, eq5 を解き，\(z_1,\ x_0,\ z_0\) を求める。
最後に，eq2 を解き \(r_1\) を求める。
include("julia-source.txt"); # julia-source.txt ソース
using SymPy
@syms a::positive, b::positive, x0::positive,
z0::positive, r1::positive, z1::positive,
r2::positive
eq1 = (2r2/√Sym(3) + r2) - b
eq2 = z1^2 + (b - r2)^2 - (r1 + r2)^2 |> expand
eq3 = x0^2/a^2 + z0^2/b^2 - 1 |> expand
eq4 = (x0 - z1)^2 + z0^2 - r1^2 |> expand
eq5 = -b^2*x0/(a^2*z0) + (x0 - z1)/z0 |> expand;
# r2 乙球の半径
ans_r2 = solve(eq1, r2)[1] |> simplify
ans_r2 |> println
b*(-3 + 2*sqrt(3))
# z1, x0, z0
res = solve([eq3, eq4, eq5], (z1, x0, z0))[4] # 4 of 4
(sqrt( (a*b^2 - a*r1^2 + b^3 - b*r1^2)/(a - b))*(a - b)/b, a^2*sqrt( (b^2 - r1^2)/(a - b))/(b*sqrt(a + b)), sqrt( (a^2*r1^2 - b^4)/(a - b))/sqrt(a + b))
eq2 に，ここまでで得られた \(r_2,\ x_1,\ x_0,\ z_0\) を代入し，\(r_1\) を求める。
eq12 = eq2(r2 => ans_r2, z1 => res[1], x0 => res[2], z0 => res[3]) |> simplify;
# r1 甲球の半径
ans_r1 = solve(eq12, r1)[1]
ans_r1 |> println
b*(a^2 - 4*sqrt(3)*b^2 + 6*b^2)/a^2
甲球の半径は \(b(a^2 - 4\sqrt{3}b^2 + 6b^2)/a^2\) である。
長径が 4，短径が 2 のとき，甲球の直径は 1.53589838486225 である。
2*ans_r1(a => 4/2, b => 2/2).evalf()|> println
1.53589838486225
術は以下のようになっている。
@syms 短径, 長径
甲球径 = 2(短径/2 - (2√3 -3)*短径^3/長径^2)
甲球径(短径 => 2, 長径 => 4) |> println
1.53589838486225
すべてのパラメータは以下の通りである。
回転楕円体の長径が 4，短径が 2 のとき，甲球の直径は 1.5359，乙球の直径は 0.928203 である。
\(a = 2; b = 1; r_1 = 0.767949; r_2 = 0.464102\)
\(z_1 = 1.1094; x_0 = 1.4792; z_0 = 0.673049\)
描画関数プログラムのソースを見る
function draw1(a, b, r1, z1, r2, x0, z0, more)
p1 = plot()
circle(0, 0, b, :green)
rotate(0, b - r2, r2, :blue)
circle(0, 0, r1)
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(0, b, "b", :green, :center, :bottom, delta=delta/2)
point(0, b - r2, "乙球:r2(0,b-r2,0)", :blue, :center, delta=-delta/2)
point(0, r1, "r1", :red, :center, :bottom, delta=delta/2)
end
end;
function draw2(a, b, r1, z1, r2, x0, z0, more)
plot()
circle2(z1, 0, r1)
circle(0, b - r2, r2, :blue)
circle(0, -(b - r2)/2, r2, :blue)
ellipse(0, 0, a, b, color=: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(0, b - r2, "乙球:r2\n(0,b-r2,0)", :blue, :center, delta=-delta)
point(a, 0, "a", :green, :left, :bottom, delta=delta)
point(0, b, "b", :green, :center, :bottom, delta=delta)
point(z1, 0, "甲球:r1,(0, 0, z1)", :red, :center, delta=-delta)
point(x0, z0, "(x0,z0)", :green, :left, :bottom, delta=delta)
end
end;
function draw(a, b, func, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = b*(a^2 - 4*sqrt(3)*b^2 + 6*b^2)/a^2
r2 = b*(-3 + 2*sqrt(3))
z1 = sqrt( (a*b^2 - a*r1^2 + b^3 - b*r1^2)/(a - b))*(a - b)/b
x0 = a^2*sqrt( (b^2 - r1^2)/(a - b))/(b*sqrt(a + b))
z0 = sqrt( (a^2*r1^2 - b^4)/(a - b))/sqrt(a + b)
@printf("回転楕円体の長径が %g，短径が %g のとき，甲球の直径は %g，乙球の直径は %g である。\n", 2a, 2b, 2r1, 2r2)
@printf("a = %g; b = %g; r1 = %g; r2 = %g; z1 = %g; x0 = %g; z0 = %g\n", a, b, r1, r2, z1, x0, z0)
func(a, b, r1, z1, r2, x0, z0, more)
end;
draw(4/2, 2/2, draw1, true)
draw(4/2, 2/2, draw2, true)
以下のアイコンをクリックして応援してください
