SymPy Plotting Backends で楕円や回転楕円体を描く

理工系の数学 B の授業で，楕円の周長や面積，回転楕円体の表面積や体積を求めているので。

必要なモジュールの import

In [1]:
from sympy.abc import *
from sympy import *
# SymPy Plotting Backends (SPB) 
from spb import *

# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']

mathtext.fontset の設定

In [2]:
# グラフを描くためではなくデフォルト設定のため
import matplotlib.pyplot as plt

plt.rcParams['mathtext.fontset'] = 'cm'

塗りつぶした楕円

SymPy Plotting Backends で長半径 $a$，単半径 $b$ の楕円を描く例。Ellipse() を使ってみる。

In [3]:
a = 2
b = 1
alim = a+0.5
blim = b+0.5

p = plot_geometry(
    Ellipse(Point(0, 0), hradius=a, vradius=b),
    {"facecolor": "yellow", "edgecolor":"darkblue", "lw":3}, 
    aspect = 'equal', size=(6.4, 4), n=400, show=False);

fig, ax = p.fig, p.ax
fig.tight_layout()

ax.grid(linestyle = 'dotted')
# 表示範囲
ax.set_xlim(-alim, alim)
ax.set_ylim(-blim, blim)
# 目盛設定
ax.set_xticks([-a, 0, a], ["$-a$", "$0$", "$a$"], fontsize=16)
ax.set_yticks([-b, 0, b], ["$-b$", "$0$", "$b$"], fontsize=16)

# x 軸，y 軸
ax.axhline(0, color='lightgray', dashes=(3, 3), linewidth=0.6)
ax.axvline(0, color='lightgray', dashes=(3, 3), linewidth=0.6);

塗りつぶし無しの楕円

In [4]:
p = plot_geometry(
    Ellipse(Point(0, 0), hradius=a, vradius=b),
    {"color":"darkblue", "lw":3}, is_filled=False,
    aspect = 'equal', size=(6.4, 4), n=400, show=False);

fig, ax = p.fig, p.ax
fig.tight_layout()

ax.grid(linestyle = 'dotted')
# 表示範囲
ax.set_xlim(-alim, alim)
ax.set_ylim(-blim, blim)
# 目盛設定
ax.set_xticks([-a, 0, a], ["$-a$", "$0$", "$a$"], fontsize=16)
ax.set_yticks([-b, 0, b], ["$-b$", "$0$", "$b$"], fontsize=16)

# x 軸，y 軸
ax.axhline(0, color='lightgray', dashes=(3, 3), linewidth=0.6)
ax.axvline(0, color='lightgray', dashes=(3, 3), linewidth=0.6);

長軸の周りに回転した回転楕円体

\begin{eqnarray}
\frac{x^2}{a^2} + \frac{y^2}{b^2} &=& 1, \quad b \equiv a \sqrt{1-e^2}\\
\therefore\ \ y &=& b \sqrt{1-\frac{x^2}{a^2}} \equiv f(x)
\end{eqnarray}

$y = f(x)$ を $x$ の周りに回転してできる回転楕円体の表面を描く。まず，

In [5]:
def f(x):
    return b*(sqrt(1-x**2/a**2))

$-a \leq x \leq a$ の範囲で $0 \leq y \leq f(x)$ の領域を塗りつぶして表示。

In [6]:
# 塗りつぶし部分
# plot_implicit は不等式で範囲指定できるが，
# 少し時間がかかる
p1 = plot_implicit(
    (0 <= y) & (y <= f(x)), (x, -a, a), (y, -b, b), 
    adaptive=True, size=(6.4, 4), 
    aspect = 'equal',
    # 塗りつぶしの色
    color='lightblue', show = False
)

p2 = plot(
    # y = f(x)
    (f(x), (x, -a, a), {"color":"darkblue", "lw":3}), 
    # 長軸
    (0, (x, -alim, alim), {"color":"blue", "lw":2}),
    show = False)

p = (p1+p2)
fig, ax = p.fig, p.ax
fig.tight_layout()

# グリッドを点線で
ax.grid(linestyle = 'dotted')
# 表示範囲
ax.set_xlim(-alim, alim)
ax.set_ylim(-blim, blim)
# 目盛設定
ax.set_xticks([-a, 0, a], ["$-a$", "$0$", "$a$"], fontsize=16)
ax.set_yticks([-b, 0, b], ["$-b$", "$0$", "$b$"], fontsize=16)

ax.set_xlabel("")
ax.set_ylabel("")

# x 軸，y 軸
ax.axhline(0, color='lightgray', dashes=(3, 3), linewidth=0.6)
ax.axvline(0, color='lightgray', dashes=(3, 3), linewidth=0.6);
# 凡例を非表示に
ax.get_legend().set_visible(False);

$y = f(x)$ を $x$ 軸の周りに $\theta$ だけ回転させると，

\begin{eqnarray}
x &=& x \\
y &=& f(x) \cos\theta \\
z &=& f(x) \sin\theta
\end{eqnarray}

となることに留意。

In [7]:
# 回転楕円体の表面
p1 = plot3d_parametric_surface(
    x, f(x)*cos(t), f(x)*sin(t), (x, -a, a), (t, 0, 2*pi),
    {'cmap':'Blues', 'lw':0.5, 'edgecolor':'lightblue'}, 
    use_cm=True, color_func=f(x)*cos(t), colorbar=False, 
    axis=False, size=(6.4, 6.4), n=40, 
    xlim=(-alim,alim), ylim=(-blim, blim), zlim=(-blim, blim), 
    aspect="equal", show=False
);

p2 = plot3d_parametric_line(
    # x 軸
    (x, 0, 0, (x, -alim, a), "", {"color":"blue","zorder":0, "lw":2}), 
    (x, 0, 0, (x, a,  alim), "", {"color":"blue","zorder":10,"lw":2}), 
    # z = z(x)
    (x, 0, f(x), (x, -a, -a+0.1), "", {"color":"darkblue","zorder":0, "lw":3}),
    (x, 0, f(x), (x, -a+0.1,  a), "", {"color":"darkblue","zorder":10,"lw":3}),
    use_cm=False, show=False
)

(p1+p2).show();

短軸の周りに回転した回転楕円体

$\displaystyle x = a \sqrt{1-\frac{y^2}{b^2}} \equiv g(y)$ を $y$ の周りに回転してできる回転楕円体の表面。

In [8]:
def g(y):
    return a*(sqrt(1-y**2/b**2))

$-b \leq y \leq b$ の範囲で $0 \leq x \leq g(y)$ の領域を塗りつぶして表示。

In [9]:
# 塗りつぶし部分
p1 = plot_implicit(
    (0 <= x) & (x <= g(y)), (x, 0, a), (y, -b, b), 
    adaptive=True, size=(6.4, 4), 
    aspect = 'equal',
    # 塗りつぶしの色
    color='lightpink', show = False
)

p2 = plot_parametric(
    # x = g(y)
    (g(y), y, (y, -b, b), {"color":"darkred", "lw":3}),
    # 短軸
    (0, y, (y, -blim, blim), {"color":"red", "lw":2}),
    use_cm=False, show = False)

p = (p1+p2)
fig, ax = p.fig, p.ax
fig.tight_layout()

ax.grid(linestyle = 'dotted')

# 表示範囲
ax.set_xlim(-alim, alim)
ax.set_ylim(-blim, blim)
# 目盛設定
ax.set_xticks([-a, 0, a], ["$-a$", "$0$", "$a$"], fontsize=16)
ax.set_yticks([-b, 0, b], ["$-b$", "$0$", "$b$"], fontsize=16)

ax.set_xlabel("")
ax.set_ylabel("")

# x 軸，y 軸
ax.axhline(0, color='lightgray', dashes=(3, 3), linewidth=0.6)
ax.axvline(0, color='lightgray', dashes=(3, 3), linewidth=0.6);
# 凡例を非表示に
ax.get_legend().set_visible(False);

$x = g(z)$ を $z$ 軸の周りに $\theta$ だけ回転させると，

\begin{eqnarray}
x &=& g(z) \cos\theta \\
y &=& g(z) \sin\theta \\
z &=& z
\end{eqnarray}

となることに留意する。

In [10]:
# 回転楕円体の表面
p1 = plot3d_parametric_surface(
    g(z)*cos(t), g(z)*sin(t), z, (z, -b, b), (t, 0, 2*pi),
    {'cmap':'Reds', 'lw':0.5, 'edgecolor':'pink'}, 
    use_cm=True, color_func=g(z)*sin(t), colorbar=False, 
    axis=False, size=(6.4, 6.4), n=80, 
    xlim=(-2.2,2.2), ylim=(-2.2, 2.2), zlim=(-1.2,1.2), 
    aspect="equal", show=False
    )

p2 = plot3d_parametric_line(
    # z 軸
    (0, 0, z, (z, -2, b), "", {"color":"red","zorder":0, "lw":2}), 
    (0, 0, z, (z, b,  2), "", {"color":"red","zorder":10,"lw":2}), 
    # x = g(z)
    (g(z), 0, z, (z, -2, -0.5), "", {"color":"darkred","zorder":0,"lw":3}), 
    (g(z), 0, z, (z, -0.5,  b), "", {"color":"darkred","zorder":10,"lw":3}), 
    use_cm=False, show=False
    )

(p1+p2).show();