はじめに

　『スタンフォード　ベクトル・行列からはじめる最適化数学』の学習ノートです。
　「数式の行間埋め」や「Pythonを使っての再現」によって理解を目指します。本と一緒に読んでください。

　この記事は3.4節「角度」の内容です。
　1つのベクトルのなす角のグラフを作成します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
ベクトルとx軸のなす角の可視化
- 2次元の場合
- 3次元の場合
参考書籍
おわりに

ベクトルとx軸のなす角の可視化

　ベクトルのなす角をグラフで確認します。
　この記事では、任意のベクトル $\mathbf{x}$ とx軸(標準単位ベクトル $\mathbf{e}_1$ )のなす角のグラフを作成します。
　「【Python】3.4：2つのベクトルのなす角の可視化【『スタンフォード線形代数入門』のノート】 - からっぽのしょこ」では、2つの任意のベクトルのなす角のグラフを作成します。

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

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

2次元の場合

　まずは、2次元空間(平面)上でベクトルのなす角を可視化します。

　2次元のベクトルと係数を指定します。また、標準単位ベクトルを作成します。

# 標準単位ベクトル(x軸)を設定:(固定)
e1 = np.array([1.0, 0.0])

# ベクトルを指定
x = np.array([-1.0, -2.0])

# 係数を指定
beta = 1.0

　 $\mathbf{e}_1 = (1, 0)^{\top}$ をe1、 $\mathbf{x} = (x_1, x_2)^{\top}$ をxとして値を指定します。ただし、Pythonでは0からインデックスが割り当てられるので、 $x_1$ の値はx[0]に対応します。

　2つのベクトルのなす角を計算します。

# ベクトルxのなす角を計算
norm_x = np.linalg.norm(beta * x)
theta = np.arccos(beta*x[0] / norm_x)
print(theta) # ラジアン
print(theta/ np.pi) # ラジアンの目安
print(theta * 180.0 / np.pi) # 角度

2.0344439357957027
0.6475836176504333
116.56505117707799

　 $\mathbf{e}_1$ と $\beta \mathbf{x}$ のなす角 $\theta$ は、次の式で定義されます。

$\displaystyle \theta = \arccos \left( \frac{ \beta \mathbf{e}_1^{\top} \mathbf{x} }{ \|\mathbf{e}_1\| \|\beta \mathbf{x}\| } \right)$

　標準単位ベクトル $\mathbf{e}_i$ とのなす角の場合は、 $\mathbf{e}_i^{\top} \mathbf{x} = x_i$ 、 $|\mathbf{e}_i^{\top}| = 1$ なので、次の式で計算できます。

$\displaystyle \theta = \arccos \left( \frac{ \beta x_1 }{ \|\beta \mathbf{x}\| } \right)$

　コサイン関数 $\cos x$ の逆関数(逆コサイン関数) $\arccos x$ はnp.arccos()、ユークリッドノルム $|\mathbf{x}|$ はnp.linalg.norm()で計算できます。
　ただし、 $\beta \mathbf{x}$ が0ベクトルだと0除算になるため定義できません(計算結果がnanになります)。

　 $\theta$ は、弧度法における角度であり、ラジアンと呼びます。度数法における角度を $\theta^{\circ}$ で表すことにします。
　 $\theta$ と $\theta^{\circ}$ は、 $\theta = \theta^{\circ} \frac{2 \pi}{360}$ 、 $\theta^{\circ} = \theta \frac{360}{2 \pi}$ で変換できます。よって、 $\theta^{\circ} = 0^{\circ}$ のとき $\theta = 0$ であり、 $\theta^{\circ} = 180^{\circ}$ のとき $\theta = \pi$ です。
　円周率 $\pi$ はnp.piで扱えます。

　2つのベクトル(線分)によってできる角は、劣角 $0^{\circ} \leq \theta^{\circ} \leq 180^{\circ}$ と優角 $360^{\circ} - \theta^{\circ}$ の2つです。ラジアンで表すと、劣角 $0 \leq \theta \leq \pi$ と優角 $2 \pi - \theta$ です。
　なす角 $\theta$ は劣角の角度に対応します。

　なす角を示す角マークを描画するための座標を作成します。なす角の計算自体には不要な処理です。

# 角マーク用の値(ラジアン)を作成
sgn = 1.0 if beta*x[1] >= 0.0 else -1.0
rad_vals = np.linspace(start=0.0, stop=sgn*theta, num=100)
print(rad_vals[:5])

# 角マークの座標を計算
r = 0.2
angle_x_vals = np.cos(rad_vals) * r
angle_y_vals = np.sin(rad_vals) * r
print(angle_x_vals[:5])
print(angle_y_vals[:5])

# 角ラベルの座標を計算
r = 0.3
angle_x = np.cos(0.5 * sgn*theta) * r
angle_y = np.sin(0.5 * sgn*theta) * r
print(angle_x)
print(angle_y)

[ 0.         -0.02054994 -0.04109988 -0.06164982 -0.08219975]
[0.2        0.19995777 0.1998311  0.19962005 0.1993247 ]
[ 0.         -0.0041097  -0.00821766 -0.01232215 -0.01642144]
0.15771933363574006
-0.25519524250561193

　 $0$ から $\theta$ までのラジアン $t$ を作成して、角マークの座標を計算します。
　ただし、第1・2象限にできる角が優角のとき(180°より大きい角度で第3象限にまで跨るとき)、言い換えると、 $\beta x_2$ が負の値のとき(ベクトル $\beta \mathbf{x}$ が第3・4象限のとき)、第3・4象限にできる劣角を $- \theta$ として扱うことにします。(完成図を見てもらえば言いたいことが伝わると思います。)
　そこでbeta * x[1]が、0.0以上であればsgnを1.0、0.0未満であればsgnを-1.0にして、0からsgn * thetaまでの値を作成してrad_valsとします。

　rad_valsの各要素を円形の座標に変換します。
　座標 $x, y$ は、半径を $r$ として、それぞれ次の式で計算できます。

$\displaystyle \left\{ \begin{aligned} x &= r \cos t \\ y &= r \sin t \end{aligned} \right.$

　半径rで角マークのサイズを調整します。

　角ラベルは、 $x = r \cos \frac{\theta}{2}$ 、 $y = r \sin \frac{\theta}{2}$ の座標(中点)に表示することにします。

　2次元空間上にベクトルのなす角を作図します。

# 作図用の値を設定
x_min = np.floor(np.min([0.0, e1[0], beta*x[0]])) - 1
x_max = np.ceil(np.max([0.0, e1[0], beta*x[0]])) + 1
y_min = np.floor(np.min([0.0, e1[1], beta*x[1]])) - 1
y_max = np.ceil(np.max([0.0, e1[1], beta*x[1]])) + 1

# 2Dベクトルのなす角を作図
fig, ax = plt.subplots(figsize=(6, 6), facecolor='white')
ax.quiver(0, 0, *e1, ec='black', linewidth=1.5, 
          scale_units='xy', scale=1, units='dots', width=0.1, headwidth=0.1, 
          label='$e_1 = ('+', '.join(map(str, e1))+')$') # 標準単位ベクトルe1
ax.quiver(0, 0, *beta*x, color='red', 
          angles='xy', scale_units='xy', scale=1, 
          label='$\\beta x = \\beta ('+', '.join(map(str, x))+')$') # ベクトルβx
ax.plot(angle_x_vals, angle_y_vals, color='black', linewidth=1) # 角マーク
ax.text(angle_x, angle_y, s='$\\theta$', ha='center', va='center') # 角ラベル
ax.set_xticks(ticks=np.arange(x_min, x_max+1))
ax.set_yticks(ticks=np.arange(y_min, y_max+1))
ax.grid()
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_title('$\\beta='+str(beta)+', ' + 
             '\\theta='+str((theta / np.pi).round(2))+'\pi, ' + 
             '\\theta^\circ='+str((theta * 180.0/np.pi).round(1))+'$', loc='left')
fig.suptitle('angle x', fontsize=20)
ax.legend(loc='upper left')
ax.set_aspect('equal')
plt.show()

　axes.quiver()でベクトルを描画します。第1・2引数に始点の座標、第3・4引数にベクトルのサイズ(変化量)を指定します。その他に、指定した値の通りに(調整せずに)矢印や線分を描画するための設定をしています。
　配列e1, beta*xの前に*を付けてアンパック(展開)して指定しています。

　係数の値を変えてみます。

　係数の大小はなす角に影響せず、係数の正負はベクトルの向きに影響するためなす角に影響するのが分かります。

　ベクトルの値を変化させたアニメーションを作成します。

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

# フレーム数を設定
frame_num = 60

# 標準単位ベクトル(x軸)を設定:(固定)
e1 = np.array([1.0, 0.0])

# ベクトルとして利用する値を指定
theta_vals = np.linspace(start=0.0, stop=2.0*np.pi, num=frame_num+1)[:frame_num]
x_vals = np.array(
    [np.cos(theta_vals), 
     np.sin(theta_vals)]
).T

# 係数を指定
beta = 1.0

# 作図用の値を設定
x_min = np.floor(np.min([0.0, e1[0], *beta*x_vals[:, 0]])) - 1
x_max = np.ceil(np.max([0.0, e1[0], *beta*x_vals[:, 0]])) + 1
y_min = np.floor(np.min([0.0, e1[1], *beta*x_vals[:, 1]])) - 1
y_max = np.ceil(np.max([0.0, e1[1], *beta*x_vals[:, 1]])) + 1

# 作図用のオブジェクトを初期化
fig, ax = plt.subplots(figsize=(6, 6), facecolor='white')
fig.suptitle('angle x', fontsize=20)

# 作図処理を関数として定義
def update(i):
    
    # 前フレームのグラフを初期化
    plt.cla()
    
    # i番目のベクトルを取得
    x = x_vals[i]
    
    # ベクトルxのなす角を計算
    norm_x = np.linalg.norm(beta * x)
    theta = np.arccos(beta*x[0] / norm_x)
    
    # 角マーク用の値(ラジアン)を作成
    sgn = 1.0 if beta*x[1] >= 0.0 else -1.0
    rad_vals = np.linspace(start=0.0, stop=sgn * theta, num=100)
    
    # 角マークの座標を計算
    r = 0.2
    angle_x_vals = np.cos(rad_vals) * r
    angle_y_vals = np.sin(rad_vals) * r
    
    # 角ラベルの座標を計算
    r = 0.3
    angle_x = np.cos(0.5 * sgn*theta) * r
    angle_y = np.sin(0.5 * sgn*theta) * r

    # 2Dベクトルのなす角を作図
    ax.quiver(0, 0, *e1, ec='black', linewidth=1.5, 
              scale_units='xy', scale=1, units='dots', width=0.1, headwidth=0.1, 
              label='$e_1 = ('+', '.join(map(str, e1))+')$') # 標準単位ベクトルe1
    ax.quiver(0, 0, *beta*x, color='red', 
              angles='xy', scale_units='xy', scale=1, 
              label='$\\beta x = \\beta ('+', '.join(map(str, x.round(2)))+')$') # ベクトルβx
    ax.plot(angle_x_vals, angle_y_vals, color='black', linewidth=1) # 角マーク
    ax.text(angle_x, angle_y, s='$\\theta$', ha='center', va='center') # 角ラベル
    ax.set_xticks(ticks=np.arange(x_min, x_max+1))
    ax.set_yticks(ticks=np.arange(y_min, y_max+1))
    ax.grid()
    ax.set_xlabel('$x_1$')
    ax.set_ylabel('$x_2$')
    ax.set_title('$\\beta='+str(beta)+', ' + 
                 '\\theta='+str((theta / np.pi).round(2))+'\pi, ' + 
                 '\\theta^\circ='+str((theta * 180.0/np.pi).round(1))+'$', loc='left')
    ax.legend(loc='upper left')
    ax.set_aspect('equal')

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

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

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

　この例では、角マークの座標の処理と同様にして、単位円上を通るように $\mathbf{x}$ の値を指定しています。
　円周は $2 \pi$ で1周期なので、 $2 \pi$ の $n$ 倍の範囲(サイズ)をフレーム数+1個の等間隔の値に切り分けて最大値を除くと、最後のフレームと最初のフレームの繋がりが良くなります。

3次元の場合

　次は、3次元空間上でベクトルのなす角を可視化します。

　3次元ベクトルと係数を指定します。また、標準単位ベクトルを作成します。

# 標準単位ベクトル(x軸)を設定:(固定)
e1 = np.array([1.0, 0.0, 0.0])

# ベクトルを指定
x = np.array([-1.0, 0.1, 1.5])

# 係数を指定
beta = 1.0

　 $\mathbf{e}_1 = (1, 0, 0)^{\top}$ をe1、 $\mathbf{x} = (x_1, x_2, x_3)^{\top}$ をxとして値を指定します。

　2つのベクトルのなす角を計算します。

# ベクトルxのなす角を計算
norm_x = np.linalg.norm(beta * x)
theta = np.arccos(beta*x[0] / norm_x)
print(theta) # ラジアン
print(theta/ np.pi) # ラジアンの目安
print(theta * 180.0 / np.pi) # 角度

2.1577759987460636
0.6868414325709747
123.63145786277545

　「2次元の場合」のコードで計算できます。

　参考として、xy面(グラフの底面)にx軸とy軸におけるなす角を示す角マークを描画するための座標を作成します。

# xy面上のなす角を計算
sgn_y = 1.0 if beta*x[1] >= 0.0 else -1.0
theta_xy = sgn_y * np.arccos(beta*x[0] / np.linalg.norm(beta * x[[0, 1]]))
print(theta_xy)

# xy面上の角マーク用の値(ラジアン)を作成
rad_xy_vals = np.linspace(start=0.0, stop=theta_xy, num=100)

# xy面上の角マーク用の座標を計算
r = 0.2
angle_xy_x_vals = np.cos(rad_xy_vals) * r
angle_xy_y_vals = np.sin(rad_xy_vals) * r

# xy面上の角ラベルの座標を計算
r = 0.3
angle_xy_x = np.cos(0.5 * theta_xy) * r
angle_xy_y = np.sin(0.5 * theta_xy) * r

3.0419240010986326

　「2次元の場合」のときと同様に処理します。

　同様に、xz面(デフォルトのグラフの奥面)にx軸とz軸におけるなす角についても作成します。

# xz面上のなす角を計算
sgn_z = 1.0 if beta*x[2] >= 0.0 else -1.0
theta_xz = sgn_z * np.arccos(beta*x[0] / np.linalg.norm(beta * x[[0, 2]]))
print(theta_xz)

# xz面上の角マーク用の値(ラジアン)を作成
rad_xz_vals = np.linspace(start=0.0, stop=theta_xz, num=100)

# xz面上の角マーク用の座標を計算
r = 0.2
angle_xz_x_vals = np.cos(rad_xz_vals) * r
angle_xz_z_vals = np.sin(rad_xz_vals) * r

# xz面上の角ラベルの座標を計算
r = 0.3
angle_xz_x = np.cos(0.5 * theta_xz) * r
angle_xz_z = np.sin(0.5 * theta_xz) * r

2.1587989303424644

　xy面のときと同様に処理します。

　空間上に、なす角を示す角マークを描画するための座標を作成します。この処理も、なす角の計算自体には不要です。

# 角マーク用の値(球面座標)を計算
psi = np.arccos(beta*x[2] / np.linalg.norm(beta * x))
sgn_y = 1.0 if beta*x[1] >= 0.0 else -1.0
phi = sgn_y * np.arccos(beta*x[0] / np.linalg.norm(beta * x[[0, 1]]))
print(psi)
print(phi)

# 角マーク用の値(ラジアン)を作成
rad_psi_vals = np.linspace(start=0.5*np.pi, stop=psi, num=100)
rad_phi_vals = np.linspace(start=0.0, stop=phi, num=100)

# 角マーク用の座標を計算
r = 0.2
angle_x_vals = np.sin(rad_psi_vals) * np.cos(rad_phi_vals) * r
angle_y_vals = np.sin(rad_psi_vals) * np.sin(rad_phi_vals) * r
angle_z_vals = np.cos(rad_psi_vals) * r
#angle_z_vals = np.linspace(start=0.0, stop=np.cos(psi), num=100) * r

# 角ラベル用の座標を計算
r = 0.4
angle_x = np.sin(np.median(rad_psi_vals)) * np.cos(0.5 * phi) * r
angle_y = np.sin(np.median(rad_psi_vals)) * np.sin(0.5 * phi) * r
angle_z = np.cos(np.median(rad_psi_vals)) * r

0.5903010240027555
3.0419240010986326

　z軸( $\mathbf{e}_3$ )と $\beta \mathbf{x}$ のなす角を $\psi$ 、x軸( $\mathbf{e}_1$ )と2次元ベクトル $\beta (x_1, x_2)^{\top}$ のなす角を $\phi$ とします。2つの角度を使って、なす角 $\theta$ の角マークの座標を計算します。(詳しくは「球面座標系」で調べてみてください。)ちなみに、 $\phi$ はtheta_xです。
　角度 $\psi, \phi$ は、それぞれ次の式で計算できます。

$\displaystyle \begin{aligned} \psi &= \arccos \left( \frac{ \beta x_3 }{ \|\beta \mathbf{x}\| } \right) \\ \phi &= \mathrm{sgn}(x_2) \arccos \left( \frac{ \beta x_1 }{ \|\beta (x_1, x_2)\| } \right) \\ \mathrm{sgn}(x) &= \begin{cases} 1 & (x \geq 0) \\ -1 & (x < 0) \end{cases} \end{aligned}$

　ただし、 $\beta \mathbf{x}$ や $\beta (x_1, x_2)$ が0ベクトルだと0除算になるため定義できません(計算結果がnanになります)。

　 $\frac{1}{2} \pi$ から $\psi$ までのラジアン $t$ と、 $0$ から $\phi$ までのラジアン $u$ を作成して、角マークの座標を計算します。 $\frac{1}{2} \pi$ と $0$ は、 $\mathbf{e}_3$ に対する $\psi, \phi$ の値です。
　 $t$ の値rad_psi_vals、 $u$ の値rad_phi_valsを作成して、各要素を球形の座標に変換します。
　座標 $x, y, z$ は、半径を $r$ として、それぞれ次の式で計算できます。

$\displaystyle \left\{ \begin{aligned} x &= \sin t \cos u \\ y &= \sin t \sin u \\ z &= \cos t \end{aligned} \right.$

　半径rで角マークのサイズを調整します。

　角ラベルは、 $x = \sin \frac{1}{2} (\frac{\pi}{2} + \psi) \cos \frac{1}{2} \phi$ 、 $y = \sin \frac{1}{2} (\frac{\pi}{2} + \psi) \sin \phi$ 、 $z = \cos \frac{\psi}{2}$ の座標(中点)に表示することにします。

　3次元空間上にベクトルのなす角を作図します。

# 作図用の値を設定
x_min = np.floor(np.min([0.0, e1[0], beta*x[0]])) - 1
x_max = np.ceil(np.max([0.0, e1[0], beta*x[0]])) + 1
y_min = np.floor(np.min([0.0, e1[1], beta*x[1]])) - 1
y_max = np.ceil(np.max([0.0, e1[1], beta*x[1]])) + 1
z_min = np.floor(np.min([0.0, e1[2], beta*x[2]])) - 1
z_max = np.ceil(np.max([0.0, e1[2], beta*x[2]])) + 1

# 3Dベクトルのなす角を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
ax.quiver(0, 0, 0, *e1, color='black', arrow_length_ratio=0.0, 
          label='$e_1 = ('+', '.join(map(str, e1))+')$') # 標準単位ベクトルe1
ax.quiver(0, 0, 0, *beta*x, color='red', arrow_length_ratio=0.1, 
          label='$\\beta x = \\beta ('+', '.join(map(str, x.round(2)))+')$') # ベクトルβx
ax.quiver(0, y_max, 0, *e1, color='black', linestyle='--', arrow_length_ratio=0.0) # xy面上のベクトルe1
ax.quiver(0, 0, z_min, *e1, color='black', linestyle='--', arrow_length_ratio=0.0) # xz面上のベクトルe1
ax.quiver(0, y_max, 0, beta*x[0], 0, beta*x[2], color='red', linestyle='--', arrow_length_ratio=0.0) # xy面上のベクトルx
ax.quiver(0, 0, z_min, beta*x[0], beta*x[1], 0, color='red', linestyle='--', arrow_length_ratio=0.0) # xz面上のベクトルx
ax.plot(angle_xy_x_vals, angle_xy_y_vals, np.repeat(z_min, len(angle_xy_x_vals)), 
        color='black', linewidth=1) # xy面上の角マーク
ax.plot(angle_xz_x_vals, np.repeat(y_max, len(angle_xz_x_vals)), angle_xz_z_vals, 
        color='black', linewidth=1) # xz面上の角マーク
ax.plot(angle_x_vals, angle_y_vals, angle_z_vals, color='black', linewidth=1) # 角マーク
ax.text(angle_x, angle_y, angle_z, s='$\\theta$', size=15, ha='center', va='center') # 角ラベル
ax.quiver([0, e1[0], beta*x[0], 0, e1[0], beta*x[0]], 
          [0, e1[1], beta*x[1], y_max, y_max, y_max], 
          [z_min, z_min, z_min, 0, e1[2], beta*x[2]], 
          [0, 0, 0, 0, 0, 0], 
          [0, 0, 0, -y_max, e1[1]-y_max, beta*x[1]-y_max], 
          [-z_min, e1[2]-z_min, beta*x[2]-z_min, 0, 0, 0], 
          color='gray', arrow_length_ratio=0, linestyle=':') # 補助線
ax.set_xticks(ticks=np.arange(x_min, x_max+1))
ax.set_yticks(ticks=np.arange(y_min, y_max+1))
ax.set_zticks(ticks=np.arange(z_min, z_max+1))
ax.set_ylim(y_min, y_max)
ax.set_zlim(z_min, z_max)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$x_3$')
ax.set_title('$\\beta='+str(beta)+', ' + 
             '\\theta='+str((theta / np.pi).round(2))+'\pi, ' + 
             '\\theta^\circ='+str((theta * 180.0 / np.pi).round(1))+'$', loc='left')
fig.suptitle('angle x', fontsize=20)
ax.legend(loc='upper left')
ax.set_aspect('equal')
#ax.view_init(elev=90, azim=270) # xy面
#ax.view_init(elev=0, azim=270) # xz面
plt.show()

　axes.quiver()の第1・2・3引数に始点の座標、第4・5・6引数にベクトルのサイズを指定します。原点を始点とします。

　係数の値を変えてみます。

　ベクトルの値を変化させたアニメーションを作成します。

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

# フレーム数を設定
frame_num = 300

# 標準単位ベクトル(x軸)を設定:(固定)
e1 = np.array([1.0, 0.0, 0.0])

# ベクトルとして利用する値を指定
psi_vals = np.linspace(start=0.5*np.pi, stop=2.5*np.pi, num=frame_num+1)[:frame_num]
phi_vals = np.linspace(start=-2.0*np.pi, stop=4.0*np.pi, num=frame_num+1)[:frame_num]
x_vals = np.array(
    [np.sin(psi_vals) * np.cos(phi_vals), 
     np.sin(psi_vals) * np.sin(phi_vals), 
     np.cos(psi_vals)]
).T

# 係数を指定
beta = 1.0

# 作図用の値を設定
x_min = np.floor(np.min([0.0, e1[0], *beta*x_vals[:, 0]])) - 1
x_max = np.ceil(np.max([0.0, e1[0], *beta*x_vals[:, 0]])) + 1
y_min = np.floor(np.min([0.0, e1[1], *beta*x_vals[:, 1]])) - 1
y_max = np.ceil(np.max([0.0, e1[1], *beta*x_vals[:, 1]])) + 1
z_min = np.floor(np.min([0.0, e1[2], *beta*x_vals[:, 2]])) - 1
z_max = np.ceil(np.max([0.0, e1[2], *beta*x_vals[:, 2]])) + 1

# 作図用のオブジェクトを初期化
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
fig.suptitle('angle x', fontsize=20)

# 作図処理を関数として定義
def update(i):
    
    # 前フレームのグラフを初期化
    plt.cla()
    
    # i番目のベクトルを取得
    x = x_vals[i]
    
    # ベクトルxのなす角を計算
    theta = np.arccos(beta*x[0] / np.linalg.norm(beta * x))
    
    # xy面・xz面上のなす角を計算
    sgn_y = 1.0 if beta*x[1] >= 0.0 else -1.0
    theta_xy = sgn_y * np.arccos(beta*x[0] / np.linalg.norm(beta * x[[0, 1]]))
    sgn_z = 1.0 if beta*x[2] >= 0.0 else -1.0
    theta_xz = sgn_z * np.arccos(beta*x[0] / np.linalg.norm(beta * x[[0, 2]]))
    
    # xy面・xz面上の角マーク用の値(ラジアン)を作成
    rad_xy_vals = np.linspace(start=0.0, stop=theta_xy, num=100)
    rad_xz_vals = np.linspace(start=0.0, stop=theta_xz, num=100)
    
    # xy面・xz面上の角マーク用の座標を計算
    r = 0.2
    angle_xy_x_vals = np.cos(rad_xy_vals) * r
    angle_xy_y_vals = np.sin(rad_xy_vals) * r
    angle_xz_x_vals = np.cos(rad_xz_vals) * r
    angle_xz_z_vals = np.sin(rad_xz_vals) * r
    
    # xy面・xz面上の角ラベルの座標を計算
    r = 0.3
    angle_xy_x = np.cos(0.5 * theta_xy) * r
    angle_xy_y = np.sin(0.5 * theta_xy) * r
    angle_xz_x = np.cos(0.5 * theta_xz) * r
    angle_xz_z = np.sin(0.5 * theta_xz) * r
    
    # 角マーク用の値(球面座標)を計算
    psi = np.arccos(beta*x[2] / np.linalg.norm(beta * x))
    sgn_y = 1.0 if beta*x[1] >= 0.0 else -1.0
    phi = sgn_y * np.arccos(beta*x[0] / np.linalg.norm(beta * x[[0, 1]]))
    
    # 角マーク用の値(ラジアン)を作成
    rad_psi_vals = np.linspace(start=0.5*np.pi, stop=psi, num=100)
    rad_phi_vals = np.linspace(start=0.0, stop=phi, num=100)

    # 角マーク用の座標を計算
    r = 0.2
    angle_x_vals = np.sin(rad_psi_vals) * np.cos(rad_phi_vals) * r
    angle_y_vals = np.sin(rad_psi_vals) * np.sin(rad_phi_vals) * r
    angle_z_vals = np.cos(rad_psi_vals) * r
    
    # 角ラベル用の座標を計算
    r = 0.4
    angle_x = np.sin(np.median(rad_psi_vals)) * np.cos(0.5 * phi) * r
    angle_y = np.sin(np.median(rad_psi_vals)) * np.sin(0.5 * phi) * r
    angle_z = np.cos(np.median(rad_psi_vals)) * r
    
    # 3Dベクトルのなす角を作図
    ax.quiver(0, 0, 0, *e1, color='black', arrow_length_ratio=0.0, 
              label='$e_1 = ('+', '.join(map(str, e1))+')$') # 標準単位ベクトルe1
    ax.quiver(0, 0, 0, *beta*x, color='red', arrow_length_ratio=0.1, 
              label='$\\beta x = \\beta ('+', '.join(map(str, x.round(2)))+')$') # ベクトルβx
    ax.quiver(0, y_max, 0, *e1, color='black', linestyle='--', arrow_length_ratio=0.0) # xy面上のベクトルe1
    ax.quiver(0, 0, z_min, *e1, color='black', linestyle='--', arrow_length_ratio=0.0) # xz面上のベクトルe1
    ax.quiver(0, y_max, 0, beta*x[0], 0, beta*x[2], color='red', linestyle='--', arrow_length_ratio=0.0) # xy面上のベクトルx
    ax.quiver(0, 0, z_min, beta*x[0], beta*x[1], 0, color='red', linestyle='--', arrow_length_ratio=0.0) # xz面上のベクトルx
    ax.plot(angle_xy_x_vals, angle_xy_y_vals, np.repeat(z_min, len(angle_xy_x_vals)), 
            color='black', linewidth=1) # xy面上の角マーク
    ax.plot(angle_xz_x_vals, np.repeat(y_max, len(angle_xz_x_vals)), angle_xz_z_vals, 
            color='black', linewidth=1) # xz面上の角マーク
    ax.plot(angle_x_vals, angle_y_vals, angle_z_vals, color='black', linewidth=1) # 角マーク
    ax.text(angle_x, angle_y, angle_z, s='$\\theta$', size=15, ha='center', va='center') # 角ラベル
    ax.quiver([0, e1[0], beta*x[0], 0, e1[0], beta*x[0]], 
              [0, e1[1], beta*x[1], y_max, y_max, y_max], 
              [z_min, z_min, z_min, 0, e1[2], beta*x[2]], 
              [0, 0, 0, 0, 0, 0], 
              [0, 0, 0, -y_max, e1[1]-y_max, beta*x[1]-y_max], 
              [-z_min, e1[2]-z_min, beta*x[2]-z_min, 0, 0, 0], 
              color='gray', arrow_length_ratio=0, linestyle=':') # 補助線
    ax.set_xticks(ticks=np.arange(x_min, x_max+1))
    ax.set_yticks(ticks=np.arange(y_min, y_max+1))
    ax.set_zticks(ticks=np.arange(z_min, z_max+1))
    ax.set_ylim(y_min, y_max)
    ax.set_zlim(z_min, z_max)
    ax.set_xlabel('$x_1$')
    ax.set_ylabel('$x_2$')
    ax.set_zlabel('$x_3$')
    ax.set_title('$\\beta='+str(beta)+', ' + 
                 '\\theta='+str((theta / np.pi).round(2))+'\pi, ' + 
                 '\\theta^\circ='+str((theta * 180.0 / np.pi).round(1))+'$', loc='left')
    ax.legend(loc='upper left')
    ax.set_aspect('equal')
    #ax.view_init(elev=90, azim=270) # xy面
    #ax.view_init(elev=0, azim=270) # xz面

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

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

　「2次元の場合」と同様にして、アニメーションを作成します。

　ちなみに、点 $\beta \mathbf{x}$ (矢印の先端)は、次の軌道をしています。

# ベクトルの軌道を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
ax.plot(*beta*x_vals.T, color='red') # ベクトルβxの軌道
ax.set_xticks(ticks=np.arange(x_min, x_max+1))
ax.set_yticks(ticks=np.arange(y_min, y_max+1))
ax.set_zticks(ticks=np.arange(z_min, z_max+1))
ax.set_ylim(y_min, y_max)
ax.set_zlim(z_min, z_max)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$x_3$')
ax.set_title('$\\beta='+str(beta)+'$', loc='left')
fig.suptitle('trace $\\beta x$', fontsize=20)
ax.set_aspect('equal')
plt.show()