Matplotlib で楕円や回転楕円体を描く

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

必要なモジュールの import

In [1]:
import matplotlib.pyplot as plt
import numpy as np

from matplotlib import patches

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

mathtext.fontset の設定

In [2]:
plt.rcParams['mathtext.fontset'] = 'cm'

塗りつぶした楕円

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

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

fig = plt.figure(figsize=[6.4, 4])
ax = fig.add_subplot(aspect='equal')
fig.tight_layout()

e1 = patches.Ellipse((0, 0), 2*a, 2*b, 
                     # 塗りつぶしの色，  線の色
                     facecolor="yellow", edgecolor="darkblue",
                     linewidth=2, fill=True)
ax.add_patch(e1)

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

# グリッド（格子線）は dotted に
ax.grid(ls = ':')
# x軸 y軸は dashed に
ax.axhline(0, c='k', ls = '--', lw=1)
ax.axvline(0, c='k', ls = '--', lw=1);

塗りつぶし無しの楕円

In [4]:
fig = plt.figure(figsize=[6.4,4])
ax = fig.add_subplot(aspect='equal')
fig.tight_layout()

e1 = patches.Ellipse((0, 0), 2*a, 2*b, 
                     edgecolor="darkblue",
                     linewidth=2, fill=False, zorder=0)
ax.add_patch(e1)

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

# グリッド（格子線）は dotted に
ax.grid(ls = ':')
# x軸 y軸は dashed に
ax.axhline(0, c='k', ls = '--', lw=1)
ax.axvline(0, c='k', ls = '--', lw=1);

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

\begin{eqnarray}
\frac{x^2}{a^2} + \frac{y^2}{b^2} &=& 1
\end{eqnarray}

を $y$ について解いて

\begin{eqnarray}
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*(np.sqrt(1-x**2/a**2))
In [6]:
fig = plt.figure(figsize=[6.4, 4])
ax = fig.add_subplot(aspect='equal')
fig.tight_layout()

# 塗りつぶし部分
x = np.linspace(-a, a, 300)
ax.fill_between(x, 0, f(x), facecolor = "lightblue")

# y = f(x)
ax.plot(x, f(x), lw = 3, c = "darkblue")

# 長軸
ax.axhline(0, c = 'blue', lw = 2)

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

# グリッド（格子線）は dotted に
ax.grid(ls = ':')
# x軸 y軸は dashed に
# ax.axhline(0, c='k', ls = '--', lw=1)
ax.axvline(0, c='k', ls = '--', lw=1);

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

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

となることに留意してデータを作成する。

In [7]:
fig = plt.figure(figsize=[6.4, 6.4])
fig.tight_layout()
ax = fig.add_subplot(projection='3d')
ax.set_box_aspect((alim, blim, blim))

# 回転楕円体の表面
Ndiv = 50
xi = np.linspace(-a, a, Ndiv)
th = np.linspace(0, 2*np.pi, Ndiv)
# np.outer でメッシュデータに
x = np.outer(xi, np.ones(len(xi)))
y = np.outer(f(xi), np.cos(th))
z = np.outer(f(xi), np.sin(th))
#
ax.plot_surface(x, y, z, 
                color="lightblue", 
                edgecolor="blue", lw=0.2)

# 長軸
X = np.linspace(-alim, alim, 100)
ax.plot(X, [0]*100, [0]*100, lw=2, c="blue")

# 奥の z = f(x)
xdiv = 30
x = np.linspace(-a+0.1, -a, xdiv)
ax.plot(x, [0]*xdiv, f(x), lw=3, c="darkblue", zorder=0)
# 手前の z = f(x)
xdiv = 300
x = np.linspace(-a+0.1, a, xdiv)
ax.plot(x, [0]*xdiv, f(x), lw=3, c="darkblue", zorder=10)


ax.set_xlim(-alim, alim)
ax.set_ylim(-blim, blim)
ax.set_zlim(-blim, blim)
ax.axis(False);

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

\begin{eqnarray}
\frac{x^2}{a^2} + \frac{y^2}{b^2} &=& 1
\end{eqnarray}

を $x$ について解いて

\begin{eqnarray}
x &=& a \sqrt{1-\frac{y^2}{b^2}} \equiv g(y)
\end{eqnarray}

$x = g(y) $ を $y$ の周りに回転してできる回転楕円体の表面。

In [8]:
def g(y):
    return a*(np.sqrt(1-y**2/b**2))
In [9]:
fig = plt.figure(figsize=[6.4, 4])
ax = fig.add_subplot(aspect='equal')
fig.tight_layout()

# 塗りつぶし部分
x = np.linspace(0, a, 200)
ax.fill_between(x, f(x), -f(x), facecolor="lightpink")

# x = g(y)
y = np.linspace(-b, b, 200)
ax.plot(g(y), y, linewidth = 3, color="darkred")

# 短軸
ax.axvline(0, color='red', linewidth=2)

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

# グリッド（格子線）は dotted に
ax.grid(ls = ':')
# x軸 y軸は dashed に
ax.axhline(0, c='k', ls = '--', lw=1);
# ax.axvline(0, c='k', ls = '--', lw=1);

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

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

となることに留意してデータを作成する。

In [10]:
fig = plt.figure(figsize=[6.4, 6.4])
ax = fig.add_subplot(projection='3d')
ax.set_box_aspect((alim, alim, blim))

# 回転楕円体の表面
Ndiv = 50
zi = np.linspace(-b, b, Ndiv)
th = np.linspace(0, 2*np.pi, Ndiv)
# np.outer でメッシュデータに
x = np.outer(g(zi), np.cos(th))
y = np.outer(g(zi), np.sin(th))
z = np.outer(zi, np.ones(len(zi)))
# 
ax.plot_surface(x, y, z, 
                color="lightpink",
                edgecolor="red", linewidth=0.2)

# 短軸
Z = np.linspace(b, 2, 30)
ax.plot([0]*len(Z), [0]*len(Z), Z, lw=2, c="red", zorder=10)
Z = np.linspace(-2, b, 30)
ax.plot([0]*len(Z), [0]*len(Z), Z, lw=2, c="red", zorder=0)

# 手前の x = g(z)
zdiv = 100
z = np.linspace(-0.3, b, zdiv)
ax.plot(g(z), [0]*zdiv, z, lw=3, c="darkred", zorder=10)
# 奥の x = g(z)
z = np.linspace(-b, -0.3, zdiv)
ax.plot(g(z), [0]*zdiv, z, lw=3, c="darkred", zorder=0)

ax.set_xlim(-alim, alim)
ax.set_ylim(-alim, alim)
ax.set_zlim(-blim, blim)
ax.axis(False);