からっぽのしょこ

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

【Python】球面上の2点を結ぶ円弧を作図したい

はじめに

 調べても分からなかったので自分なりにやってみる黒魔術シリーズです。もっといい方法があれば教えてください。

 この記事では、球面上の2点をいい感じに繋ぐ曲線のグラフを作成します。

【前の内容】

www.anarchive-beta.com

【目次】

球面上の2点を結ぶ円弧を作図したい

 Matplotlibライブラリを利用して、3次元空間における球面(spherical surface)上の2点を結ぶ良い感じの曲線を作成したい。また、その円弧(circular arc)の半径を調整して、2つの3次元ベクトルのなす角(angle between two vectors)を示したい。

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

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


なす角の描画

 まずは、完成図(やりたいこと)を確認する。
 2つの3次元ベクトルのなす角を円弧(角マーク)で示すグラフを作成する。
 なす角については「【Python】3.4:ベクトルとx軸のなす角の可視化【『スタンフォード線形代数入門』のノート】 - からっぽのしょこ」を参照のこと。

 3次元ベクトル \mathbf{a}, \mathbf{b}を作成して、なす角を計算する。

# 3次元ベクトルを指定
a = np.array([1.0, 1.0, 1.0])
b = np.array([-2.0, -0.5, 1.0])

# ベクトルa,bのなす角を計算
cos_theta = np.dot(a, b) / np.linalg.norm(a) / np.linalg.norm(b)
theta = np.arccos(cos_theta)
print(theta) # 0からπの値
print(theta / np.pi) # 0から1の値
1.9583930134500773
0.623375857214425

 2つのベクトル \mathbf{a} = (a_1, a_2, a_3), \mathbf{b} = (b_1, b_2, b_3)a, bとして値を指定する。Pythonではインデックスが0から割り当てられるので、 a_1a[0]に対応することに注意する。

  \mathbf{a}, \mathbf{b}のなす角 \thetaは、次の式で計算できる。

 \displaystyle
\theta
    = \arccos \left(
          \frac{
              \mathbf{a}^{\top} \mathbf{b}
          }{
              \|\mathbf{a}\| \|\mathbf{b}\|
          }
      \right)

  \thetaはラジアン(弧度法における角度)であり、 0 \leq \theta \leq \piの値(度数法による角度 0^{\circ} \leq \theta^{\circ} \leq 180^{\circ}に対応)をとる。ただし、 \mathbf{a}または \mathbf{b}が0ベクトルだと0除算になるため定義できない(計算結果が nan になる)。
 コサイン関数 \cos xの逆関数(逆コサイン関数) \arccos xnp.arccos()、内積 \mathbf{x}^{\top} \mathbf{y}np.dot()、ユークリッドノルム \|\mathbf{x}\|np.linalg.norm()で計算できる。

 角マークと角ラベルの座標を計算するための係数を作成する。

# 角マーク用のラジアンを作成
t = np.linspace(start=0.0, stop=0.5*np.pi, num=100)
print(t[:5])
print(t.shape)

# ラジアンを複製
T = np.stack([t, t, t], axis=1)
print(T[:5])
print(T.shape)

# 角マーク用の係数を作成
Beta1 = np.cos(T)
Beta2 = np.sin(T)

# 角ラベル用のラジアンを作成
t = 0.25 * np.pi

# 角ラベル用の係数を作成
beta1 = np.cos(t)
beta2 = np.sin(t)
print(beta1)
print(beta2)
[0.         0.01586663 0.03173326 0.04759989 0.06346652]
(100,)
[[0.         0.         0.        ]
 [0.01586663 0.01586663 0.01586663]
 [0.03173326 0.03173326 0.03173326]
 [0.04759989 0.04759989 0.04759989]
 [0.06346652 0.06346652 0.06346652]]
(100, 3)
0.7071067811865476
0.7071067811865476

  0 \leq t \leq \frac{\pi}{2}の範囲のラジアン( 0^{\circ} \leq t^{\circ} \leq 90^{\circ}に対応)をtとして、3列分を結合した配列をTとする。
 Tを用いて、線形結合に使う係数 \beta_1 = \cos t, \beta_2 = \sin tを計算してBeta1, Beta2とする。

 角ラベルは角マークの中点に配置することにする。
 そこで、中点のラジアン t = \frac{\pi}{4}tとして、 \beta_1 = \cos t, \beta_2 = \sin tを計算してbeta1, beta2とする。

 角マークと角ラベルの座標を計算する。

# 角マークのサイズ(なす角からの半径)を指定
r = 0.3

# ベクトルa,b間を結ぶ曲線の座標を計算(ベクトルa,bを線形結合)
X = Beta1 * a/np.linalg.norm(a) + Beta2 * b/np.linalg.norm(b)

# 角マークの座標を計算(サイズを調整)
X *= r / np.linalg.norm(X, axis=1, keepdims=True)
print(X[:5])
print(X.shape)
print(np.linalg.norm(X[:5], axis=1))

# 角ラベルの表示位置(なす角からのノルム)を指定
r = 0.5

# ベクトルa,b間を結ぶ曲線の中点の座標を計算(ベクトルa,bを線形結合)
x = beta1 * a/np.linalg.norm(a) + beta2 * b/np.linalg.norm(b)

# 角ラベルの座標を計算(表示位置を調整)
x *= r / np.linalg.norm(x)
print(x)
print(x.shape)
print(np.linalg.norm(x))
[[0.17320508 0.17320508 0.17320508]
 [0.1700513  0.17318617 0.17632103]
 [0.1668212  0.1731285  0.1794358 ]
 [0.16351323 0.17303067 0.18254811]
 [0.16012584 0.17289119 0.18565654]]
(100, 3)
[0.3 0.3 0.3 0.3 0.3]
[-0.13247569  0.16099114  0.45445797]
(3,)
0.5

  \mathbf{a}, \mathbf{b}を正規化し線形結合したベクトル \mathbf{x} = (x_1, x_2, x_3)を計算して、さらに \mathbf{x}を正規化し正の実数 rを掛けたベクトルを \tilde{\mathbf{x}} = (\tilde{x}_1, \tilde{x}_2, \tilde{x}_3)とする。

 \displaystyle
\begin{aligned}
\mathbf{x}
   &= \beta_1 \frac{\mathbf{a}}{\|\mathbf{a}\|}
      + \beta_2 \frac{\mathbf{b}}{\|\mathbf{b}\|}
\\
\tilde{\mathbf{x}}
   &= r \frac{\mathbf{x}}{\|\mathbf{x}\|}
\end{aligned}

 ベクトル \mathbf{x}をそのベクトルのノルム \|\mathbf{x}\|で割るとノルムが1になる \left\|\frac{\mathbf{x}}{\|\mathbf{x}\|}\right\| = 1。この計算をベクトルの正規化と呼ぶ。
 さらに、定数 rを掛けるとノルムが rの絶対値になる \left\|r \frac{\mathbf{x}}{\|\mathbf{x}\|}\right\| = |r|

 曲線用の係数Beta1, Beta2を用いて求めた座標(2次元配列)をX、ラベル用の係数beta1, beta2を用いて求めた座標(1次元配列)をxとする。Xの各行が1つの点 \mathbf{x}に対応する。

 ベクトル \mathbf{a}, \mathbf{b}となす角のグラフを作成する。

# グラフサイズ用の値を設定
x_min = np.floor(np.min([0.0, a[0], b[0]])) - 1.0
x_max =  np.ceil(np.max([0.0, a[0], b[0]])) + 1.0
y_min = np.floor(np.min([0.0, a[1], b[1]])) - 1.0
y_max =  np.ceil(np.max([0.0, a[1], b[1]])) + 1.0
z_min = np.floor(np.min([0.0, a[2], b[2]])) - 1.0
z_max =  np.ceil(np.max([0.0, a[2], b[2]])) + 1.0

# 矢のサイズを指定
l = 0.2

# なす角を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
ax.quiver(0, 0, 0, *a, 
          color='red', arrow_length_ratio=l/np.linalg.norm(a), 
          label='$a=('+', '.join(map(str, a))+')$') # ベクトルa
ax.quiver(0, 0, 0, *b, 
          color='blue', arrow_length_ratio=l/np.linalg.norm(b), 
          label='$b=('+', '.join(map(str, b.round(2)))+')$') # ベクトルb
ax.plot(*X.T, 
        color='black', linewidth=1, zorder=10) # 角マーク
ax.text(*x, 
        s='$\\theta$', size=15, ha='center', va='center', zorder=15) # 角ラベル
ax.quiver([0, a[0], b[0]], 
          [0, a[1], b[1]], 
          [0, a[2], b[2]], 
          [0, 0, 0], 
          [0, 0, 0], 
          [z_min, z_min-a[2], z_min-b[2]], 
          color=['gray', 'red', 'blue'], 
          arrow_length_ratio=0, linestyle=':') # 座標用の補助線
ax.set_xlim(left=x_min, right=x_max)
ax.set_ylim(bottom=y_min, top=y_max)
ax.set_zlim(bottom=z_min, top=z_max)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$x_3$')
ax.set_title('$\\theta=' + str((theta/np.pi).round(2)) + '\pi, ' + 
             '\\theta^{\circ}=' + str((theta/np.pi*180.0).round(1)) + '^{\circ}$', loc='left')
fig.suptitle('$\\theta = \\arccos(\\frac{a \cdot b}{\|a\| \|b\|})$', fontsize=20)
ax.legend()
ax.set_aspect('equal')
#ax.view_init(elev=90, azim=270) # xy面
#ax.view_init(elev=0, azim=270) # xz面
#ax.view_init(elev=0, azim=0) # yz面
plt.show()

ベクトルのなす角のグラフ

 axes.quiver()の第1・2・3引数に始点の座標、第4・5・6引数にベクトルの移動量を指定する。
 各配列の形状に応じて、配列の前に*を付けてアンパック(展開)して座標などを指定している。原点を始点とする。

 グラフを回転させて確認する。

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

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

# 水平方向の角度として利用する値を作成
h_n = np.linspace(start=0.0, stop=360.0, num=frame_num+1)[:frame_num]

# グラフオブジェクトを初期化
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
fig.suptitle('$\\theta = \\arccos(\\frac{a \cdot b}{\|a\| \|b\|})$', fontsize=20)

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

    # i番目の角度を取得
    h = h_n[i]
    
    # なす角を作図
    ax.quiver(0, 0, 0, *a, 
              color='red', arrow_length_ratio=l/np.linalg.norm(a), 
              label='$a=('+', '.join(map(str, a))+')$') # ベクトルa
    ax.quiver(0, 0, 0, *b, 
              color='blue', arrow_length_ratio=l/np.linalg.norm(b), 
              label='$b=('+', '.join(map(str, b.round(2)))+')$') # ベクトルb
    ax.plot(*X.T, 
            color='black', linewidth=1, zorder=10) # 角マーク
    ax.text(*x, 
            s='$\\theta$', size=15, ha='center', va='center', zorder=15) # 角ラベル
    ax.quiver([0, a[0], b[0]], 
              [0, a[1], b[1]], 
              [0, a[2], b[2]], 
              [0, 0, 0], 
              [0, 0, 0], 
              [z_min, z_min-a[2], z_min-b[2]], 
              color=['gray', 'red', 'blue'], 
              arrow_length_ratio=0, linestyle=':') # 座標用の補助線
    ax.set_xlim(left=x_min, right=x_max)
    ax.set_ylim(bottom=y_min, top=y_max)
    ax.set_zlim(bottom=z_min, top=z_max)
    ax.set_xlabel('$x_1$')
    ax.set_ylabel('$x_2$')
    ax.set_zlabel('$x_3$')
    ax.set_title('$\\theta=' + str((theta/np.pi).round(2)) + '\pi, ' + 
                 '\\theta^{\circ}=' + str((theta/np.pi*180.0).round(1)) + '^{\circ}$', loc='left')
    ax.legend()
    ax.set_aspect('equal')
    ax.view_init(elev=30, azim=h) # 表示角度

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

# gif画像を保存
ani.save(filename='angle360_3d.gif')

 作図処理をupdate()として定義して、FuncAnimation()でgif画像を作成する。作成したグラフオブジェクトのsave()メソッドでgifファイルとして書き出す。

ベクトルのなす角のグラフ

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

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

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

# 固定するベクトルを指定
a = np.array([1.0, 1.0, 1.0])

# 変化するベクトル用のラジアンを作成
t_n = np.linspace(start=0.0, stop=2.0*np.pi, num=frame_num+1)[:frame_num]
u_n = np.linspace(start=2.0/3.0*np.pi, stop=5.0/3.0*np.pi, num=frame_num+1)[:frame_num]

# 変化するベクトルのノルムを指定
d = 2.0

# 変化するベクトルを作成
b_n = np.array(
    [d * np.sin(t_n)*np.cos(u_n), d * np.sin(t_n)*np.sin(u_n), d * np.cos(t_n)]
).T

# 矢のサイズを指定
l = 0.2

# グラフサイズ用の値を設定
axis_size = d
z_min = -axis_size

# グラフオブジェクトを初期化
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
fig.suptitle('$\\theta = \\arccos(\\frac{a \cdot b}{\|a\| \|b\|})$', fontsize=20)

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

    # i番目のベクトルを取得
    b = b_n[i]
    
    # ベクトルa,bのなす角を計算
    theta = np.arccos(np.dot(a, b) / np.linalg.norm(a) / np.linalg.norm(b))
    
    # 角マークの座標を計算
    r = 0.3
    X = Beta1 * a/np.linalg.norm(a) + Beta2 * b/np.linalg.norm(b) # 2点間の曲線
    X *= r / np.linalg.norm(X, axis=1, keepdims=True) # サイズ調整
    
    # 角ラベルの座標を計算
    r = 0.5
    x = beta1 * a/np.linalg.norm(a) + beta2 * b/np.linalg.norm(b) # 2点間の中点
    x *= r / np.linalg.norm(x) # サイズ調整
    
    # なす角を作図
    ax.quiver(0, 0, 0, *a, 
              color='red', arrow_length_ratio=l/np.linalg.norm(a), 
              label='$a=('+', '.join(map(str, a))+')$') # ベクトルa
    ax.quiver(0, 0, 0, *b, 
              color='blue', arrow_length_ratio=l/np.linalg.norm(b), 
              label='$b=('+', '.join(map(str, b.round(2)))+')$') # ベクトルb
    ax.plot(*X.T, 
            color='black', linewidth=1, zorder=10) # 角マーク
    ax.text(*x, 
            s='$\\theta$', size=15, ha='center', va='center', zorder=15) # 角ラベル
    ax.quiver([0, a[0], b[0]], 
              [0, a[1], b[1]], 
              [0, a[2], b[2]], 
              [0, 0, 0], 
              [0, 0, 0], 
              [z_min, z_min-a[2], z_min-b[2]], 
              color=['gray', 'red', 'blue'], 
              arrow_length_ratio=0, linestyle=':') # 座標用の補助線
    ax.set_xlim(left=-axis_size, right=axis_size)
    ax.set_ylim(bottom=-axis_size, top=axis_size)
    ax.set_zlim(bottom=-axis_size, top=axis_size)
    ax.set_xlabel('$x_1$')
    ax.set_ylabel('$x_2$')
    ax.set_zlabel('$x_3$')
    ax.set_title('$\\theta=' + str((theta/np.pi).round(2)) + '\pi, ' + 
                 '\\theta^{\circ}=' + str((theta/np.pi*180.0).round(1)) + '^{\circ}$', 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(filename='angle_3d.gif')

ベクトルのなす角のグラフ


 以上で、なす角を示すグラフを作成できた。次からは、(どの説明が分かりやすいのか分からないので)3つの方法でこの処理(計算)を確認する。

円弧の描画

 2つの3次元ベクトルを繋ぐ円弧(角マーク)の座標計算とノルム操作をグラフで確認する。

 3次元ベクトル \mathbf{a}, \mathbf{b}を作成する。

# 3次元ベクトルを指定
a = np.array([1.0, 1.0, 1.0])
b = np.array([-1.0, -1.0, 1.0])

 2つのベクトル \mathbf{a} = (a_1, a_2, a_3), \mathbf{b} = (b_1, b_2, b_3)a, bとして値を指定する。

 ベクトル \mathbf{a}, \mathbf{b}のグラフを作成する。 

# グラフサイズ用の値を設定
x_min = np.floor(np.min([0.0, a[0], b[0]])) - 1.0
x_max =  np.ceil(np.max([0.0, a[0], b[0]])) + 1.0
y_min = np.floor(np.min([0.0, a[1], b[1]])) - 1.0
y_max =  np.ceil(np.max([0.0, a[1], b[1]])) + 1.0
z_min = np.floor(np.min([0.0, a[2], b[2]])) - 1.0
z_max =  np.ceil(np.max([0.0, a[2], b[2]])) + 1.0

# 矢のサイズを指定
l = 0.2

# ベクトルを作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
ax.scatter(0, 0, 0, 
           color='orange', s=50, label='$O=(0, 0, 0)$') # 原点
ax.quiver(0, 0, 0, *a, 
          color='red', arrow_length_ratio=l/np.linalg.norm(a), 
           label='$a=('+', '.join(map(str, a))+')$') # ベクトルa
ax.quiver(0, 0, 0, *b, 
          color='blue', arrow_length_ratio=l/np.linalg.norm(b), 
           label='$b=('+', '.join(map(str, b))+')$') # ベクトルb
ax.quiver([0, a[0], b[0]], 
          [0, a[1], b[1]], 
          [0, a[2], b[2]], 
          [0, 0, 0], 
          [0, 0, 0], 
          [z_min, z_min-a[2], z_min-b[2]], 
          color=['orange', 'red', 'blue'], 
          arrow_length_ratio=0, linestyle=':') # 座標用の補助線
ax.set_xticks(ticks=np.arange(x_min, x_max+0.1, 0.5))
ax.set_yticks(ticks=np.arange(y_min, y_max+0.1, 0.5))
ax.set_zticks(ticks=np.arange(z_min, z_max+0.1, 0.5))
ax.set_zlim(bottom=z_min, top=z_max)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('$\|a\|=' + str(np.linalg.norm(a).round(2)) + ', ' + 
             '\|b\|=' + str(np.linalg.norm(b).round(2)) + '$', loc='left')
fig.suptitle('vectors', fontsize=20)
ax.legend()
ax.set_aspect('equal')
plt.show()

ベクトル

 axes.quiver()の第1・2・3引数に始点の座標、第4・5・6引数にベクトルの移動量を指定する。
 各配列の形状に応じて、配列の前に*を付けてアンパック(展開)して座標などを指定している。原点を始点とする。

 原点からのノルム rを指定して、ベクトル上の点(ノルムが rのベクトル) \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を計算する。

# 円弧の半径を指定
r = 0.5

# ベクトルのノルムを調整
tilde_a = r * a / np.linalg.norm(a)
tilde_b = r * b / np.linalg.norm(b)
print(tilde_a)
print(np.linalg.norm(tilde_a))
print(tilde_b)
print(np.linalg.norm(tilde_b))
[0.28867513 0.28867513 0.28867513]
0.5
[-0.28867513 -0.28867513  0.28867513]
0.5

  \mathbf{a}, \mathbf{b}をそれぞれ正規化し正の実数 rを掛けて、ノルムを rに変更したベクトル \tilde{\mathbf{a}} = (\tilde{a}_1, \tilde{a}_2, \tilde{a}_3), \tilde{\mathbf{b}} = (\tilde{b}_1, \tilde{b}_2, \tilde{b}_3)tilde_a, tilde_bとする。

 \displaystyle
\begin{aligned}
\tilde{\mathbf{a}}
   &= r \frac{\mathbf{a}}{\|\mathbf{a}\|}
\\
\tilde{\mathbf{b}}
   &= r \frac{\mathbf{b}}{\|\mathbf{b}\|}
\end{aligned}

 ノルムの操作については「なす角の描画」を参照のこと。

 ベクトル \mathbf{a}, \mathbf{b}上の点 \tilde{\mathbf{a}}, \tilde{\mathbf{b}}のグラフを作成する。 

# ベクトル上の点を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
ax.scatter(0, 0, 0, 
           color='orange', s=50, label='$O=(0, 0, 0)$') # 原点
ax.scatter(*tilde_a, 
           color='red', s=50, 
           label='$\\tilde{a}=('+', '.join(map(str, tilde_a.round(2)))+')$') # ベクトルa上の点
ax.scatter(*tilde_b, 
           color='blue', s=50, 
           label='$\\tilde{b}=('+', '.join(map(str, tilde_b.round(2)))+')$') # ベクトルb上の点
ax.quiver(0, 0, 0, *a, 
          color='red', arrow_length_ratio=l/np.linalg.norm(a)) # ベクトルa
ax.quiver(0, 0, 0, *b, 
          color='blue', arrow_length_ratio=l/np.linalg.norm(b)) # ベクトルb
ax.quiver([0, a[0], b[0], tilde_a[0], tilde_b[0]], 
          [0, a[1], b[1], tilde_a[1], tilde_b[1]], 
          [0, a[2], b[2], tilde_a[2], tilde_b[2]], 
          [0, 0, 0, 0, 0], 
          [0, 0, 0, 0, 0], 
          [z_min, z_min-a[2], z_min-b[2], z_min-tilde_a[2], z_min-tilde_b[2]], 
          color=['orange', 'red', 'blue', 'red', 'blue'], 
          arrow_length_ratio=0, linestyle=':') # 座標用の補助線
ax.set_xticks(ticks=np.arange(x_min, x_max+0.1, 0.5))
ax.set_yticks(ticks=np.arange(y_min, y_max+0.1, 0.5))
ax.set_zticks(ticks=np.arange(z_min, z_max+0.1, 0.5))
ax.set_zlim(bottom=z_min, top=z_max)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('$r=' + str(r) + ', ' + 
             '\|\\tilde{a}\|=' + str(np.linalg.norm(tilde_a).round(2)) + ', ' + 
             '\|\\tilde{b}\|=' + str(np.linalg.norm(tilde_b).round(2)) + '$', loc='left')
fig.suptitle('$\\tilde{a} = r \\frac{a}{\|a\|}, ' + 
             '\\tilde{b} = r \\frac{b}{\|b\|}$', fontsize=20)
ax.legend()
ax.set_aspect('equal')
plt.show()

ベクトル上の点

 点 \tilde{\mathbf{a}}, \tilde{\mathbf{b}}がそれぞれ、ベクトル \mathbf{a}, \mathbf{b}と平行な直線上で、原点からのノルムが rの点になるのを確認できる。

 ベクトル \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を線形結合して、点 \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を結ぶ曲線を計算する。

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

# 曲線用の係数を作成
beta1 = np.cos(t)
beta2 = np.sin(t)
print(beta1[:5])
print(beta2[:5])

# 2点間の曲線の座標を計算
x = beta1 * tilde_a[0] + beta2 * tilde_b[0]
y = beta1 * tilde_a[1] + beta2 * tilde_b[1]
z = beta1 * tilde_a[2] + beta2 * tilde_b[2]
print(x[:5])
print(y[:5])
print(z[:5])
[0.         0.01586663 0.03173326 0.04759989 0.06346652]
[1.         0.99987413 0.99949654 0.99886734 0.99798668]
[0.         0.01586596 0.03172793 0.04758192 0.06342392]
[0.28867513 0.28405869 0.27937073 0.27461245 0.26978503]
[0.28867513 0.28405869 0.27937073 0.27461245 0.26978503]
[0.28867513 0.29321891 0.29768886 0.30208388 0.30640285]

  \cos 0 = 1, \sin 0 = 0 \cos \frac{\pi}{2} = 0, \sin \frac{\pi}{2} = 1であるのを利用して、線形結合の係数を \beta_1 = \cos t, \beta_2 = \sin tとする。 0 \leq t \leq \frac{\pi}{2}の範囲のラジアンを用いると、係数は 1 \geq \beta_1 \geq 0, 0 \leq \beta_2 \leq 1の値をとる。
 ベクトル \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を線形結合したベクトルを \mathbf{x} = (x_1, x_2, x_3)とする。

 \displaystyle
\mathbf{x}
    = \beta_1 \tilde{\mathbf{a}}
      + \beta_2 \tilde{\mathbf{b}}

  x_1, x_2, x_3をそれぞれ計算してx, y, zとする。x, y, zの各インデックスが1つの点 \mathbf{x}に対応する。

  \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を結ぶ曲線のグラフを作成する。

# 2点間の曲線を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
ax.scatter(0, 0, 0, 
           color='orange', s=50, label='$O=(0, 0, 0)$') # 原点
ax.scatter(*tilde_a, 
           color='red', s=50, 
           label='$\\tilde{a}=('+', '.join(map(str, tilde_a.round(2)))+')$') # ベクトルa上の点
ax.scatter(*tilde_b, 
           color='blue', s=50, 
           label='$\\tilde{b}=('+', '.join(map(str, tilde_b.round(2)))+')$') # ベクトルb上の点
ax.quiver(0, 0, 0, *a, 
          color='red', arrow_length_ratio=l/np.linalg.norm(a)) # ベクトルa
ax.quiver(0, 0, 0, *b, 
          color='blue', arrow_length_ratio=l/np.linalg.norm(b)) # ベクトルb
ax.plot(x, y, z, 
        color='black', zorder=10) # 2点間の曲線
ax.quiver([0, a[0], b[0], tilde_a[0], tilde_b[0]], 
          [0, a[1], b[1], tilde_a[1], tilde_b[1]], 
          [0, a[2], b[2], tilde_a[2], tilde_b[2]], 
          [0, 0, 0, 0, 0], 
          [0, 0, 0, 0, 0], 
          [z_min, z_min-a[2], z_min-b[2], z_min-tilde_a[2], z_min-tilde_b[2]], 
          color=['orange', 'red', 'blue', 'red', 'blue'], 
          arrow_length_ratio=0, linestyle=':') # 座標用の補助線
ax.set_xticks(ticks=np.arange(x_min, x_max+0.1, 0.5))
ax.set_yticks(ticks=np.arange(y_min, y_max+0.1, 0.5))
ax.set_zticks(ticks=np.arange(z_min, z_max+0.1, 0.5))
ax.set_zlim(bottom=z_min, top=z_max)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('$0 \leq t \leq \\frac{\pi}{2}$', loc='left')
fig.suptitle('$x = \cos t\ \\tilde{a} + \sin t\ \\tilde{b}$', fontsize=20)
ax.legend()
ax.set_aspect('equal')
plt.show()

ベクトル上の点を結ぶ曲線

 係数 \beta_1, \beta_2を変化させた線形結合の点 \mathbf{x}を繋いだ曲線が、点 \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を結ぶ曲線になるのを確認できる。

 曲線のノルムを rに変更する。

# 曲線の座標を結合
xyz = np.stack([x, y, z], axis=1)
print(xyz[:5])
print(xyz.shape)

# 点ごとにノルムを計算
norm_xyz = np.linalg.norm(xyz, axis=1)
print(norm_xyz[:5])
print(norm_xyz.shape)

# 曲線のノルムを調整
tilde_x = r * x / norm_xyz
tilde_y = r * y / norm_xyz
tilde_z = r * z / norm_xyz
print(tilde_x[:5])
print(tilde_y[:5])
print(tilde_z[:5])
print(np.linalg.norm(np.stack([tilde_x, tilde_y, tilde_z], axis=1), axis=1)[:5])
[[0.28867513 0.28867513 0.28867513]
 [0.28405869 0.28405869 0.29321891]
 [0.27937073 0.27937073 0.29768886]
 [0.27461245 0.27461245 0.30208388]
 [0.26978503 0.26978503 0.30640285]]
(100, 3)
[0.5        0.49734898 0.49468644 0.4920149  0.48933693]
(100,)
[0.28867513 0.28557281 0.28237153 0.27906924 0.27566388]
[0.28867513 0.28557281 0.28237153 0.27906924 0.27566388]
[0.28867513 0.29478185 0.30088642 0.30698651 0.31307963]
[0.5 0.5 0.5 0.5 0.5]

  \mathbf{x}を正規化し正の実数 rを掛けて、ノルムを rに変更したベクトルを \tilde{\mathbf{x}} = (\tilde{x}_1, \tilde{x}_2, \tilde{x}_3)とする。

 \displaystyle
\tilde{\mathbf{x}}
    = r \frac{\mathbf{x}}{\|\mathbf{x}\|}

  \tilde{x}_1, \tilde{x}_2, \tilde{x}_3をそれぞれ計算してtilde_x, tilde_y, tilde_zとする。tilde_x, tilde_y, tilde_zの各インデックスが1つの点 \tilde{\mathbf{x}}に対応する。

  \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を結ぶ円弧(原点からのノルムが rの曲線)のグラフを作成する。

# 2点間の円弧を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
ax.scatter(0, 0, 0, 
           color='orange', s=50, label='$O=(0, 0, 0)$') # 原点
ax.scatter(*tilde_a, 
           color='red', s=50, 
           label='$\\tilde{a}=('+', '.join(map(str, tilde_a.round(2)))+')$') # ベクトルa上の点
ax.scatter(*tilde_b, 
           color='blue', s=50, 
           label='$\\tilde{b}=('+', '.join(map(str, tilde_b.round(2)))+')$') # ベクトルb上の点
ax.quiver(0, 0, 0, *a, 
          color='red', arrow_length_ratio=l/np.linalg.norm(a)) # ベクトルa
ax.quiver(0, 0, 0, *b, 
          color='blue', arrow_length_ratio=l/np.linalg.norm(b)) # ベクトルb
ax.plot(tilde_x, tilde_y, tilde_z, 
        color='black', zorder=10) # 2点間の円弧
ax.quiver([0, a[0], b[0], tilde_a[0], tilde_b[0]], 
          [0, a[1], b[1], tilde_a[1], tilde_b[1]], 
          [0, a[2], b[2], tilde_a[2], tilde_b[2]], 
          [0, 0, 0, 0, 0], 
          [0, 0, 0, 0, 0], 
          [z_min, z_min-a[2], z_min-b[2], z_min-tilde_a[2], z_min-tilde_b[2]], 
          color=['orange', 'red', 'blue', 'red', 'blue'], 
          arrow_length_ratio=0, linestyle=':') # 座標用の補助線
ax.set_xticks(ticks=np.arange(x_min, x_max+0.1, 0.5))
ax.set_yticks(ticks=np.arange(y_min, y_max+0.1, 0.5))
ax.set_zticks(ticks=np.arange(z_min, z_max+0.1, 0.5))
ax.set_zlim(bottom=z_min, top=z_max)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('$r=' + str(r) + '$', loc='left')
fig.suptitle('$\\tilde{x} = r \\frac{x}{\|x\|}$', fontsize=20)
ax.legend()
ax.set_aspect('equal')
#ax.view_init(elev=90, azim=270) # xy面
#ax.view_init(elev=0, azim=270) # xz面
#ax.view_init(elev=0, azim=0) # yz面
plt.show()

ベクトル上の点を結ぶ円弧

 ベクトル \mathbf{a}, \mathbf{b}を結ぶいい感じの曲線が得られた。

 グラフを回転させて確認する。

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

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

# 水平方向の角度として利用する値を作成
h_n = np.linspace(start=0.0, stop=360.0, num=frame_num+1)[:frame_num]

# グラフオブジェクトを初期化
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
fig.suptitle('circular arc', fontsize=20)

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

    # i番目の角度を取得
    h = h_n[i]
    
    # 2点間の円弧を作図
    ax.scatter(0, 0, 0, 
               color='orange', s=50, label='$O=(0, 0, 0)$') # 原点
    ax.scatter(*tilde_a, 
               color='red', s=50, 
               label='$\\tilde{a}=('+', '.join(map(str, tilde_a.round(2)))+')$') # ベクトルa上の点
    ax.scatter(*tilde_b, 
               color='blue', s=50, 
               label='$\\tilde{b}=('+', '.join(map(str, tilde_b.round(2)))+')$') # ベクトルb上の点
    ax.quiver(0, 0, 0, *a, 
              color='red', arrow_length_ratio=l/np.linalg.norm(a)) # ベクトルa
    ax.quiver(0, 0, 0, *b, 
              color='blue', arrow_length_ratio=l/np.linalg.norm(b)) # ベクトルb
    ax.plot(x, y, z, 
            color='purple', zorder=10) # 2点間の曲線
    ax.plot(tilde_x, tilde_y, tilde_z, 
            color='black', zorder=10) # 2点間の円弧
    ax.quiver([0, a[0], b[0], tilde_a[0], tilde_b[0]], 
              [0, a[1], b[1], tilde_a[1], tilde_b[1]], 
              [0, a[2], b[2], tilde_a[2], tilde_b[2]], 
              [0, 0, 0, 0, 0], 
              [0, 0, 0, 0, 0], 
              [z_min, z_min-a[2], z_min-b[2], z_min-tilde_a[2], z_min-tilde_b[2]], 
              color=['orange', 'red', 'blue', 'red', 'blue'], 
              arrow_length_ratio=0, linestyle=':') # 座標用の補助線
    ax.set_xticks(ticks=np.arange(x_min, x_max+0.1, 0.5))
    ax.set_yticks(ticks=np.arange(y_min, y_max+0.1, 0.5))
    ax.set_zticks(ticks=np.arange(z_min, z_max+0.1, 0.5))
    ax.set_zlim(bottom=z_min, top=z_max)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_title('$r=' + str(r) + '$', loc='left')
    ax.legend()
    ax.set_aspect('equal')
    ax.view_init(elev=30, azim=h) # 表示角度

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

# gif画像を保存
ani.save(filename='arc360_3d.gif')

ベクトル上の点を結ぶ曲線

 正規化したベクトル \tilde{\mathbf{a}}, \tilde{\mathbf{b}}の線形結合 \mathbf{x}を紫色の曲線、ノルムを調整した円弧 \tilde{\mathbf{x}}を黒色の曲線で描画する。

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

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

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

# 円弧の半径を指定
r = 0.5

# 固定するベクトルを指定
a = np.array([1.0, 1.0, 1.0])
tilde_a = r * a / np.linalg.norm(a) # ノルムを調整

# 変化するベクトル用のラジアンを作成
t_n = np.linspace(start=-1.0*np.pi, stop=1.0*np.pi, num=frame_num+1)[:frame_num]
u_n = np.linspace(start=2.0/3.0*np.pi, stop=5.0/3.0*np.pi, num=frame_num+1)[:frame_num]

# 変化するベクトルのサイズを指定
d = 2.0

# 変化するベクトルを作成
b_n = np.array(
    [d * np.sin(t_n)*np.cos(u_n), d * np.sin(t_n)*np.sin(u_n), d * np.cos(t_n)]
).T

# 矢のサイズを指定
l = 0.2

# グラフサイズ用の値を設定
axis_size = np.ceil(d)
z_min = -axis_size

# グラフオブジェクトを初期化
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
fig.suptitle('circular arc', fontsize=20)

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

    # i番目のベクトルを取得
    b = b_n[i]
    tilde_b = r * b / np.linalg.norm(b) # ノルムを調整
    
    # ベクトルa,bのなす角を計算
    theta = np.arccos(np.dot(a, b) / np.linalg.norm(a) / np.linalg.norm(b))
    
    # 2点間の曲線の座標を計算
    x = beta1 * tilde_a[0] + beta2 * tilde_b[0]
    y = beta1 * tilde_a[1] + beta2 * tilde_b[1]
    z = beta1 * tilde_a[2] + beta2 * tilde_b[2]
    
    # 点ごとにノルムを計算
    norm_xyz = np.linalg.norm(np.stack([x, y, z], axis=1), axis=1)
    
    # 曲線のノルムを調整
    tilde_x = r * x / norm_xyz
    tilde_y = r * y / norm_xyz
    tilde_z = r * z / norm_xyz
    
    # 2点間の円弧を作図
    ax.scatter(0, 0, 0, 
               color='orange', s=50, label='$O=(0, 0, 0)$') # 原点
    ax.scatter(*tilde_a, 
               color='red', s=50, 
               label='$\\tilde{a}=('+', '.join(map(str, tilde_a.round(2)))+')$') # ベクトルa上の点
    ax.scatter(*tilde_b, 
               color='blue', s=50, 
               label='$\\tilde{b}=('+', '.join(map(str, tilde_b.round(2)))+')$') # ベクトルb上の点
    ax.quiver(0, 0, 0, *a, 
              color='red', arrow_length_ratio=l/np.linalg.norm(a)) # ベクトルa
    ax.quiver(0, 0, 0, *b, 
              color='blue', arrow_length_ratio=l/np.linalg.norm(b)) # ベクトルb
    ax.quiver(0, 0, z_min, a[0], a[1], 0, 
              color='red', linestyle='--', arrow_length_ratio=l/np.linalg.norm(a)) # 底面上のベクトルa
    ax.quiver(0, 0, z_min, b[0], b[1], 0, 
              color='blue', linestyle='--', arrow_length_ratio=l/np.linalg.norm(b)) # 底面上のベクトルb
    ax.plot(x, y, z, 
            color='purple', zorder=10) # 2点間の曲線
    ax.plot(tilde_x, tilde_y, tilde_z, 
            color='black', zorder=15) # 2点間の円弧
    ax.quiver([0, a[0], b[0], tilde_a[0], tilde_b[0]], 
              [0, a[1], b[1], tilde_a[1], tilde_b[1]], 
              [0, a[2], b[2], tilde_a[2], tilde_b[2]], 
              [0, 0, 0, 0, 0], 
              [0, 0, 0, 0, 0], 
              [z_min, z_min-a[2], z_min-b[2], z_min-tilde_a[2], z_min-tilde_b[2]], 
              color=['orange', 'red', 'blue', 'red', 'blue'], 
              arrow_length_ratio=0, linestyle=':') # 座標用の補助線
    ax.set_xticks(ticks=np.arange(-axis_size, axis_size+0.1, 0.5))
    ax.set_yticks(ticks=np.arange(-axis_size, axis_size+0.1, 0.5))
    ax.set_zticks(ticks=np.arange(-axis_size, axis_size+0.1, 0.5))
    ax.set_xlim(left=-axis_size, right=axis_size)
    ax.set_ylim(bottom=-axis_size, top=axis_size)
    ax.set_zlim(bottom=-axis_size, top=axis_size)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_title('$r=' + str(r) + ', ' + 
                 '\\theta=' + str((theta/np.pi).round(2)) + '\pi$', 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(filename='arc_3d.gif')

2ベクトル上の各点を結ぶ曲線

 なす角が \theta \lt \frac{\pi}{2}だと \|\mathbf{x}\| \geq r(紫色の曲線が黒色の曲線の外側)になり、 \theta \gt \frac{\pi}{2}だと \|\mathbf{x}\| \leq r(紫色の曲線が黒色の曲線の内側)になる模様?(理由が気になるけど今回は保留)。

 ここまでは、曲線の計算を確認した。次からは、2つの曲線の関係を確認する。

球面上で可視化

 「円弧の描画」では、一定のノルムで2つのベクトルを結ぶ曲線の計算を確認した。
 ここでは、その曲線が球面上の2点を結ぶ曲線になることをグラフで確認する。

 半径 rを指定して、球面の座標を計算する。

# 球面の半径を指定
r = 1.2

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

# ラジアンの格子点を作成
T, U = np.meshgrid(t, t)
print(T.shape)
print(U.shape)

# 球面の座標を計算
X = r * np.sin(T) * np.cos(U)
Y = r * np.sin(T) * np.sin(U)
Z = r * np.cos(T)
print(X[:5, :5])
print(Y[:5, :5])
print(Z[:5, :5])
[0.         0.06346652 0.12693304 0.19039955 0.25386607]
(100,)
(100, 100)
(100, 100)
[[0.         0.0761087  0.15191094 0.22710149 0.30137758]
 [0.         0.07595547 0.1516051  0.22664426 0.30077081]
 [0.         0.07549639 0.15068879 0.22527442 0.29895295]
 [0.         0.07473332 0.14916572 0.22299747 0.2959313 ]
 [0.         0.07366932 0.147042   0.2198226  0.29171804]]
[[0.         0.         0.         0.         0.        ]
 [0.         0.00482711 0.00963479 0.01440367 0.01911455]
 [0.         0.00963479 0.01923078 0.02874934 0.03815213]
 [0.         0.01440367 0.02874934 0.04297924 0.05703608]
 [0.         0.01911455 0.03815213 0.05703608 0.07569037]]
[[1.2        1.19758401 1.19034578 1.17831444 1.16153844]
 [1.2        1.19758401 1.19034578 1.17831444 1.16153844]
 [1.2        1.19758401 1.19034578 1.17831444 1.16153844]
 [1.2        1.19758401 1.19034578 1.17831444 1.16153844]
 [1.2        1.19758401 1.19034578 1.17831444 1.16153844]]

 詳しくは「球面の作図【Matplotlib】 - からっぽのしょこ」を参照のこと。

 3次元ベクトル \mathbf{a}, \mathbf{b}を作成して、半径が rの球面上の点 \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を計算する。

# 3次元ベクトルを指定
a = np.array([1.0, 1.0, 1.0])
b = np.array([-1.0, -1.0, 1.0])

# 球面上の点の座標を計算
tilde_a = r * a / np.linalg.norm(a)
tilde_b = r * b / np.linalg.norm(b)
print(tilde_a)
print(tilde_b)
[0.69282032 0.69282032 0.69282032]
[-0.69282032 -0.69282032  0.69282032]

 「円弧の描画」のときと同様にして、ベクトル \mathbf{a} = (a_1, a_2, a_3), \mathbf{b} = (b_1, b_2, b_3)a, bとして値を指定して、ノルムが rのベクトル \tilde{\mathbf{a}} = (\tilde{a}_1, \tilde{a}_2, \tilde{a}_3), \tilde{\mathbf{b}} = (\tilde{b}_1, \tilde{b}_2, \tilde{b}_3)を計算してtilde_a, tilde_bとする。

 球面上の点 \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を結ぶ曲線を計算する。

# 曲線用のラジアンを作成
#t = np.linspace(start=0.0, stop=0.5*np.pi, num=100) # なす角用
t = np.linspace(start=0.0, stop=2.0*np.pi, num=300) # 1周用

# 曲線用の係数を作成
beta1 = np.cos(t)
beta2 = np.sin(t)

# 2点間の曲線の座標を計算
x = beta1 * tilde_a[0] + beta2 * tilde_b[0]
y = beta1 * tilde_a[1] + beta2 * tilde_b[1]
z = beta1 * tilde_a[2] + beta2 * tilde_b[2]
print(x[:5])
print(y[:5])
print(z[:5])

# 点ごとにノルムを計算
xyz = np.stack([x, y, z], axis=1) # 座標を結合
norm_xyz = np.linalg.norm(xyz, axis=1)

# 球面上の曲線の座標を計算
tilde_x = r * x / norm_xyz
tilde_y = r * y / norm_xyz
tilde_z = r * z / norm_xyz
print(tilde_x[:5])
print(tilde_y[:5])
print(tilde_z[:5])
[0.69282032 0.6781095  0.66309925 0.64779619 0.63220709]
[0.69282032 0.6781095  0.66309925 0.64779619 0.63220709]
[0.69282032 0.70722521 0.72131781 0.73509189 0.74854138]
[0.69282032 0.6829085  0.6725773  0.66181524 0.65061089]
[0.69282032 0.6829085  0.6725773  0.66181524 0.65061089]
[0.69282032 0.71223026 0.73162801 0.75100012 0.77033172]

 「円弧の描画」のときと同様にして、曲線の座標を計算する。
 係数用のラジアンtに関して、範囲を 0 \leq t \leq \frac{\pi}{2}とすると、2点のなす角 \thetaを結ぶ曲線を描画できた。 0 \leq t \leq 2 \piとすると、なす角の反対側の角 2 \pi - \theta(度数法だと 360^{\circ} - \theta^{\circ})も含めた曲線(円)になる。

 球面上の点 \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を結ぶ曲線のグラフを作成する。

# グラフサイズ用の値を設定
axis_size = np.ceil(r)# + 1.0
z_min = -axis_size

# 矢のサイズを指定
l = 0.2

# 球面上の2点を結ぶ曲線を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
ax.plot_wireframe(X, Y, Z, alpha=0.1) # 球面
ax.quiver(0, 0, 0, *a, 
          color='red', arrow_length_ratio=l/np.linalg.norm(a)) # ベクトルa
ax.quiver(0, 0, 0, *b, 
          color='blue', arrow_length_ratio=l/np.linalg.norm(b)) # ベクトルb
ax.scatter(*tilde_a, 
           color='red', s=50, 
           label='$\\tilde{a}=('+', '.join(map(str, tilde_a.round(2)))+')$') # ベクトルaの球面上の点
ax.scatter(*tilde_b, 
           color='blue', s=50, 
           label='$\\tilde{b}=('+', '.join(map(str, tilde_b.round(2)))+')$') # ベクトルbの球面上の点
ax.plot(x, y, z, 
        color='purple', zorder=0) # 2点間の曲線
ax.plot(tilde_x, tilde_y, tilde_z, 
        color='black', zorder=0) # 球面上の2点間の曲線
ax.quiver([0, tilde_a[0], tilde_b[0]], 
          [0, tilde_a[1], tilde_b[1]], 
          [0, tilde_a[2], tilde_b[2]], 
          [0, 0, 0], 
          [0, 0, 0], 
          [z_min, z_min-tilde_a[2], z_min-tilde_b[2]], 
          color=['gray', 'red', 'blue'], 
          arrow_length_ratio=0, linestyle=':') # 座標用の補助線
ax.set_xticks(ticks=np.arange(-axis_size, axis_size+0.1, 0.5))
ax.set_yticks(ticks=np.arange(-axis_size, axis_size+0.1, 0.5))
ax.set_zticks(ticks=np.arange(-axis_size, axis_size+0.1, 0.5))
ax.set_zlim(bottom=-axis_size, top=axis_size)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('$r=' + str(r) + '$', loc='left')
fig.suptitle('curve on the sphere', fontsize=20)
ax.legend()
ax.set_aspect('equal')
#ax.view_init(elev=90, azim=270) # xy面
#ax.view_init(elev=0, azim=270) # xz面
#ax.view_init(elev=0, azim=0) # yz面
plt.show()

球面上の2点を結ぶ曲線

 (各グラフを2次元的に重ねて描画するので分かりにくいかも知れない。特に、点が球の表側なのか裏側なのかを判別できないと見誤る。紫色の曲線が傾いた円に見えるかもしれないが、楕円(正確には楕円ではないのかもしれないが円ではない)である。表示角度を調整すると、紫色の曲線と黒色の曲線が同一平面上なのが分かる。)

 グラフを回転させて確認する。

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

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

# 水平方向の角度として利用する値を作成
h_n = np.linspace(start=0.0, stop=360.0, num=frame_num+1)[:frame_num]

# グラフオブジェクトを初期化
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
fig.suptitle('curve on the sphere', fontsize=20)

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

    # i番目の角度を取得
    h = h_n[i]
    
    # 球面上の2点を結ぶ曲線を作図
    ax.plot_wireframe(X, Y, Z, alpha=0.1) # 球面
    ax.quiver(0, 0, 0, *a, 
              color='red', arrow_length_ratio=l/np.linalg.norm(a)) # ベクトルa
    ax.quiver(0, 0, 0, *b, 
              color='blue', arrow_length_ratio=l/np.linalg.norm(b)) # ベクトルb
    ax.scatter(*tilde_a, 
               color='red', s=50, 
               label='$\\tilde{a}=('+', '.join(map(str, tilde_a.round(2)))+')$') # ベクトルaの球面上の点
    ax.scatter(*tilde_b, 
               color='blue', s=50, 
               label='$\\tilde{b}=('+', '.join(map(str, tilde_b.round(2)))+')$') # ベクトルbの球面上の点
    ax.plot(x, y, z, 
            color='purple', zorder=0) # 2点間の曲線
    ax.plot(tilde_x, tilde_y, tilde_z, 
            color='black', zorder=0) # 球面上の2点間の曲線
    ax.quiver([0, tilde_a[0], tilde_b[0]], 
              [0, tilde_a[1], tilde_b[1]], 
              [0, tilde_a[2], tilde_b[2]], 
              [0, 0, 0], 
              [0, 0, 0], 
              [z_min, z_min-tilde_a[2], z_min-tilde_b[2]], 
              color=['gray', 'red', 'blue'], 
              arrow_length_ratio=0, linestyle=':') # 座標用の補助線
    ax.set_xticks(ticks=np.arange(-axis_size, axis_size+0.1, 0.5))
    ax.set_yticks(ticks=np.arange(-axis_size, axis_size+0.1, 0.5))
    ax.set_zticks(ticks=np.arange(-axis_size, axis_size+0.1, 0.5))
    ax.set_zlim(bottom=-axis_size, top=axis_size)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_title('$r=' + str(r) + '$', loc='left')
    ax.legend()
    ax.set_aspect('equal')
    ax.view_init(elev=30, azim=h) # 表示角度

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

# gif画像を保存
ani.save(filename='CurveOnSphere360_3d.gif')

球面上の2点を結ぶ曲線

 (ところで、黒色の曲線は球面上の直線?)

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

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

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

# 固定するベクトルを指定
a = np.array([1.0, 1.0, 1.0])
tilde_a = r * a / np.linalg.norm(a) # ノルムを調整

# 変化するベクトル用のラジアンを作成
t_n = np.linspace(start=-1.0*np.pi, stop=1.0*np.pi, num=frame_num+1)[:frame_num]
u_n = np.linspace(start=2.0/3.0*np.pi, stop=5.0/3.0*np.pi, num=frame_num+1)[:frame_num]

# 変化するベクトルのサイズを指定
d = 2.0

# 変化するベクトルを作成
b_n = np.array(
    [d * np.sin(t_n)*np.cos(u_n), d * np.sin(t_n)*np.sin(u_n), d * np.cos(t_n)]
).T

# 矢のサイズを指定
l = 0.2

# グラフサイズ用の値を設定
axis_size = np.ceil(np.max([r, d]))# + 1.0
z_min = -axis_size

# グラフオブジェクトを初期化
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
fig.suptitle('curve on the sphere', fontsize=20)

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

    # i番目のベクトルを取得
    b = b_n[i]
    tilde_b = r * b / np.linalg.norm(b) # ノルムを調整
    
    # ベクトルa,bのなす角を計算
    theta = np.arccos(np.dot(a, b) / np.linalg.norm(a) / np.linalg.norm(b))
    
    # 2点間の曲線の座標を計算
    x = beta1 * tilde_a[0] + beta2 * tilde_b[0]
    y = beta1 * tilde_a[1] + beta2 * tilde_b[1]
    z = beta1 * tilde_a[2] + beta2 * tilde_b[2]
    
    # 点ごとにノルムを計算
    norm_xyz = np.linalg.norm(np.stack([x, y, z], axis=1), axis=1)
    
    # 球面上の曲線の座標を計算
    tilde_x = r * x / norm_xyz
    tilde_y = r * y / norm_xyz
    tilde_z = r * z / norm_xyz
    
    # 球面上の2点を結ぶ曲線を作図
    ax.plot_wireframe(X, Y, Z, alpha=0.1) # 球面
    ax.quiver(0, 0, 0, *a, 
              color='red', arrow_length_ratio=l/np.linalg.norm(a)) # ベクトルa
    ax.quiver(0, 0, 0, *b, 
              color='blue', arrow_length_ratio=l/np.linalg.norm(b)) # ベクトルb
    ax.scatter(*tilde_a, 
               color='red', s=50, 
               label='$\\tilde{a}=('+', '.join(map(str, tilde_a.round(2)))+')$') # ベクトルaの球面上の点
    ax.scatter(*tilde_b, 
               color='blue', s=50, 
               label='$\\tilde{b}=('+', '.join(map(str, tilde_b.round(2)))+')$') # ベクトルbの球面上の点
    ax.plot(x, y, z, 
            color='purple', zorder=0) # 2点間の曲線
    ax.quiver(*np.tile(0, reps=(3, len(x))), 
              x, y, z, 
              color='purple', alpha=0.1, arrow_length_ratio=0.0) # 2点間の曲線の簡易塗りつぶし
    ax.plot(tilde_x, tilde_y, tilde_z, 
            color='black', zorder=0) # 球面上の2点間の曲線
    ax.quiver(*np.tile(0, reps=(3, len(x))), 
              tilde_x, tilde_y, tilde_z, 
              color='black', alpha=0.2, arrow_length_ratio=0.0) # 球面上の2点間の曲線の簡易塗りつぶし
    ax.quiver([0, tilde_a[0], tilde_b[0]], 
              [0, tilde_a[1], tilde_b[1]], 
              [0, tilde_a[2], tilde_b[2]], 
              [0, 0, 0], 
              [0, 0, 0], 
              [z_min, z_min-tilde_a[2], z_min-tilde_b[2]], 
              color=['gray', 'red', 'blue'], 
              arrow_length_ratio=0, linestyle=':') # 座標用の補助線
    ax.set_xticks(ticks=np.arange(-axis_size, axis_size+0.1, 0.5))
    ax.set_yticks(ticks=np.arange(-axis_size, axis_size+0.1, 0.5))
    ax.set_zticks(ticks=np.arange(-axis_size, axis_size+0.1, 0.5))
    ax.set_xlim(left=-axis_size, right=axis_size)
    ax.set_ylim(bottom=-axis_size, top=axis_size)
    ax.set_zlim(bottom=-axis_size, top=axis_size)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_title('$r=' + str(r) + ', ' + 
                 '\\theta=' + str((theta/np.pi).round(2)) + '\pi$', 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(filename='CurveOnSphere_3d.gif')

球面上の2点を結ぶ曲線

 曲線の内側をquiver()を使って簡易的に塗りつぶした。
 2つの曲線は、どちらも原点と点 \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を通る平面であり、平行なのを確認できる。 また、黒色の曲線は球の断面である。

平面上で可視化

 「球面上で可視化」では、「一定のノルムで2つのベクトルを結ぶ曲線」が「球面上の2点を結ぶ曲線」になるのを確認した。
 ここでは、その曲線が平面上の2点を結ぶ曲線になることをグラフで確認する。

 ベクトルや曲線(円・円弧)のサイズ(原点からのノルム) rを指定して、3次元ベクトル \mathbf{a}, \mathbf{b}を作成して、ノルムが rのベクトル \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を計算する。

# ノルム(サイズ)を指定
r = 1.2

# 3次元ベクトルを指定
a = np.array([1.0, 1.0, 1.0])
b = np.array([-1.0, -1.0, 1.0])

# ノルムを一定に調整
tilde_a = r * a / np.linalg.norm(a)
tilde_b = r * b / np.linalg.norm(b)
print(tilde_a)
print(tilde_b)
[0.69282032 0.69282032 0.69282032]
[-0.69282032 -0.69282032  0.69282032]

 「球面上で可視化」のときのコードで処理できる。

 平面上の点 \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を結ぶ曲線を計算する。

# 曲線用のラジアンを作成
#t = np.linspace(start=0.0, stop=0.5*np.pi, num=100) # なす角用
t = np.linspace(start=0.0, stop=2.0*np.pi, num=300) # 1周用

# 曲線用の係数を作成
beta1 = np.cos(t)
beta2 = np.sin(t)

# 2点間の曲線の座標を計算
x = beta1 * tilde_a[0] + beta2 * tilde_b[0]
y = beta1 * tilde_a[1] + beta2 * tilde_b[1]
z = beta1 * tilde_a[2] + beta2 * tilde_b[2]
print(x[:5])
print(y[:5])
print(z[:5])

# 点ごとにノルムを計算
xyz = np.stack([x, y, z], axis=1)
norm_xyz = np.linalg.norm(xyz, axis=1)

# ノルムを一定に調整
tilde_x = r * x / norm_xyz
tilde_y = r * y / norm_xyz
tilde_z = r * z / norm_xyz
print(tilde_x[:5])
print(tilde_y[:5])
print(tilde_z[:5])
[0.69282032 0.6781095  0.66309925 0.64779619 0.63220709]
[0.69282032 0.6781095  0.66309925 0.64779619 0.63220709]
[0.69282032 0.70722521 0.72131781 0.73509189 0.74854138]
[0.69282032 0.6829085  0.6725773  0.66181524 0.65061089]
[0.69282032 0.6829085  0.6725773  0.66181524 0.65061089]
[0.69282032 0.71223026 0.73162801 0.75100012 0.77033172]

 「球面上で可視化」のときのコードで処理できる。

 ベクトル \tilde{\mathbf{a}}, \tilde{\mathbf{b}}による平面(原点と点 \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を通る平面)の座標を計算する。

# 平面用の係数を指定
alpha = np.arange(start=-2.0, stop=2.1, step=0.5)
print(alpha)
print(alpha.shape)

# 係数の格子点を作成
Alpha1, Alpha2 = np.meshgrid(alpha, alpha)
print(Alpha1.shape)
print(Alpha2.shape)

# 平面の座標を計算
X = Alpha1 * tilde_a[0] + Alpha2 * tilde_b[0]
Y = Alpha1 * tilde_a[1] + Alpha2 * tilde_b[1]
Z = Alpha1 * tilde_a[2] + Alpha2 * tilde_b[2]
print(X[:5, :5])
print(Y[:5, :5])
print(Z[:5, :5])
[-2.  -1.5 -1.  -0.5  0.   0.5  1.   1.5  2. ]
(9,)
(9, 9)
(9, 9)
[[ 0.          0.34641016  0.69282032  1.03923048  1.38564065]
 [-0.34641016  0.          0.34641016  0.69282032  1.03923048]
 [-0.69282032 -0.34641016  0.          0.34641016  0.69282032]
 [-1.03923048 -0.69282032 -0.34641016  0.          0.34641016]
 [-1.38564065 -1.03923048 -0.69282032 -0.34641016  0.        ]]
[[ 0.          0.34641016  0.69282032  1.03923048  1.38564065]
 [-0.34641016  0.          0.34641016  0.69282032  1.03923048]
 [-0.69282032 -0.34641016  0.          0.34641016  0.69282032]
 [-1.03923048 -0.69282032 -0.34641016  0.          0.34641016]
 [-1.38564065 -1.03923048 -0.69282032 -0.34641016  0.        ]]
[[-2.77128129 -2.42487113 -2.07846097 -1.73205081 -1.38564065]
 [-2.42487113 -2.07846097 -1.73205081 -1.38564065 -1.03923048]
 [-2.07846097 -1.73205081 -1.38564065 -1.03923048 -0.69282032]
 [-1.73205081 -1.38564065 -1.03923048 -0.69282032 -0.34641016]
 [-1.38564065 -1.03923048 -0.69282032 -0.34641016  0.        ]]

 線形結合の係数 \alpha_1, \beta_2として格子状の点Alpha1, Alpha2を作成して、 \mathbf{x} = \alpha_1 \tilde{\mathbf{a}} + \alpha_2 \tilde{\mathbf{b}}を計算する。
 詳しくは「【Python】5.1:ベクトルの線形従属の可視化【『スタンフォード線形代数入門』のノート】 - からっぽのしょこ」を参照のこと。

 平面上の点 \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を結ぶ曲線のグラフを作成する。

# グラフサイズ用の値を設定
x_min = np.floor(X.min())
x_max = np.ceil(X.max())
y_min = np.floor(Y.min())
y_max = np.ceil(Y.max())
z_min = np.floor(Z.min())
z_max = np.ceil(Z.max())

# 矢のサイズを指定
l = 0.2

# 平面上の2点を結ぶ曲線を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
ax.plot_wireframe(X, Y, Z, zorder=-100) # 平面
ax.quiver(0, 0, 0, *a, 
          color='red', linewidth=4, arrow_length_ratio=l/np.linalg.norm(a)) # ベクトルa
ax.quiver(0, 0, 0, *b, 
          color='blue', linewidth=4, arrow_length_ratio=l/np.linalg.norm(b)) # ベクトルb
ax.scatter(*tilde_a, 
           color='red', s=100, 
           label='$\\tilde{a}=('+', '.join(map(str, tilde_a.round(2)))+')$', zorder=100) # ノルムを調整したベクトルaの点
ax.scatter(*tilde_b, 
           color='blue', s=100, 
           label='$\\tilde{b}=('+', '.join(map(str, tilde_b.round(2)))+')$') # ノルムを調整したベクトルbの点
ax.plot(x, y, z, 
        color='purple', linewidth=2, zorder=100) # 2点間の曲線
ax.plot(tilde_x, tilde_y, tilde_z, 
        color='black', linewidth=2, zorder=100) # ノルムを調整した2点間の曲線
ax.quiver([0, tilde_a[0], tilde_b[0]], 
          [0, tilde_a[1], tilde_b[1]], 
          [0, tilde_a[2], tilde_b[2]], 
          [0, 0, 0], 
          [0, 0, 0], 
          [z_min, z_min-tilde_a[2], z_min-tilde_b[2]], 
          color=['gray', 'red', 'blue'], 
          arrow_length_ratio=0, linestyle=':') # 座標用の補助線
ax.set_xticks(ticks=np.arange(x_min, x_max+0.1))
ax.set_yticks(ticks=np.arange(y_min, y_max+0.1))
ax.set_zticks(ticks=np.arange(z_min, z_max+0.1))
ax.set_zlim(bottom=z_min, top=z_max)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('$r=' + str(r) + '$', loc='left')
fig.suptitle('curve on the plane', fontsize=20)
ax.legend()
ax.set_aspect('equal')
#ax.view_init(elev=90, azim=270) # xy面
#ax.view_init(elev=0, azim=270) # xz面
#ax.view_init(elev=0, azim=0) # yz面
plt.show()

平面上の2点を結ぶ曲線


 グラフを回転させて確認する。

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

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

# 水平方向の角度として利用する値を作成
h_n = np.linspace(start=0.0, stop=360.0, num=frame_num+1)[:frame_num]

# グラフオブジェクトを初期化
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
fig.suptitle('curve on the plane', fontsize=20)

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

    # i番目の角度を取得
    h = h_n[i]
    
    # 平面上の2点を結ぶ曲線を作図
    ax.plot_wireframe(X, Y, Z, zorder=-100) # 平面
    ax.quiver(0, 0, 0, *a, 
              color='red', linewidth=4, arrow_length_ratio=l/np.linalg.norm(a)) # ベクトルa
    ax.quiver(0, 0, 0, *b, 
              color='blue', linewidth=4, arrow_length_ratio=l/np.linalg.norm(b)) # ベクトルb
    ax.scatter(*tilde_a, 
               color='red', s=100, 
               label='$\\tilde{a}=('+', '.join(map(str, tilde_a.round(2)))+')$', zorder=100) # ノルムを調整したベクトルaの点
    ax.scatter(*tilde_b, 
               color='blue', s=100, 
               label='$\\tilde{b}=('+', '.join(map(str, tilde_b.round(2)))+')$') # ノルムを調整したベクトルbの点
    ax.plot(x, y, z, 
            color='purple', linewidth=2, zorder=100) # 2点間の曲線
    ax.plot(tilde_x, tilde_y, tilde_z, 
            color='black', linewidth=2, zorder=100) # ノルムを調整した2点間の曲線
    ax.quiver([0, tilde_a[0], tilde_b[0]], 
              [0, tilde_a[1], tilde_b[1]], 
              [0, tilde_a[2], tilde_b[2]], 
              [0, 0, 0], 
              [0, 0, 0], 
              [z_min, z_min-tilde_a[2], z_min-tilde_b[2]], 
              color=['gray', 'red', 'blue'], 
              arrow_length_ratio=0, linestyle=':') # 座標用の補助線
    ax.set_xticks(ticks=np.arange(x_min, x_max+0.1))
    ax.set_yticks(ticks=np.arange(y_min, y_max+0.1))
    ax.set_zticks(ticks=np.arange(z_min, z_max+0.1))
    ax.set_zlim(bottom=z_min, top=z_max)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_title('$r=' + str(r) + '$', loc='left')
    ax.legend()
    ax.set_aspect('equal')
    ax.view_init(elev=30, azim=h) # 表示角度

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

# gif画像を保存
ani.save(filename='CurveOnPlane360_3d.gif')

平面上の2点を結ぶ曲線

 数式からも分かるように、平面も曲線もベクトル \tilde{\mathbf{a}}, \tilde{\mathbf{b}}を線形結合したものなので、どちらも平行になっているのを確認できる。

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

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

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

# ノルム(サイズ)を指定
r = 1.2

# 固定するベクトルを指定
a = np.array([1.0, 1.0, 1.0])
tilde_a = r * a / np.linalg.norm(a) # ノルムを調整

# 変化するベクトル用のラジアンを作成
t_n = np.linspace(start=-1.0*np.pi, stop=1.0*np.pi, num=frame_num+1)[:frame_num]
u_n = np.linspace(start=2.0/3.0*np.pi, stop=5.0/3.0*np.pi, num=frame_num+1)[:frame_num]

# 変化するベクトルのサイズを指定
d = 2.0

# 変化するベクトルを作成
b_n = np.array(
    [d * np.sin(t_n) * np.cos(u_n), d * np.sin(t_n) * np.sin(u_n), d * np.cos(t_n)]
).T

# 矢のサイズを指定
l = 0.2

# グラフサイズ用の値を設定
x_min = np.floor(np.min([np.min(Alpha1*tilde_a[0] + Alpha2*b[0]/d) for b in b_n]))
x_max =  np.ceil(np.max([np.max(Alpha1*tilde_a[0] + Alpha2*b[0]/d) for b in b_n]))
y_min = np.floor(np.min([np.min(Alpha1*tilde_a[1] + Alpha2*b[1]/d) for b in b_n]))
y_max =  np.ceil(np.max([np.max(Alpha1*tilde_a[1] + Alpha2*b[1]/d) for b in b_n]))
z_min = np.floor(np.min([np.min(Alpha1*tilde_a[2] + Alpha2*b[2]/d) for b in b_n]))
z_max =  np.ceil(np.max([np.max(Alpha1*tilde_a[2] + Alpha2*b[2]/d) for b in b_n]))

# グラフオブジェクトを初期化
fig, ax = plt.subplots(figsize=(9, 9), facecolor='white', 
                       subplot_kw={'projection': '3d'})
fig.suptitle('curve on the plane', fontsize=20)

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

    # i番目のベクトルを取得
    b = b_n[i]
    tilde_b = r * b / np.linalg.norm(b) # ノルムを調整
    
    # ベクトルa,bのなす角を計算
    theta = np.arccos(np.dot(a, b) / np.linalg.norm(a) / np.linalg.norm(b))
    
    # 2点間の曲線の座標を計算
    x = beta1 * tilde_a[0] + beta2 * tilde_b[0]
    y = beta1 * tilde_a[1] + beta2 * tilde_b[1]
    z = beta1 * tilde_a[2] + beta2 * tilde_b[2]
    
    # 点ごとにノルムを計算
    norm_xyz = np.linalg.norm(np.stack([x, y, z], axis=1), axis=1)
    
    # ノルムを一定に調整
    tilde_x = r * x / norm_xyz
    tilde_y = r * y / norm_xyz
    tilde_z = r * z / norm_xyz
    
    # 平面の座標を計算
    X = Alpha1 * tilde_a[0] + Alpha2 * tilde_b[0]
    Y = Alpha1 * tilde_a[1] + Alpha2 * tilde_b[1]
    Z = Alpha1 * tilde_a[2] + Alpha2 * tilde_b[2]
    
    # 平面上の2点を結ぶ曲線を作図
    ax.plot_wireframe(X, Y, Z, zorder=-100) # 平面
    ax.quiver(0, 0, 0, *a, 
              color='red', linewidth=4, arrow_length_ratio=l/np.linalg.norm(a)) # ベクトルa
    ax.quiver(0, 0, 0, *b, 
              color='blue', linewidth=4, arrow_length_ratio=l/np.linalg.norm(b)) # ベクトルb
    ax.scatter(*tilde_a, 
               color='red', s=100, 
               label='$\\tilde{a}=('+', '.join(map(str, tilde_a.round(2)))+')$', zorder=100) # ノルムを調整したベクトルaの点
    ax.scatter(*tilde_b, 
               color='blue', s=100, 
               label='$\\tilde{b}=('+', '.join(map(str, tilde_b.round(2)))+')$') # ノルムを調整したベクトルbの点
    ax.plot(x, y, z, 
            color='purple', linewidth=2, zorder=100) # 2点間の曲線
    ax.plot(tilde_x, tilde_y, tilde_z, 
            color='black', linewidth=2, zorder=100) # ノルムを調整した2点間の曲線
    ax.quiver([0, tilde_a[0], tilde_b[0]], 
              [0, tilde_a[1], tilde_b[1]], 
              [0, tilde_a[2], tilde_b[2]], 
              [0, 0, 0], 
              [0, 0, 0], 
              [z_min, z_min-tilde_a[2], z_min-tilde_b[2]], 
              color=['gray', 'red', 'blue'], 
              arrow_length_ratio=0, linestyle=':') # 座標用の補助線
    ax.set_xticks(ticks=np.arange(x_min, x_max+0.1))
    ax.set_yticks(ticks=np.arange(y_min, y_max+0.1))
    ax.set_zticks(ticks=np.arange(z_min, z_max+0.1))
    ax.set_xlim(left=x_min, right=x_max)
    ax.set_ylim(bottom=y_min, top=y_max)
    ax.set_zlim(bottom=z_min, top=z_max)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_title('$r=' + str(r) + ', ' + 
                 '\\theta=' + str((theta/np.pi).round(2)) + '\pi$', 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(filename='CurveOnPlane_3d.gif')

平面上の2点を結ぶ曲線

 この平面は、「球面上で可視化」における球の断面を拡張した平面である。

 この記事では、3次元空間において良い感じの曲線を引けた。

おわりに

 一度諦めたことができるようになりました。その嬉しさをなんとか伝えたい思いで(?)、4節も使って(ほぼ同じ内容を)解説しています。

 この記事の投稿日に公開された歌唱動画を聴きましょう。

 まーちゃんの供給が盛りだくさんで浮かれ気分な今日この頃です。