からっぽのしょこ

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

3次元空間における原点と2点を通る平面上の点を求めたい

はじめに

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

 この記事では、3次元空間上の平面上の点の座標を計算します。

【目次】

3次元空間における原点と2点を通る平面上の点を求めたい

 この記事では、3次元空間において、原点と任意の2点を通る平面(2つのベクトルに平行な平面)上の点に関して、x軸とy軸の値が与えられたときのz軸の値を計算したい。
 「3次元空間における3点を通る平面上の点を求めたい - からっぽのしょこ」では、任意の3点を通る平面上の点を考える。

 検索したら色々出てきて分からなかったので、自分で考えてみた。もっといい方法があれば教えてほしい。

平面上の点の計算式の導出

 まずは、3次元空間における平面上の点について、x軸とy軸の値が与えられたときのz軸の値の計算式を導出する。

 原点を O、任意の2点を A, Bとする。
 点 A, Bの座標を、それぞれ次の3次元ベクトルで表す。

 \displaystyle
\mathbf{a}
    = \begin{bmatrix}
          a_1 \\ a_2 \\ a_3
      \end{bmatrix}
,\ 
\mathbf{b}
    = \begin{bmatrix}
          b_1 \\ b_2 \\ b_3
      \end{bmatrix}

 ベクトル \mathbf{a}, \mathbf{b}は、平行でなく、どちらも0ベクトルでないとする。
 原点 Oと点 A = \mathbf{a}, B = \mathbf{b}を通る平面(ベクトル \mathbf{a}, \mathbf{b}によって張られる平面)上の点を考える。

 実数(スカラ) \alpha, \betaを係数として、ベクトル \mathbf{a}, \mathbf{b}を線形結合したベクトルを \mathbf{p}とする。

 \displaystyle
\mathbf{p}
    = \alpha \mathbf{a}
      + \beta \mathbf{b}
\tag{1}

  \mathbf{p} = (x, y, z)^{\top}とすると、各成分は次の式で表わせる。

 \displaystyle
\begin{bmatrix}
    x \\ y \\ z
\end{bmatrix}
    = \alpha
      \begin{bmatrix}
          a_1 \\ a_2 \\ a_3
      \end{bmatrix}
      + \beta
        \begin{bmatrix}
          b_1 \\ b_2 \\ b_3
        \end{bmatrix}
    = \begin{bmatrix}
          \alpha a_1 + \beta b_1 \\ 
          \alpha a_2 + \beta b_2 \\ 
          \alpha a_3 + \beta b_3
      \end{bmatrix}
\tag{1'}

  \mathbf{p}を点 Pの座標とすると、点 P = \mathbf{p}はベクトル \mathbf{a}, \mathbf{b}による平面上の点になり、式(1)は次の式で表せる。

 \displaystyle
\overrightarrow{OP}
    = \alpha \overrightarrow{OA}
      + \beta \overrightarrow{OB}

 点 Pのx軸とy軸の値 x, yが与えられたときのz軸の値 zを求める。

 平面上の点 Pの座標(ベクトル \mathbf{p}の成分)は、式(1')よりそれぞれ次の式で計算できる。

 \displaystyle
\begin{cases}
    x = \alpha a_1 + \beta b_1
    &\qquad (2) \\
    y = \alpha a_2 + \beta b_2
    &\qquad (3) \\
    z = \alpha a_3 + \beta b_3
    &\qquad (4)
\end{cases}

 x軸の式(2)を \alpha, \betaそれぞれについて整理する。

 \displaystyle
\begin{cases}
    \alpha
        = \frac{x}{a_1}
          - \frac{\beta b_1}{a_1} \\
    \beta
        = \frac{x}{b_1}
          - \frac{\alpha a_1}{b_1}
\end{cases}
\tag{2'}

 同様に、y軸の式(3)を整理する。

 \displaystyle
\begin{cases}
    \alpha
        = \frac{y}{a_2}
          - \frac{\beta b_2}{a_2} \\
    \beta
        = \frac{y}{b_2}
          - \frac{\alpha a_2}{b_2}
\end{cases}
\tag{3'}

  \mathbf{a}, \mathbf{b}が0ベクトルだと0除算になってしまう。

  \alphaに関する式について式(2')を式(3')に代入して、 \betaについて整理する。

 \displaystyle
\begin{align}
\frac{x}{a_1}
- \frac{\beta b_1}{a_1}
   &= \frac{y}{a_2}
      - \frac{\beta b_2}{a_2}
\\
\Rightarrow
a_2 x
- \beta a_2 b_1
   &= a_1 y
      - \beta a_1 b_2
\\
\Rightarrow
\beta a_1 b_2
- \beta a_2 b_1
   &= a_1 y
      - a_2 x
\\
\Rightarrow
\beta
(a_1 b_2 - a_2 b_1)
   &= a_1 y - a_2 x
\\
\Rightarrow
\beta
   &= \frac{
          a_1 y - a_2 x
      }{
          a_1 b_2 - a_2 b_1
      }
\tag{5}
\end{align}

  \betaの計算式が得られた。1行目から2行目では、両辺に a_1 a_2を掛けて分母を払った。

 同様に、 \betaに関する式について式(2')を式(3')の代入して、 \alphaについて整理する。

 \displaystyle
\begin{align}
\frac{x}{b_1}
- \frac{\alpha a_1}{b_1}
   &= \frac{y}{b_2}
      - \frac{\alpha a_2}{b_2}
\\
\Rightarrow
b_2 x
- \alpha a_1 b_2
   &= b_1 y
      - \alpha a_2 b_1
\\
\Rightarrow
\alpha a_1 b_2
- \alpha a_2 b_1
   &= b_2 x
      - b_1 y
\\
\Rightarrow
\alpha
( a_1b_2 - a_2 b_1)
   &= b_2 x - b_1 y
\\
\Rightarrow
\alpha
   &= \frac{
          b_2 x - b_1 y
      }{
          a_1 b_2 - a_2 b_1
      }
\tag{6}
\end{align}

  \alphaの計算式が得られた。1行目から2行目では、両辺に b_1 b_2を掛けて分母を払った。

  \alphaの式(5)と \betaの式(6)をz軸の式(4)に代入して、式を整理する。

 \displaystyle
\begin{align}
z
   &= \alpha a_3 + \beta b_3
\tag{4}\\
   &= \frac{
          b_2 x - b_1 y
      }{
          a_1 b_2 - a_2 b_1
      }
      a_3
      + \frac{
          a_1 y - a_2 x
        }{
          a_1 b_2 - a_2 b_1
        }
        b_3
\\
   &= \frac{
          a_3 b_2 x - a_3 b_1 y
          + a_1 b_3 y - a_2 b_3 x
      }{
          a_1 b_2 - a_2 b_1
      }
\\
   &= \frac{
          (a_3 b_2 - a_2 b_3)
          x
          + (a_1 b_3 - a_3 b_1)
            y
      }{
          a_1 b_2 - a_2 b_1
      }
\\
   &= \frac{
          a_3 b_2 - a_2 b_3
      }{
          a_1 b_2 - a_2 b_1
      }
      x
      + \frac{
          a_1 b_3 - a_3 b_1
        }{
          a_1 b_2 - a_2 b_1
        }
        y
\tag{7}
\end{align}

 点 Pのx軸とy軸の値が与えられたときのz軸の値の計算式が得られた。
  \mathbf{a}, \mathbf{b}が0ベクトルだと0除算になってしまう。また、 \mathbf{a}, \mathbf{b}が平行だと、 a_i = \gamma b_iが成り立つので分母が a_1 b_2 - a_2 b_1 = \gamma a_1 a_2 - \gamma a_2 a_1 = 0になり、0除算になってしまう。
 逆に、 \mathbf{a}, \mathbf{b}が平行でなく0ベクトルでもないと分子は0にならないので、 x, yが0でないとき \mathbf{p}は0ベクトルにならない。 x, yが0のとき点 Pは原点( \mathbf{p}は0ベクトル)になる。

平面上の点の計算の可視化

 次は、3次元空間における線形結合(平面上の点の計算)と、x軸とy軸の値が与えられたときのz軸の値の計算をグラフで可視化する。
 線形結合については「【Python】1.3:ベクトルの線形結合の可視化【『スタンフォード線形代数入門』のノート】 - からっぽのしょこ」も参照のこと。

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

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


平面の計算

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

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

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

 3次元空間におけるベクトル \overrightarrow{OA} = \mathbf{a}, \overrightarrow{OB} = \mathbf{b}をグラフで確認する。

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

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

# 基底ベクトルa,bを作図
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.round(2)))+')$') # ベクトル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.text(0, 0, 0, s='O', size=15, ha='left', va='top') # 原点ラベル
ax.text(*a, s='A', size=15, ha='center', va='bottom') # 点aラベル
ax.text(*b, s='B', size=15, ha='right', va='bottom') # 点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='gray', 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')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('$O = (0, 0, 0)$', loc='left')
fig.suptitle('$A = a, B = 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面
plt.show()

基底ベクトル

 axes.quiver()でベクトルを描画する。第1・2・3引数に始点の座標、第4・5・6引数にベクトルのサイズ(移動量)を指定する。この例では、原点を始点とする。
 配列a, bの前に*を付けてアンパック(展開)して座標を指定している。

 arrow_length_ratio引数で矢のサイズ(全体に対する矢の割合)を指定できる。この例では、各ベクトルのノルムの逆数を指定することで、ベクトルのサイズに関わらず、全ての矢のサイズを統一している。
 さらに、全てのベクトルで同じ値lを掛けることで、サイズを調整する。

 係数 \alpha, \betaを指定して、ベクトル \mathbf{a}, \mathbf{b}の線形結合(平面上の点 P = \mathbf{p})を計算する。

# 係数を指定
alpha = 2.0
beta  = -0.5

# 線形結合を計算:式(1)
p = alpha * a + beta * b
print(p)
[-3.  4.  3.]

 各ベクトルのスカラー倍の和を計算する。

 ベクトルの線形結合(ベクトル \overrightarrow{OP} = \mathbf{p})をグラフで確認する。

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

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

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

# ベクトルa,bの線形結合を作図
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.round(2)))+')$') # ベクトル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.quiver(0, 0, 0, *alpha*a, 
          color='deeppink', arrow_length_ratio=l/np.linalg.norm(alpha*a), 
          label='$\\alpha a=('+', '.join(map(str, (alpha*a).round(2)))+')$') # ベクトルαa
ax.quiver(*beta*b, *alpha*a, 
          color='deeppink', linestyle=':', arrow_length_ratio=0.0) # ベクトルαaの平行線
ax.quiver(0, 0, 0, *beta*b, 
          color='deepskyblue', arrow_length_ratio=l/np.linalg.norm(beta*b), 
          label='$\\beta b=('+', '.join(map(str, (beta*b).round(2)))+')$') # ベクトルβb
ax.quiver(*alpha*a, *beta*b, 
          color='deepskyblue', linestyle=':', arrow_length_ratio=0.0) # ベクトルβbの平行線
ax.quiver(0, 0, 0, *p, 
          color='purple', arrow_length_ratio=l/np.linalg.norm(p), 
          label='$p=('+', '.join(map(str, p.round(2)))+')$') # ベクトルp
ax.text(0, 0, 0, s='O', size=15, ha='left', va='top') # 原点ラベル
ax.text(*a, s='A', size=15, ha='left', va='center') # 点aラベル
ax.text(*b, s='B', size=15, ha='right', va='bottom') # 点bラベル
ax.text(*p, s='P', size=15, ha='center', va='bottom') # 点pラベル
ax.quiver([0, a[0], b[0], alpha*a[0], beta*b[0], p[0]], 
          [0, a[1], b[1], alpha*a[1], beta*b[1], p[1]], 
          [0, a[2], b[2], alpha*a[2], beta*b[1], p[2]], 
          [0, 0, 0, 0, 0, 0], 
          [0, 0, 0, 0, 0, 0], 
          [z_min, z_min-a[2], z_min-b[2], z_min-alpha*a[2], z_min-beta*b[2], z_min-p[2]], 
          color='gray', 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')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('$\\alpha=' + str(alpha) + ', ' + 
             '\\beta=' + str(beta) + '$', loc='left')
fig.suptitle('$p = \\alpha a + \\beta b$', fontsize=20)
ax.legend(prop={'size': 9})
ax.set_aspect('equal')
#ax.view_init(elev=90, azim=270) # xy面
#ax.view_init(elev=0, azim=270) # xz面
plt.show()

線形結合

 点 \alpha Aから \beta \mathbf{b}移動、または点 \beta Bから \alpha \mathbf{a}移動した点が Pになるのを確認できる。

 ベクトル \mathbf{a}, \mathbf{b}と平行な平面の座標を計算する。

# (平面の描画用の)係数を作成
alpha_vals = np.linspace(start=-1.0, stop=3.0, num=9)
beta_vals  = np.linspace(start=-1.5, stop=2.0, num=8)
print(alpha_vals)
print(beta_vals)

# 格子状の点を作成
alpha_grid, beta_grid = np.meshgrid(alpha_vals, beta_vals)
print(alpha_grid.shape)
print(beta_grid.shape)

# ベクトルa,bと平行な平面の座標を計算:式(1')
x_grid = alpha_grid * a[0] + beta_grid * b[0]
y_grid = alpha_grid * a[1] + beta_grid * b[1]
z_grid = alpha_grid * a[2] + beta_grid * b[2]
[-1.  -0.5  0.   0.5  1.   1.5  2.   2.5  3. ]
[-1.5 -1.  -0.5  0.   0.5  1.   1.5  2. ]
(8, 9)
(8, 9)

 係数 \alpha, \betaの値を等間隔で作成してalpha_vals, beta_valsとする。
 alpha_vals, beta_valsの格子状の点(全ての組み合わせ)をnp.meshgrid()で作成してalpha_grid, beta_gridとする。
 alpha_vals, beta_valsを用いて、成分ごとに線形結合(1')の計算を行う。

 ベクトル \mathbf{a}, \mathbf{b}によって張られる平面をグラフで確認する。

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

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

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

# タイトル用の符号を設定
sgn_beta = '+' if beta >= 0.0 else ''

# ベクトルa,bと平行な平面を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
ax.plot_wireframe(x_grid, y_grid, z_grid, 
                  alpha=0.5, label='$\\alpha a + \\beta b$') # ベクトルa,bによる平面
ax.quiver(0, 0, 0, *a, 
          color='red', linewidth=3, arrow_length_ratio=l/np.linalg.norm(a), 
          label='$a=('+', '.join(map(str, a.round(2)))+')$') # ベクトルa
ax.quiver(0, 0, 0, *b, 
          color='blue', linewidth=3, arrow_length_ratio=l/np.linalg.norm(b), 
          label='$b=('+', '.join(map(str, b.round(2)))+')$') # ベクトルb
ax.quiver(0, 0, 0, *alpha*a, 
          color='deeppink', arrow_length_ratio=l/np.linalg.norm(alpha*a)) # ベクトルαa
ax.quiver(*beta*b, *alpha*a, 
          color='deeppink', linestyle=':', arrow_length_ratio=0.0) # ベクトルαaの平行線
ax.quiver(0, 0, 0, *beta*b, 
          color='deepskyblue', arrow_length_ratio=l/np.linalg.norm(beta*b)) # ベクトルβb
ax.quiver(*alpha*a, *beta*b, 
          color='deepskyblue', linestyle=':', arrow_length_ratio=0.0) # ベクトルβbの平行線
ax.quiver(0, 0, 0, *p, 
          color='purple', linewidth=3, arrow_length_ratio=l/np.linalg.norm(p), 
          label='$p=('+', '.join(map(str, p.round(2)))+')$') # ベクトルp
ax.text(0, 0, 0, s='O', size=15, ha='left', va='top') # 原点ラベル
ax.text(*a, s='A', size=15, ha='right', va='center') # 点aラベル
ax.text(*b, s='B', size=15, ha='right', va='bottom') # 点bラベル
ax.text(*p, s='P', size=15, ha='center', va='bottom') # 点pラベル
ax.quiver([0, a[0], b[0], p[0]], 
          [0, a[1], b[1], p[1]], 
          [0, a[2], b[2], p[2]], 
          [0, 0, 0, 0], 
          [0, 0, 0, 0], 
          [z_min, z_min-a[2], z_min-b[2], z_min-p[2]], 
          color='gray', 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')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('$p=' + str(alpha) + 'a' + sgn_beta + str(beta) + 'b$', loc='left')
fig.suptitle('linear combination of two vectors', fontsize=20)
ax.legend(prop={'size': 9})
ax.set_aspect('equal')
#ax.view_init(elev=90, azim=270) # xy面
#ax.view_init(elev=0, azim=270) # xz面
plt.show()

3次元空間上の平面

3次元空間上の平面

 平面(水色のグリッド線)は、係数 \alpha, \betaを変更したときの線形結合(点 P)が取り得る座標を表す。

平面上の点の計算

 点 \mathbf{p}のx軸とy軸の値 x, yを指定して、z軸の値 zを計算する。

# x軸とy軸の値を指定
x = 4.0
y = 2.0
#x, y = p[0], p[1]

# z軸の値を計算:式(7)
z = (a[2]*b[1] - a[1]*b[2]) * x
z += (a[0]*b[2] - a[2]*b[0]) * y
z /= (a[0]*b[1] - a[1]*b[0])

# 座標を格納
p = np.array([x, y, z])
print(p)
[4.         2.         2.28571429]

  x, yx, yとして値を指定する。先ほど計算したpの0番目の要素をx、1番目の要素をyとすると、z軸の値p[2]を計算できるのを確認できる。
 式(7)により zを計算してzとする。
 x, y, zを配列に格納してpとする。

  \mathbf{p}から係数 \alpha, \betaを計算する。

# 係数を計算:式(6)(5)
alpha = (b[1]*x - b[0]*y) / (a[0]*b[1] - a[1]*b[0])
beta  = (a[0]*y - a[1]*x) / (a[0]*b[1] - a[1]*b[0])
print(alpha)
print(beta)
-0.5714285714285714
-1.4285714285714286

 式(6)(5)により係数を計算してalpha, betaとする。

 確認のため、求めた係数を用いて線形結合を計算する。

# 線形結合を計算:式(1)
p = alpha * a + beta * b
print(p)
[4.         2.         2.28571429]

 値が変わらないのを確認できる。

 ベクトル \mathbf{a}, \mathbf{b}, \mathbf{p}と平行な平面の座標を計算する。

# (平面の描画用の)係数の値を作成
alpha_min = np.min([-2.0, np.floor(alpha)])
alpha_max = np.max([2.0, np.ceil(alpha)])
beta_min  = np.min([-2.0, np.floor(beta)])
beta_max  = np.max([2.0, np.ceil(beta)])
alpha_vals = np.linspace(start=alpha_min, stop=alpha_max, num=int(2*(alpha_max-alpha_min)+1))
beta_vals  = np.linspace(start=beta_min, stop=beta_max, num=int(2*(beta_max-beta_min)+1))
print(alpha_vals)
print(beta_vals)

# 格子状の点を作成
alpha_grid, beta_grid = np.meshgrid(alpha_vals, beta_vals)

# ベクトルa,bと平行な平面の座標を計算:式(1')
x_grid = alpha_grid * a[0] + beta_grid * b[0]
y_grid = alpha_grid * a[1] + beta_grid * b[1]
z_grid = alpha_grid * a[2] + beta_grid * b[2]
[-2.  -1.5 -1.  -0.5  0.   0.5  1.   1.5  2. ]
[-2.  -1.5 -1.  -0.5  0.   0.5  1.   1.5  2. ]

 先ほどと同様にして、平面の座標を計算する。こちらは、係数の範囲をいい感じに設定している。

 平面上の点 Pをグラフで確認する。

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

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

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

# タイトル用の文字列を指定
form_label = 'z = '
form_label += '\\frac{a_3 b_2 - a_2 b_3}{a_1 b_2 - a_2 b_1} x'
form_label += '+ \\frac{a_1 b_3 - a_3 b_1}{a_1 b_2 - a_2 b_1} y'

# タイトル用の符号を設定
sgn_beta = '+' if beta >= 0.0 else ''

# ベクトルa,bと平行な平面上の点を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
ax.plot_wireframe(x_grid, y_grid, z_grid, 
                  alpha=0.5, label='$\\alpha a + \\beta b$') # ベクトルa,bによる平面
ax.scatter(*p, s=50, zorder=-100) # 点p
ax.quiver(0, 0, 0, *a, 
          color='red', linewidth=3, arrow_length_ratio=l/np.linalg.norm(a), 
          label='$a=('+', '.join(map(str, a.round(2)))+')$') # ベクトルa
ax.quiver(0, 0, 0, *b, 
          color='blue', linewidth=3, arrow_length_ratio=l/np.linalg.norm(b), 
          label='$b=('+', '.join(map(str, b.round(2)))+')$') # ベクトルb
ax.quiver(0, 0, 0, *alpha*a, 
          color='deeppink', arrow_length_ratio=l/np.linalg.norm(alpha*a)) # ベクトルαa
ax.quiver(*beta*b, *alpha*a, 
          color='deeppink', linestyle=':', arrow_length_ratio=0.0) # ベクトルαaの平行線
ax.quiver(0, 0, 0, *beta*b, 
          color='deepskyblue', arrow_length_ratio=l/np.linalg.norm(beta*b)) # ベクトルβb
ax.quiver(*alpha*a, *beta*b, 
          color='deepskyblue', linestyle=':', arrow_length_ratio=0.0) # ベクトルβbの平行線
ax.quiver(0, 0, 0, *p, 
          color='purple', linewidth=3, arrow_length_ratio=l/np.linalg.norm(p), 
          label='$p=('+', '.join(map(str, p.round(2)))+')$') # ベクトルp
ax.text(0, 0, 0, s='O', size=15, ha='right', va='bottom') # 原点ラベル
ax.text(*a, s='A', size=15, ha='right', va='bottom') # 点aラベル
ax.text(*b, s='B', size=15, ha='right', va='bottom') # 点bラベル
ax.text(*p, s='P', size=15, ha='left', va='top', zorder=100) # 点pラベル
ax.quiver([0, a[0], b[0], p[0]], 
          [0, a[1], b[1], p[1]], 
          [0, a[2], b[2], p[2]], 
          [0, 0, 0, 0], 
          [0, 0, 0, 0], 
          [z_min, z_min-a[2], z_min-b[2], z_min-p[2]], 
          color='gray', 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')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('$p=(x, y, z)=' + 
             str(alpha.round(2)) + 'a' + sgn_beta + str(beta.round(2)) + 'b$', loc='left')
fig.suptitle('$' + form_label + '$', fontsize=20)
ax.legend(prop={'size': 9})
ax.set_aspect('equal')
#ax.view_init(elev=90, azim=270) # xy面
#ax.view_init(elev=0, azim=270) # xz面
plt.show()

3次元空間における平面上の点

 点 Pとベクトル \mathbf{p}が平面上の点であるのを確認できる。

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

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

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

# 水平方向の角度として利用する値を作成
h_vals = 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('$' + form_label + '$', fontsize=20)

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

    # i番目の角度を取得
    h = h_vals[i]
    
    # ベクトルa,bと平行な平面上の点を作図
    ax.plot_wireframe(x_grid, y_grid, z_grid, 
                      alpha=0.5, label='$\\alpha a + \\beta b$') # ベクトルa,bによる平面
    ax.scatter(*p, s=100) # 点p
    ax.quiver(0, 0, 0, *a, 
              color='red', linewidth=3, arrow_length_ratio=l/np.linalg.norm(a), 
              label='$a=('+', '.join(map(str, a.round(2)))+')$') # ベクトルa
    ax.quiver(0, 0, 0, *b, 
              color='blue', linewidth=3, arrow_length_ratio=l/np.linalg.norm(b), 
              label='$b=('+', '.join(map(str, b.round(2)))+')$') # ベクトルb
    ax.quiver(0, 0, 0, *alpha*a, 
              color='deeppink', arrow_length_ratio=l/np.linalg.norm(alpha*a)) # ベクトルαa
    ax.quiver(*beta*b, *alpha*a, 
              color='deeppink', linestyle=':', arrow_length_ratio=0.0) # ベクトルαaの平行線
    ax.quiver(0, 0, 0, *beta*b, 
              color='deepskyblue', arrow_length_ratio=l/np.linalg.norm(beta*b)) # ベクトルβb
    ax.quiver(*alpha*a, *beta*b, 
              color='deepskyblue', linestyle=':', arrow_length_ratio=0.0) # ベクトルβbの平行線
    ax.quiver(0, 0, 0, *p, 
              color='purple', linewidth=3, arrow_length_ratio=l/np.linalg.norm(p), 
              label='$p=('+', '.join(map(str, p.round(2)))+')$') # ベクトルp
    ax.text(0, 0, 0, s='O', size=15, ha='center', va='bottom') # 原点ラベル
    ax.text(*a, s='A', size=15, ha='center', va='bottom') # 点aラベル
    ax.text(*b, s='B', size=15, ha='center', va='bottom') # 点bラベル
    ax.text(*p, s='P', size=15, ha='center', va='bottom') # 点pラベル
    ax.quiver([0, a[0], b[0], p[0]], 
              [0, a[1], b[1], p[1]], 
              [0, a[2], b[2], p[2]], 
              [0, 0, 0, 0], 
              [0, 0, 0, 0], 
              [z_min, z_min-a[2], z_min-b[2], z_min-p[2]], 
              color='gray', 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')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_title('$p=(x, y, z)=' + 
                 str(alpha.round(2)) + 'a' + sgn_beta + str(beta.round(2)) + 'b$', loc='left')
    ax.legend(prop={'size': 9})
    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('OPonPlane_3d.gif')

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

3次元空間における平面上の点


 以上で、3次元空間における原点と任意の2点を通る平面上の点の計算を確認できた。次の記事では、任意の3点を通る平面上の点を計算を確認する。

おわりに

 線形代数の勉強のためのほんのひとネタとしてやりたかったことができなくて、調べたら色々出てきて分からなかった(面倒になった)ので、自分なりにやってみました。
 本などを読まずに我流で書いたので、変な表現や間違い・勘違いなどあると思います。教えてください。

 さて、数学が辛いときはこの曲を聴きましょう。

 プラマイパラッパー➕➖✖️➗

【次の内容】

www.anarchive-beta.com