からっぽのしょこ

読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

円周の作図【Matplotlib】

はじめに

 Matplotlibライブラリを利用して、円のグラフを作成します。

【目次】

円周の作図

 Matplotlibライブラリを利用して、2次元空間(平面)上に円周(circumference)のグラフを作成します。また、円座標系(circular coordinates)をグラフで確認します。

 利用するライブラリを読み込みます。

# 利用ライブラリ
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation


座標の計算

 まずは、円周上の点の座標計算を確認します。

 角度と半径を指定して、円周上の点(ノルムが指定した半径になる2次元ベクトル)を計算します。

# 半径を指定
r = 14.0

# 度数法の角度を指定
t = 180.0

# 弧度法の角度(ラジアン)を計算
t *= np.pi / 180.0
print(t)

# 円周上の点の座標を計算
x = r * np.cos(t)
y = r * np.sin(t)

# 座標を格納
p = np.array([x, y])
print(p)
print(np.linalg.norm(p))
3.141592653589793
[-1.40000000e+01  1.71450552e-15]
14.0

 度数法における角度を \theta^{\circ}、弧度法における角度(ラジアン)を \thetaで表すことにします。度数法と弧度法は \theta = \frac{2 \pi}{360^{\circ}} \theta^{\circ} \theta^{\circ} = \frac{360}{2 \pi} \thetaの関係です。つまり、 \theta^{\circ} = 0^{\circ}のとき \theta = 0 \theta^{\circ} = 180^{\circ}のとき \theta = \pi \theta^{\circ} = 360^{\circ}のとき \theta = 2 \piです。
  \piは円周率で、np.piで扱えます。

 半径を r、角度(x軸とのなす角)を \thetaとします。オブジェクト名はr, tとします。
 半径が rの円周上の点(ノルムが rの2次元ベクトル)を \mathbf{p} = (x, y)とすると、x軸・y軸の座標は、次の式で計算できます。

 \displaystyle
\begin{cases}
    x = r \cos \theta \\
    y = r \sin \theta
\end{cases}
\tag{1}

 直交座標系の座標を (x, y)、円座標系の座標を (r, \theta)で表します。
 式(1)は、円座標系から直交座標系への変換式です。

  \mathbf{p}のユークリッドノルム \|\mathbf{p}\| = \sqrt{x^2 + y^2} rになります。ユークリッドノルムは、np.linalg.norm()で計算できます。

 あるいは、2次元ベクトルを指定して、球面上の点(ノルムが指定した半径になるベクトル)に変換します。

# ベクトルを指定
p = np.array([10.0, 17.0])
#p = np.array([x, y]) # 確認用
print(np.linalg.norm(p))

# 半径を指定
r = 14.0

# ノルムをrに変換(円周上の座標を計算)
p *= r / np.linalg.norm(p)
print(p)
print(np.linalg.norm(p))
19.72308292331602
[ 7.09828177 12.06707901]
14.0

 任意の2次元ベクトル \tilde{\mathbf{p}} = (\tilde{x}, \tilde{y})を指定して、ノルムで割り正規化(ノルムを1に変換)します。さらに、半径 rを掛けるとノルムが rのベクトルになります。

 \displaystyle
\mathbf{p}
    = r
      \frac{
          \tilde{\mathbf{p}}
      }{
          \|\tilde{\mathbf{p}}\|
      }
    = r
      \frac{
          \tilde{\mathbf{p}}
      }{
          \sqrt{\tilde{x}^2 + \tilde{y}^2}
      }

 この計算により、2次元空間上の点 \tilde{\mathbf{p}}と角度 \thetaが同じである球面上の点(ベクトル \tilde{\mathbf{p}}と同じ方向であるノルムが rのベクトル) \mathbf{p}が得られます。ベクトルのスカラ倍なのでベクトルの向きは変わりません。

 円周上の点から半径と角度を計算します。

# 半径を計算
r = np.linalg.norm(p)
print(r)

# 弧度法の角度(ラジアン)を計算
sgn_y = 1.0 if p[1] >= 0.0 else -1.0
t = sgn_y * np.arccos(p[0] / np.linalg.norm(p))
print(t)

# 度数法の角度を計算
t *= 180.0 / np.pi
print(t)
14.0
1.0390722595360908
59.53445508054011

 半径 rと角度 \thetaは、次の式で計算できます。

 \displaystyle
\begin{align*}
&\begin{cases}
    r   = \|\mathbf{p}\|
        = \sqrt{x^2 + y^2} \\
    \theta
        = \mathrm{sgn}(y)
          \arccos \Bigl(
              \frac{z}{\|\mathbf{p}\|}
          \Bigr)
        = \mathrm{sgn}(y)
          \arccos \Bigl(
              \frac{x}{\sqrt{x^2 + y^2}}
          \Bigr) \\
\end{cases}
\tag{2}\\
&\mathrm{sgn}(y)
    = \begin{cases}
          1 & (y \geq 0) \\
          -1 & (y < 0)
      \end{cases}
\end{align*}

  \arccos xは逆コサイン関数(コサイン関数 \cos xの逆関数)で、np.arccos()で計算できます。y軸の値 yの符号を \mathrm{sgn}(y)で表しています。
 式(2)は、直交座標系から円座標系への変換式です。

 ここでは、2通りの座標計算方法を確認しました。以降は、1つ目の方法を用います。

円の描画

 次は、円周全体の座標を計算して、円のグラフを作成します。

 1周期分のラジアンを作成して、円周上の座標を計算します。周期については「変数と各軸の値の関係」で確認します。

# 半径を指定
r = 1.0

# ラジアンを作成
t = np.linspace(start=0.0, stop=2.0*np.pi, num=81)
print(t[:5])
print(len(t))

# 円周の座標を計算
x = r * np.cos(t)
y = r * np.sin(t)
print(x[:5].round(3))
print(y[:5].round(3))
[0.         0.07853982 0.15707963 0.23561945 0.31415927]
81
[1.    0.997 0.988 0.972 0.951]
[0.    0.078 0.156 0.233 0.309]

  0 \leq \theta \leq 2 \piの範囲のラジアンを作成してtとします。
  \thetaの範囲(最大値と最小値の差)が 2 \piであれば、他の値であっても同じグラフになります。
 tを用いて、座標の計算式(1)を計算して、x軸・y軸の値をそれぞれx, yとします。

 確認のため、点ごとにノルムを計算します。また、色付け用の配列を作成します。

# 2つの軸の値を結合
xy = np.stack([x, y], axis=1)
print(xy[:5])
print(xy.shape)
print(np.linalg.norm(xy, axis=1)[:5])

# 2つの軸の値と0ベクトルを結合
xyo = (np.stack([x, y, np.zeros_like(x)], axis=1) + r) * 0.5 / r
print(xyo[:5])
[[1.         0.        ]
 [0.99691733 0.0784591 ]
 [0.98768834 0.15643447]
 [0.97236992 0.23344536]
 [0.95105652 0.30901699]]
(81, 2)
[1. 1. 1. 1. 1.]
[[1.         0.5        0.5       ]
 [0.99845867 0.53922955 0.5       ]
 [0.99384417 0.57821723 0.5       ]
 [0.98618496 0.61672268 0.5       ]
 [0.97552826 0.6545085  0.5       ]]

 ノルムの計算用に、x, yをそれぞれ1列として結合した配列xyを作成します。
 np.linalg.norm()axix引数に1を指定して行(点)ごとにノルムを計算すると、どの点も値がrになるのを確認できます。

 散布図の装飾に関して、各点の座標 (x, y)を正規化した値をRGB値として使うことで、点ごとの座標に応じて色付けられます。
 x, yと全てが0の列を結合した配列xyoを作成します。どの列(軸)を0の列にするかによって、配色が変わります。
 各軸の値は最小値が -r、最大値が rです。そこで、各点に rを足すことで 0から 2 rの値になり、さらに 2 rで割ることで 0から 1の値に変換できます。x軸(1つ目の列)の値が大きい( x rに近く y 0に近い)ほど赤色、y軸(2つ目の列)の値が大きいほど緑色、3つ目の列の値が大きいほど青色に近い色になります。

 円のグラフを作成します。

# タイトル用の文字列(座標の計算式)を作成
coord_label = '$x = r\ \\cos \\theta$\n'
coord_label += '$y = r\ \\sin \\theta$'

# 円を作図
fig, ax = plt.subplots(figsize=(7, 7), facecolor='white')
ax.plot(x, y, color='black') # 曲線
#ax.scatter(x, y, c=y) # 散布図:(y軸の値により色付け)
#ax.scatter(x, y, color=xyo) # 散布図:(xy座標の値により色付け)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('$r=' + str(r) + '$\n' + coord_label, loc='left')
fig.suptitle('circle', fontsize=20)
ax.grid()
ax.set_aspect('equal')
plt.show()

円のグラフ

円のグラフ:散布図

 Axes.plot()で曲線(折れ線グラフ)、Axes.scatter()で散布図として円を描画できます。

 綺麗な円を描画するには、Axes.set_aspect('equal')を使ってアスペクト比を1に設定します。

 cmap引数にカラーマップ名を指定して、y軸の値などに応じてグラデーションで色付けられます。
 散布図に関して、color引数に「各点の座標 (x, y)を正規化した値をRGB値」として指定して、点ごとの座標に応じて色付けられます。

 円は、原点からの距離(ノルム)が rの点なのが分かります。

 以上で、円を描画できました。次からは、変数(角度)と2つの軸の値の関係を確認します。

変数と各軸の値の関係

 角度とx軸・y軸それぞれの値の関係をグラフで確認します。

 1周期分以上のラジアンを作成して、円周の座標を計算します。

# 半径を指定
r = 1.0

# ラジアンを作成
t = np.linspace(start=-2.0*np.pi, stop=2.0*np.pi, num=200)
print(t[:5])
print(len(t))

# 円周の座標を計算
x = r * np.cos(t)
y = r * np.sin(t)
print(x[:5].round(3))
print(y[:5].round(3))
[-6.28318531 -6.22003772 -6.15689013 -6.09374253 -6.03059494]
200
[1.    0.998 0.992 0.982 0.968]
[0.    0.063 0.126 0.188 0.25 ]

 「円の描画」のときと同様にして、円周の座標を計算します。ただしこの例では、周期性が分かるように、ラジアン \theta, \phiの範囲をそれぞれ -2 \piから 2 \piとします。

  \theta軸のラジアン目盛用の配列を作成します。グラフ自体には不要な処理です。

# 半円の目盛の数(分母の値)を指定
denom = 6.0

# 目盛の通し番号(分子の値)を作成
tick_vals = np.arange(start=-2.0*denom, stop=2.0*denom+0.1)
print(tick_vals[:5])

# 目盛用のラジアンを計算
tick_t_vals = tick_vals / denom * np.pi
print(tick_t_vals[:5].round(2))

# 目盛ラベルを作成
tick_label_vals = [
    '$' + ('-' if nomer < 0.0 else '') + '\\frac{' + str(abs(nomer)) + '}{' + str(int(denom)) + '} \pi$' for nomer in tick_vals.astype('int16')
]
print(tick_label_vals[:3])
[-12. -11. -10.  -9.  -8.]
[-6.28 -5.76 -5.24 -4.71 -4.19]
['$-\\frac{12}{6} \\pi$', '$-\\frac{11}{6} \\pi$', '$-\\frac{10}{6} \\pi$']

 軸目盛ラベルを \frac{i}{n} \piの形で表示することにします。この例では2周期分を描画するので、 i -2 nから 2 nの整数 i = -2 n, -2 (n-1), \dots , -1, 0, 1, \dots, 2 (n-1), 2 nとします。 ndenom itick_valsとして値を指定します。
 軸目盛の座標用のラジアン t = \frac{i}{n} \piを計算してtick_t_valsとします。

  \theta xの関係のグラフを作成します。

# x軸の値を曲線で作図
fig, ax = plt.subplots(figsize=(9, 3), facecolor='white')
ax.plot(t, x)
ax.set_xticks(ticks=tick_t_vals, labels=tick_label_vals) # θ目盛ラベル
ax.set_xlabel('$\\theta$')
ax.set_ylabel('$x$')
ax.set_title('$r=' + str(r) + '$', loc='left')
fig.suptitle('$x = r\ \\cos \\theta$', fontsize=20)
ax.grid()
ax.set_aspect('equal')
plt.show()

角度とx軸の値の関係:コサイン関数のグラフ

  \theta = t, t + 2 \pi xの値が一致するのが分かります。 2 \piごとに周期的に推移します。

  \theta yの関係のグラフを作成します。

# y軸の値を曲線で作図
fig, ax = plt.subplots(figsize=(9, 3), facecolor='white')
ax.plot(t, y)
ax.set_xticks(ticks=tick_t_vals, labels=tick_label_vals) # θ目盛ラベル
ax.set_xlabel('$\\theta$')
ax.set_ylabel('$y$')
ax.set_title('$r=' + str(r) + '$', loc='left')
fig.suptitle('$y = r\ \\sin \\theta$', fontsize=20)
ax.grid()
ax.set_aspect('equal')
plt.show()

角度とy軸の値の関係:サイン関数のグラフ

 こちらも 2 \piごとに周期します。

 ここでは、変数(角度)と2つの軸それぞれの値の関係を確認しました。次は、全ての軸の値との関係を確認します。

変数と座標の関係

 最後に、円周上に変数(角度)に関する軸目盛を描画して、変数と座標の関係をグラフで確認します。

 円周の座標を計算します。

# 半径を指定
r = 1.0

# ラジアンを作成
t = np.linspace(start=0.0, stop=2.0*np.pi, num=100)
print(t[:5])
print(len(t))

# 円周の座標を計算
x = r * np.cos(t)
y = r * np.sin(t)
[0.         0.06346652 0.12693304 0.19039955 0.25386607]
100

 「円の描画」のときのコードで、円周の座標を計算します。

  \theta軸目盛用の配列を作成します。

# 半円の目盛の数(分母の値)を指定
denom = 6.0

# 目盛の通し番号(分子の値)を作成
tick_vals = np.arange(start=0.0*denom, stop=2.0*denom)
print(tick_vals)

# 目盛用のラジアンを計算
tick_t_vals = tick_vals / denom * np.pi
print(tick_t_vals.round(2))

# 目盛ラベルを作成
tick_label_vals = [
    '$' + ('-' if nomer < 0.0 else '') + '\\frac{' + str(abs(nomer)) + '}{' + str(int(denom)) + '} \pi$' for nomer in tick_vals.astype('int16')
]
print(tick_label_vals[:3])

# 目盛用の座標を計算
tick_x_vals = r * np.cos(tick_t_vals)
tick_y_vals = r * np.sin(tick_t_vals)
print(tick_x_vals.round(2))
print(tick_y_vals.round(2))
[ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11.]
[0.   0.52 1.05 1.57 2.09 2.62 3.14 3.67 4.19 4.71 5.24 5.76]
['$\\frac{0}{6} \\pi$', '$\\frac{1}{6} \\pi$', '$\\frac{2}{6} \\pi$']
[ 1.    0.87  0.5   0.   -0.5  -0.87 -1.   -0.87 -0.5  -0.    0.5   0.87]
[ 0.    0.5   0.87  1.    0.87  0.5   0.   -0.5  -0.87 -1.   -0.87 -0.5 ]

 「変数と各軸の値の関係」のときと同様にして、ラジアン目盛用の値を作成します。この例では1周期分を描画するので、 i 0から 2 nの整数 i = 0, 1, \dots, 2 nとします。
 さらに、軸目盛用のラジアンtick_t_valsを用いて、目盛の描画位置の座標を計算してtick_*_valsとします。

 なす角 \thetaの角マークと角ラベル用の配列を作成します。

# 角マーク用のラジアンを作成
angle_t_vals = np.linspace(start=0.0, stop=1.0/denom*np.pi, num=50)

# 角マーク用の座標を計算
d = 0.1
angle_x_vals = d * np.cos(angle_t_vals)
angle_y_vals = d * np.sin(angle_t_vals)

# 角ラベル用の座標を計算
d = 0.15
angle_x = d * np.cos(0.5/denom*np.pi)
angle_y = d * np.sin(0.5/denom*np.pi)

 この例では、なす角 \thetaを示すために、軸目盛1つ分の角マークを表示することにします。
 そこで、角マーク用のラジアン 0 \leq t \leq \frac{\pi}{n}angle_t_valsとして、円弧の座標を計算してangle_***_valsとします。

 角ラベルは、角マークの中点に表示することにします。
 そこで、 r \sin \frac{\pi}{2 n} r \cos \frac{\pi}{2 n}を計算してangle_***とします。

 円周上に変数(角度)軸を重ねて表示したグラフを作成します。

# タイトル用の文字列(座標の計算式)を作成
coord_label = '$x = r\ \\cos \\theta$\n'
coord_label += '$y = r\ \\sin \\theta$'

# グラフサイズ用の値を指定
axis_size = r * 1.5

# 円周上の変数(ラジアン)目盛を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white')
ax.plot(x, y, color='black') # 円
ax.plot([[0, 0], [1, 0]], [[0, 0], [0, 1]],
        color='black', linewidth=2) # x,y軸
ax.plot(np.stack([np.zeros_like(tick_vals), tick_x_vals], axis=1).T, 
        np.stack([np.zeros_like(tick_vals), tick_y_vals], axis=1).T, 
        color='gray', linestyle=':') # θ目盛グリッド線
for i in range(len(tick_vals)):
    ax.annotate(xy=[tick_x_vals[i], tick_y_vals[i]], 
                text='|', rotation=tick_vals[i]*30.0+90.0, ha='center', va='center') # θ目盛線
    d = 1.15
    ax.text(tick_x_vals[i]*d, tick_y_vals[i]*d, 
            s='$\\frac{' + str(int(tick_vals[i])) + '}{' + str(int(denom)) + '} \pi$', 
            size=15, ha='center', va='center') # θ目盛ラベル
ax.plot(angle_x_vals, angle_y_vals, 
        color='black', lw=1) # θの角マーク
ax.text(angle_x, angle_y, 
        s='$\\theta$', size=15, ha='center', va='center') # θの角ラベル
ax.text(0.5*r, 0.0, 
        s='$r$', size=15, ha='center', va='bottom') # 半径rラベル
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('$r=' + str(r) + '$\n' + coord_label, loc='left')
fig.suptitle('circular coordinates', fontsize=20)
ax.grid()
ax.set_xlim(left=-axis_size, right=axis_size)
ax.set_ylim(bottom=-axis_size, top=axis_size)
ax.set_aspect('equal')
plt.show()

円座標系のグラフ

 目盛ラベルや角ラベルをAxes.text()で描画します。
 x軸・y軸・z軸を黒色の線分として描画します。

 x軸線から反時計回りの角度(なす角)が \thetaに対応しているのが分かります。時計回りの角度は負の値で表現できます。

 円周上の点を変化させたアニメーションで確認します。

・作図コード(クリックで展開)

# フレーム数を指定
frame_num = 120

# ラジアンとして利用する値を作成
t_n = np.linspace(start=-2.0*np.pi, stop=2.0*np.pi, num=frame_num+1)[:frame_num]

# グラフオブジェクトを初期化
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white')
fig.suptitle('circular coordinates', fontsize=20)

# 作図処理を関数として定義
def update(i):
    
    # 前フレームのグラフを初期化
    plt.cla()

    # i番目のラジアンを取得
    t = t_n[i]
    
    # 円周上の点の座標を計算
    point_x = r * np.cos(t)
    point_y = r * np.sin(t)
    
    # 角マーク用のラジアンを作成
    angle_t_vals = np.linspace(start=0.0, stop=t, num=100)

    # 角マーク用の座標を計算
    d = 0.1
    angle_x_vals = d * np.cos(angle_t_vals)
    angle_y_vals = d * np.sin(angle_t_vals)

    # 角ラベル用の座標を計算
    d = 0.15
    angle_x = d * np.cos(0.5*t)
    angle_y = d * np.sin(0.5*t)
    
    # タイトル用の文字列(座標の計算式)を作成
    coord_label = '$r=' + str(r) + '$\n'
    coord_label += '$\\theta^{\circ}='+str((t*180.0/np.pi).round(1))+'^{\circ}, ' + '\\theta='+str((t/np.pi).round(2))+'\pi$\n'
    coord_label += '$x = r\ \\cos \\theta, y = r\ \\sin \\theta$'
    
    # 円周上の点を作図
    ax.plot(x, y, color='black') # 円
    ax.plot([[0.0, 0.0], [1.0, 0.0]], [[0.0, 0.0], [0.0, 1.0]],
            color='black', linewidth=2) # x,y軸
    ax.plot(np.stack([np.zeros_like(tick_vals), tick_x_vals], axis=1).T, 
            np.stack([np.zeros_like(tick_vals), tick_y_vals], axis=1).T, 
            color='gray', linestyle=':') # θ目盛グリッド線
    for i in range(len(tick_vals)):
        ax.annotate(xy=[tick_x_vals[i], tick_y_vals[i]], 
                    text='|', rotation=tick_vals[i]*30.0+90.0, ha='center', va='center') # θ目盛線
        d = 1.15
        ax.text(tick_x_vals[i]*d, tick_y_vals[i]*d, 
                s='$\\frac{' + str(int(tick_vals[i])) + '}{' + str(int(denom)) + '} \pi$', 
                size=15, ha='center', va='center') # θ目盛ラベル
    ax.scatter(point_x, point_y, 
               s=50, 
               label='('+str(point_x.round(2))+', '+str(point_y.round(2))+')', zorder=10) # 円周上の点
    ax.plot([0.0, point_x], [0.0, point_y]) # 原点と円周上の点の線分
    ax.plot(angle_x_vals, angle_y_vals, 
            color='black', lw=1) # θの角マーク
    ax.text(angle_x, angle_y, 
            s='$\\theta$', size=15, ha='center', va='center') # θの角ラベル
    ax.text(0.5*r, 0.0, 
            s='$r$', size=15, ha='center', va='bottom') # 半径rラベル
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_title(coord_label, loc='left')
    ax.grid()
    ax.legend(loc='upper left')
    ax.set_xlim(left=-axis_size, right=axis_size)
    ax.set_ylim(bottom=-axis_size, top=axis_size)
    ax.set_aspect('equal')

# gif画像を作成
ani = FuncAnimation(fig=fig, func=update, frames=frame_num, interval=100)

# gif画像を保存
ani.save('circle_axis_2d.gif')

 作図処理をupdate()として定義して、FuncAnimation()でgif画像を作成します。

円座標系

 弧度法の角度 \thetaと度数法の角度 \theta^{\circ}の関係 \theta = \frac{\pi}{180^{\circ}} \theta^{\circ}を確認できます。

 この記事では、円のグラフを作成しました。次の記事では、球面のグラフを作成します。

おわりに

 しばらく前にRでやったことをいい機会なのでPythonでもやっておきました。
 円周上のラジアン目盛線を'|'で書くのが個人的萌えポイント。なんか少しズレてるけど。

 サイン・コサイン・なんだっけ?となった方はこちらの記事をどうぞ。他の内容もPythonでやってみる機会はなさそうですが。

www.anarchive-beta.com

 あるいはこちらのドラマをどうぞ。

video.unext.jp

 (タイトルとサムネが表示されてほしい…)
 そして4月12日は、元℃-ute・Buono!の鈴木愛理さんの29歳のお誕生日です!

 なんかもう色々凄いのでぜひ聴いてみてください。ちなみに上のドラマのエンディングテーマだったりします。

【次の内容】

www.anarchive-beta.com