人工衛星の力学と制御に必要な数学の基礎：座標系・ベクトル

　JAXAのリポジトリで無料公開されている「人工衛星の力学と制御宇宙機ダイナミクス・姿勢制御技術ユニット」というPDFが気になりました。WATLABブログとの相性も良く、振動関係の実践的な知識も学べそうだと思い、記事にしながら勉強を進めていきます。ここでは基礎数学の章をまとめてみます。

こんにちは。wat（@watlablog）です。ここでは体系的に振動制御を学べそうなJAXAの資料から、数学の基礎部分をまとめてみます！

目次（項目クリックでジャンプできます）

モチベーション

　ここ数年宇宙関係の書籍を読み漁っており、衛星の力学に興味を持ちました。どうやらこの分野には振動関係の知識が非常に役にたつそうです。この記事の筆者は振動分野の計算力学技術者1級ですが、これまで制御関係が弱かったと思うので、これを機に理解を深めたいところです。…と相棒のChatGPTに相談したところ、JAXAのリポジトリを紹介してもらいました。なんと以下リンク先のPDFは無料公開されています。
・人工衛星の力学と制御宇宙機ダイナミクス・姿勢制御技術ユニット
　https://jaxa.repo.nii.ac.jp/records/6101

　数学の基礎から始まっているので、この資料に沿って勉強していけばこの分野に入門できると思いました。筆者の本業ではこのレベルまではやりませんが、振動現象の問題は周りで散見され、日々相談されることも多いので勉強してみようと思います。必要に応じてPythonで原理確認しつつ進めることで、理解を深めることができると考えます。

　本記事は、人工衛星の力学と制御を学ぶシリーズの第1回です。まずは数学の基礎から始めます。理論の詳しい説明は参考資料のPDFに譲り、この記事では数式を整理するとともに、Pythonで実装・検証した結果を紹介することを目的におきます。

　それでは早速はじめていきましょう。

座標系と座標

Matplotlibで座標系を可視化してみよう！

Pythonコードを見る

　次のコードは参考文献^[1]に載っている座標系をPython/Matplotlibで書くコードです。参考までに。

from itertools import product

import matplotlib.pyplot as plt
import numpy as np


def draw_figure() -> None:
    plt.rcParams.update(
        {
            "font.family": "Times New Roman",
            "mathtext.fontset": "custom",
            "mathtext.rm": "Times New Roman",
            "mathtext.it": "Times New Roman:italic",
            "mathtext.bf": "Times New Roman:bold",
            "mathtext.cal": "Times New Roman:italic",
        }
    )

    # Oblique 2-D projections of the three mutually perpendicular axes.
    # Their 3-D meaning is conveyed by the right-handed coordinate convention.
    basis = np.array(
        [
            [-0.48, -0.50],  # axis 1: down-left
            [0.68, 0.00],  # axis 2: right
            [0.00, 0.67],  # axis 3: up
        ]
    )
    p = np.array([2.0, 2.5, 2.5])
    components = p[:, None] * basis
    vertices = {
        bits: np.sum(np.array(bits)[:, None] * components, axis=0)
        for bits in product((0, 1), repeat=3)
    }

    _, ax = plt.subplots(figsize=(5.2, 4.6), constrained_layout=True)
    ax.set_aspect("equal")
    ax.axis("off")

    # Dashed projection box around P.
    for bits in product((0, 1), repeat=3):
        for dim in range(3):
            if bits[dim] == 0:
                next_bits = list(bits)
                next_bits[dim] = 1
                a, b = vertices[bits], vertices[tuple(next_bits)]
                ax.plot(
                    [a[0], b[0]],
                    [a[1], b[1]],
                    color="0.30",
                    linewidth=0.9,
                    linestyle=(0, (2.2, 2.2)),
                    zorder=1,
                )

    # Coordinate axes, with arrowheads extending past the third tick.
    origin = np.zeros(2)
    axis_length = 3.3
    for i, direction in enumerate(basis):
        end = axis_length * direction
        ax.annotate(
            "",
            xy=end,
            xytext=origin,
            arrowprops={
                "arrowstyle": "-|>",
                "lw": 1.25,
                "color": "black",
                "mutation_scale": 13,
            },
            zorder=3,
        )

        # Unit ticks and their numeric labels.
        normal = np.array([-direction[1], direction[0]])
        normal /= np.linalg.norm(normal)
        for n in (1, 2, 3):
            point = n * direction
            half_tick = 0.075 * normal
            ax.plot(
                [point[0] - half_tick[0], point[0] + half_tick[0]],
                [point[1] - half_tick[1], point[1] + half_tick[1]],
                color="black",
                linewidth=0.8,
                zorder=4,
            )
            offset = 0.19 * normal
            ax.text(
                *(point + offset),
                str(n),
                ha="center",
                va="center",
                fontsize=11,
                zorder=4,
            )

        axis_label_pos = end + 0.18 * direction / np.linalg.norm(direction)
        if i == 0:
            axis_label_pos += np.array([-0.02, -0.02])
        elif i == 1:
            axis_label_pos += np.array([0.02, -0.02])
        else:
            axis_label_pos += np.array([0.00, 0.02])
        ax.text(
            *axis_label_pos,
            str(i + 1),
            ha="center",
            va="center",
            fontsize=14,
            fontweight="bold",
            zorder=5,
        )

    # Origin and point P.
    P = vertices[(1, 1, 1)]
    ax.scatter(*origin, s=13, color="black", zorder=6)
    ax.scatter(*P, s=16, color="black", zorder=6)
    ax.text(0.12, -0.12, r"$O$", ha="left", va="top", fontsize=14)
    ax.text(P[0] + 0.15, P[1] + 0.01, r"$P$", ha="left", va="center", fontsize=14)

    # Coordinate labels are placed beside the corresponding axis projections.
    labels = [r"$p_1$", r"$p_2$", r"$p_3$"]
    offsets = [
        np.array([-0.10, -0.25]),
        np.array([0.05, -0.22]),
        np.array([-0.26, 0.02]),
    ]
    for i, label in enumerate(labels):
        q = p[i] * basis[i]
        ax.text(*(q + offsets[i]), label, ha="center", va="center", fontsize=14)

    ax.set_xlim(-2.1, 2.6)
    ax.set_ylim(-1.95, 2.55)

    plt.show()


if __name__ == "__main__":
    draw_figure()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

from itertools import product

import matplotlib.pyplot as plt

import numpy as np

def draw_figure() -> None:

plt.rcParams.update(

{

"font.family": "Times New Roman",

"mathtext.fontset": "custom",

"mathtext.rm": "Times New Roman",

"mathtext.it": "Times New Roman:italic",

"mathtext.bf": "Times New Roman:bold",

"mathtext.cal": "Times New Roman:italic",

}

)

# Oblique 2-D projections of the three mutually perpendicular axes.

# Their 3-D meaning is conveyed by the right-handed coordinate convention.

basis = np.array(

[

[-0.48, -0.50], # axis 1: down-left

[0.68, 0.00], # axis 2: right

[0.00, 0.67], # axis 3: up

]

)

p = np.array([2.0, 2.5, 2.5])

components = p[:, None] * basis

vertices = {

bits: np.sum(np.array(bits)[:, None] * components, axis=0)

for bits in product((0, 1), repeat=3)

}

_, ax = plt.subplots(figsize=(5.2, 4.6), constrained_layout=True)

ax.set_aspect("equal")

ax.axis("off")

# Dashed projection box around P.

for bits in product((0, 1), repeat=3):

for dim in range(3):

if bits[dim] == 0:

next_bits = list(bits)

next_bits[dim] = 1

a, b = vertices[bits], vertices[tuple(next_bits)]

ax.plot(

[a[0], b[0]],

[a[1], b[1]],

color="0.30",

linewidth=0.9,

linestyle=(0, (2.2, 2.2)),

zorder=1,

)

# Coordinate axes, with arrowheads extending past the third tick.

origin = np.zeros(2)

axis_length = 3.3

for i, direction in enumerate(basis):

end = axis_length * direction

ax.annotate(

"",

xy=end,

xytext=origin,

arrowprops={

"arrowstyle": "-|>",

"lw": 1.25,

"color": "black",

"mutation_scale": 13,

zorder=3,

)

# Unit ticks and their numeric labels.

normal = np.array([-direction[1], direction[0]])

normal /= np.linalg.norm(normal)

for n in (1, 2, 3):

point = n * direction

half_tick = 0.075 * normal

ax.plot(

[point[0] - half_tick[0], point[0] + half_tick[0]],

[point[1] - half_tick[1], point[1] + half_tick[1]],

color="black",

linewidth=0.8,

zorder=4,

)

offset = 0.19 * normal

ax.text(

*(point + offset),

str(n),

ha="center",

va="center",

fontsize=11,

zorder=4,

)

axis_label_pos = end + 0.18 * direction / np.linalg.norm(direction)

if i == 0:

axis_label_pos += np.array([-0.02, -0.02])

elif i == 1:

axis_label_pos += np.array([0.02, -0.02])

else:

axis_label_pos += np.array([0.00, 0.02])

ax.text(

*axis_label_pos,

str(i + 1),

ha="center",

va="center",

fontsize=14,

fontweight="bold",

zorder=5,

)

# Origin and point P.

P = vertices[(1, 1, 1)]

ax.scatter(*origin, s=13, color="black", zorder=6)

ax.scatter(*P, s=16, color="black", zorder=6)

ax.text(0.12, -0.12, r"$O$", ha="left", va="top", fontsize=14)

ax.text(P[0] + 0.15, P[1] + 0.01, r"$P$", ha="left", va="center", fontsize=14)

# Coordinate labels are placed beside the corresponding axis projections.

labels = [r"$p_1$", r"$p_2$", r"$p_3$"]

offsets = [

np.array([-0.10, -0.25]),

np.array([0.05, -0.22]),

np.array([-0.26, 0.02]),

]

for i, label in enumerate(labels):

q = p[i] * basis[i]

ax.text(*(q + offsets[i]), label, ha="center", va="center", fontsize=14)

ax.set_xlim(-2.1, 2.6)

ax.set_ylim(-1.95, 2.55)

plt.show()

if __name__ == "__main__":

draw_figure()

　3次元の軸方向を$x, y, z$ではなく$1, 2, 3$とおいた座標系$F$において、点$P$の位置は式(1)となります。

\[ \chi_F(P)=(p_1,p_2,p_3) \tag{1} \]

　参考文献^[1]にこの表現はありませんが、ここで$\chi_F(P)$は写像を表現しており、式全体は「座標系$F$における点$P$の座標は$(p_1,p_2,p_3)$である、と読みます。

ベクトル

　座標系の軸方向を表す単位ベクトル$\vec e_1,\vec e_2,\vec e_3$を基底ベクトルと呼びます。これらを並べた配列を$\{\vec e\}$（式(2)）、ベクトル$\vec p$ の座標成分を並べた列行列を $\mathbf p$ （式(3)）と表します。

\[ \{\vec e\} = \begin{Bmatrix} \vec e_1\\ \vec e_2\\ \vec e_3 \end{Bmatrix} = \begin{Bmatrix} \vec e_1 & \vec e_2 & \vec e_3 \end{Bmatrix}^{T} \tag{2} \]

\[ \mathbf p = \begin{bmatrix} p_1\\ p_2\\ p_3 \end{bmatrix} = \begin{bmatrix} p_1 & p_2 & p_3 \end{bmatrix}^{T} \tag{3} \]

　任意のベクトル $\vec p$ は、基底ベクトルと座標成分を使って式(4)と表せます。

\[ \vec p = p_1\vec e_1+p_2\vec e_2+p_3\vec e_3 = \{\vec e\}^{T}\mathbf p = \mathbf p^{T}\{\vec e\} \tag{4} \]

　参考文献^[1]では基底ベクトル配列を波括弧 $\{\ \}$、座標成分を角括弧 $[\ ]$ 、太字で書いたベクトルは成分を表すといった使い分けがされているようです。

例題1-1

　ここで参考文献^[1]の例題1-1を解いてみます（問題文と図はPDFを確認してください）。

　まず、元の基底でベクトル $\vec p$ を表すと、式(5)になります。

\[ \vec p=p_1\vec e_1+p_2\vec e_2 \tag{5} \]

　問題の図の角度に従うと、元の基底ベクトルは回転後の基底で次の図のように成分分解できます。

　そのため式(6)で元の基底ベクトルを表せます。

\[ \begin{aligned} \vec e_1&=\cos\theta\,\vec e'_1+\sin\theta\,\vec e'_2 \\ \vec e_2&=-\sin\theta\,\vec e'_1+\cos\theta\,\vec e'_2 \end{aligned} \tag{6} \]

　式(5)に式(6)を代入すると式(7)となります。これで答えとなります。

\[ \begin{aligned} \vec p &=p_1\left(\cos\theta\,\vec e'_1+\sin\theta\,\vec e'_2\right) +p_2\left(-\sin\theta\,\vec e'_1+\cos\theta\,\vec e'_2\right)\\ &=\left(\cos\theta\,p_1-\sin\theta\,p_2\right)\vec e'_1 +\left(\sin\theta\,p_1+\cos\theta\,p_2\right)\vec e'_2 \end{aligned} \tag{7} \]

ベクトルの内積

　任意の2つのベクトル $\vec a,\vec b$ のなす角を $\theta$ とすると、内積は式(8)で定義されます。

\[ \vec a\cdot\vec b = |\vec a|\,|\vec b|\cos\theta \tag{8} \]

　内積の結果はベクトルではなくスカラーとなり、幾何学的には、$\vec a$ を $\vec b$ の方向に射影した長さ $|\vec a|\cos\theta$ に、$|\vec b|$ を掛けた量です（式(9)）。

\[ \vec a\cdot\vec b = \underbrace{|\vec a|\cos\theta}_{\vec b\text{ 方向への射影}} \,|\vec b| \tag{9} \]

　内積は「向きの揃い具合」を $\cos\theta$ で表しつつ、2つの長さの積でスケールした量です。長さが同じなら、なす角が $0^\circ$ に近いほど正に大きく、$90^\circ$ で0、$180^\circ$ に近いほど負に大きくなります。この辺の話は「Pythonでベクトルと関数の相関を計算してみる」等でやっているので、当WATLABブログではお馴染みですね。

ベクトルの外積

　外積 $\vec a\times\vec b$ は、内積とは異なり、結果がベクトルになる演算です。大きさは式(10)と定まります。

\[ |\vec a\times\vec b| = |\vec a|\,|\vec b|\sin\theta \tag{10} \]

　向きは、外積ベクトルをその大きさで割った単位ベクトルとして式(11)で表せます。

\[ \hat{\vec n} = \frac{\vec a\times\vec b}{|\vec a\times\vec b|} \tag{11} \]

　そのため外積は式(12)と書くことができます。このベクトルの向きが $\hat{\vec n}$ の向きです。なお、$\vec a,\vec b$ が平行なら外積は零ベクトルなので、向きを定めることはできません。

\[ \vec a\times\vec b = |\vec a|\,|\vec b|\sin\theta\;\hat{\vec n} \tag{12} \]

直交基底ベクトル同士の関係

　ここでは、右手系の直交座標系を構成する基底ベクトル $\vec e_1,\vec e_2,\vec e_3$ に、内積と外積を適用したときの関係を整理します。基底ベクトルは大きさが1で、互いに直交しています。

　まず、内積についてです。同じ基底ベクトル同士の内積は1、異なる基底ベクトルどうしの内積は0になります（式(13), 式(14)）。

\[ \vec e_1\cdot\vec e_1 = \vec e_2\cdot\vec e_2 = \vec e_3\cdot\vec e_3 = 1 \tag{13} \]

\[ \vec e_1\cdot\vec e_2 = \vec e_2\cdot\vec e_3 = \vec e_3\cdot\vec e_1 = 0 \tag{14} \]

　次に、外積についてです。同じ基底ベクトル同士の外積は零ベクトルになります。異なる基底ベクトル同士の外積は、右手系の規則に従って、残りの基底ベクトルと同じ向きになります（式(15), 式(16)）。

\[ \vec e_1\times\vec e_1 = \vec e_2\times\vec e_2 = \vec e_3\times\vec e_3 = \vec 0 \tag{15} \]

\[ \vec e_1\times\vec e_2=\vec e_3,\qquad \vec e_2\times\vec e_3=\vec e_1,\qquad \vec e_3\times\vec e_1=\vec e_2 \tag{16} \]

　外積の順序を逆にすると、向きが反対になるため符号が反転します（式(17)）。

\[ \vec e_2\times\vec e_1=-\vec e_3,\qquad \vec e_3\times\vec e_2=-\vec e_1,\qquad \vec e_1\times\vec e_3=-\vec e_2 \tag{17} \]

　基底ベクトル同士の内積を各要素に並べると、式 (18) の行列になります。また、各組み合わせの外積を並べると、式 (19) のベクトル配列になります。

\[ \begin{aligned} \{\vec e\}\cdot\{\vec e\}^{T} &= \begin{Bmatrix} \vec e_1\\ \vec e_2\\ \vec e_3 \end{Bmatrix} \cdot \begin{Bmatrix} \vec e_1 & \vec e_2 & \vec e_3 \end{Bmatrix}\\ &= \begin{bmatrix} \vec e_1\cdot\vec e_1 & \vec e_1\cdot\vec e_2 & \vec e_1\cdot\vec e_3\\ \vec e_2\cdot\vec e_1 & \vec e_2\cdot\vec e_2 & \vec e_2\cdot\vec e_3\\ \vec e_3\cdot\vec e_1 & \vec e_3\cdot\vec e_2 & \vec e_3\cdot\vec e_3 \end{bmatrix}\\ &= \begin{bmatrix} 1&0&0\\ 0&1&0\\ 0&0&1 \end{bmatrix} =\mathbf U_3 \end{aligned} \tag{18} \]

\[ \begin{aligned} \{\vec e\}\times\{\vec e\}^{T} &= \begin{Bmatrix} \vec e_1\\ \vec e_2\\ \vec e_3 \end{Bmatrix} \times \begin{Bmatrix} \vec e_1 & \vec e_2 & \vec e_3 \end{Bmatrix}\\ &= \begin{Bmatrix} \vec e_1\times\vec e_1 & \vec e_1\times\vec e_2 & \vec e_1\times\vec e_3\\ \vec e_2\times\vec e_1 & \vec e_2\times\vec e_2 & \vec e_2\times\vec e_3\\ \vec e_3\times\vec e_1 & \vec e_3\times\vec e_2 & \vec e_3\times\vec e_3 \end{Bmatrix}\\ &= \begin{Bmatrix} \vec 0 & \vec e_3 & -\vec e_2\\ -\vec e_3 & \vec 0 & \vec e_1\\ \vec e_2 & -\vec e_1 & \vec 0 \end{Bmatrix} \end{aligned} \tag{19} \]

スカラー三重積

　任意の3つのベクトル$\vec a, \vec b, \vec c $ に関して、式(20)に示すスカラー三重積の関係があります。

\[ \vec a\cdot(\vec b\times\vec c) = \vec b\cdot(\vec c\times\vec a) = \vec c\cdot(\vec a\times\vec b) \tag{20} \]

　幾何学的には、まず $\vec b\times\vec c$ の大きさが、$\vec b$ と $\vec c$ が張る平行四辺形の面積になります。さらに $\vec a$ をその面に垂直な方向へ射影すると、平行六面体の高さが得られます。したがって、スカラー三重積の絶対値は「底面積×高さ」、つまり平行六面体の体積です。符号は、3つのベクトルの並び方によって決まります。

Pythonでスカラー三重積を体験しよう！

Pythonコードを見る

　スカラー三重積が持つ便利な特性、計算値は平行六面体の体積となる、という内容をPythonで検証してみましょう。次のコードを実行するとMatplotlibで可視化結果とともに計算結果も描画されます。

import itertools

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d.art3d import Poly3DCollection


# 平行六面体をつくる3つのベクトル
a = np.array([2.0, 0.0, 0.0])
b = np.array([1.0, 3.0, 0.0])
c = np.array([1.0, 1.0, 4.0])

# スカラー三重積と体積
cross_bc = np.cross(b, c)
scalar_triple_product = np.dot(a, cross_bc)
volume = abs(scalar_triple_product)

# b, c が張る底面の面積と、a の底面に垂直な高さ
base_area = np.linalg.norm(cross_bc)
normal = cross_bc / base_area
signed_height = np.dot(a, normal)
height = abs(signed_height)
foot = a - signed_height * normal

# 頂点 i*a + j*b + k*c（i, j, k は 0 または 1）
vertices = {
    bits: bits[0] * a + bits[1] * b + bits[2] * c
    for bits in itertools.product((0, 1), repeat=3)
}

# 6つの面を構成する頂点
faces = [
    [(0, 0, 0), (0, 1, 0), (0, 1, 1), (0, 0, 1)],
    [(1, 0, 0), (1, 1, 0), (1, 1, 1), (1, 0, 1)],
    [(0, 0, 0), (1, 0, 0), (1, 0, 1), (0, 0, 1)],
    [(0, 1, 0), (1, 1, 0), (1, 1, 1), (0, 1, 1)],
    [(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)],
    [(0, 0, 1), (1, 0, 1), (1, 1, 1), (0, 1, 1)],
]

plt.rcParams.update({"font.size": 11, "mathtext.fontset": "stix"})

fig = plt.figure(figsize=(12.5, 7.2))
ax = fig.add_axes([0.02, 0.08, 0.65, 0.82], projection="3d")

# 平行六面体の面を半透明で描画
polygons = [[vertices[index] for index in face] for face in faces]
solid = Poly3DCollection(
    polygons,
    facecolors="#79a9d1",
    edgecolors="#34566f",
    linewidths=1.0,
    alpha=0.16,
)
ax.add_collection3d(solid)

# 12本の辺を描画
for bits in itertools.product((0, 1), repeat=3):
    for axis in range(3):
        if bits[axis] == 0:
            other = list(bits)
            other[axis] = 1

            start = vertices[bits]
            end = vertices[tuple(other)]

            ax.plot(
                [start[0], end[0]],
                [start[1], end[1]],
                [start[2], end[2]],
                color="#34566f",
                linewidth=1.0,
            )

# 原点から出る3つのベクトルを描画
vector_specs = [
    (a, "#c23b32", r"$\vec{a}$"),
    (b, "#27864b", r"$\vec{b}$"),
    (c, "#245b9b", r"$\vec{c}$"),
]

for vector, color, label in vector_specs:
    ax.quiver(
        0, 0, 0,
        vector[0], vector[1], vector[2],
        color=color,
        linewidth=2.4,
        arrow_length_ratio=0.10,
    )
    ax.text(*(vector * 1.05), label, color=color, fontsize=14)

# a の先端から底面に下ろした垂線を描画
ax.plot(
    [a[0], foot[0]],
    [a[1], foot[1]],
    [a[2], foot[2]],
    color="#8a4f9e",
    linestyle="--",
    linewidth=2.0,
)
midpoint = 0.5 * (a + foot)
ax.text(*midpoint, r"$h$", color="#8a4f9e", fontsize=14)
ax.scatter(*foot, color="#8a4f9e", s=24, depthshade=False)

# 原点を表示
ax.scatter(0, 0, 0, color="black", s=18, depthshade=False)
ax.text(0, 0, 0, "  O", color="black")

# 入力条件と計算結果を図の右側に表示
info = (
    "Input vectors\n"
    f"a = {np.array2string(a, precision=1)}\n"
    f"b = {np.array2string(b, precision=1)}\n"
    f"c = {np.array2string(c, precision=1)}\n\n"
    f"b × c = {np.array2string(cross_bc, precision=1)}\n"
    f"a · (b × c) = {scalar_triple_product:.1f}\n"
    f"Volume = |a · (b × c)| = {volume:.1f}\n"
    f"Base area × height = {base_area:.2f} × {height:.2f}"
    f" = {base_area * height:.1f}"
)
fig.text(
    0.70, 0.82, info,
    va="top",
    family="monospace",
    fontsize=10.5,
    linespacing=1.55,
    bbox={
        "boxstyle": "round,pad=0.7",
        "facecolor": "#f7f9fb",
        "edgecolor": "0.6",
    },
)

# 軸・視点などの設定
fig.suptitle("Scalar triple product and parallelepiped volume", y=0.97, fontsize=15)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.set_xlim(-0.4, 4.5)
ax.set_ylim(-0.4, 4.5)
ax.set_zlim(-0.4, 4.5)
ax.set_box_aspect((1, 1, 1))
ax.view_init(elev=24, azim=-55)
ax.grid(True, alpha=0.25)

plt.show()

print(f"Scalar triple product: {scalar_triple_product:.1f}")
print(f"Parallelepiped volume: {volume:.1f}")
print(
    f"Base area * height: "
    f"{base_area:.2f} * {height:.2f} = {base_area * height:.1f}"
)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

import itertools

import matplotlib.pyplot as plt

import numpy as np

from mpl_toolkits.mplot3d.art3d import Poly3DCollection

# 平行六面体をつくる3つのベクトル

a = np.array([2.0, 0.0, 0.0])

b = np.array([1.0, 3.0, 0.0])

c = np.array([1.0, 1.0, 4.0])

# スカラー三重積と体積

cross_bc = np.cross(b, c)

scalar_triple_product = np.dot(a, cross_bc)

volume = abs(scalar_triple_product)

# b, c が張る底面の面積と、a の底面に垂直な高さ

base_area = np.linalg.norm(cross_bc)

normal = cross_bc / base_area

signed_height = np.dot(a, normal)

height = abs(signed_height)

foot = a - signed_height * normal

# 頂点 i*a + j*b + k*c（i, j, k は 0 または 1）

vertices = {

bits: bits[0] * a + bits[1] * b + bits[2] * c

for bits in itertools.product((0, 1), repeat=3)

}

# 6つの面を構成する頂点

faces = [

[(0, 0, 0), (0, 1, 0), (0, 1, 1), (0, 0, 1)],

[(1, 0, 0), (1, 1, 0), (1, 1, 1), (1, 0, 1)],

[(0, 0, 0), (1, 0, 0), (1, 0, 1), (0, 0, 1)],

[(0, 1, 0), (1, 1, 0), (1, 1, 1), (0, 1, 1)],

[(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)],

[(0, 0, 1), (1, 0, 1), (1, 1, 1), (0, 1, 1)],

]

plt.rcParams.update({"font.size": 11, "mathtext.fontset": "stix"})

fig = plt.figure(figsize=(12.5, 7.2))

ax = fig.add_axes([0.02, 0.08, 0.65, 0.82], projection="3d")

# 平行六面体の面を半透明で描画

polygons = [[vertices[index] for index in face] for face in faces]

solid = Poly3DCollection(

polygons,

facecolors="#79a9d1",

edgecolors="#34566f",

linewidths=1.0,

alpha=0.16,

)

ax.add_collection3d(solid)

# 12本の辺を描画

for bits in itertools.product((0, 1), repeat=3):

for axis in range(3):

if bits[axis] == 0:

other = list(bits)

other[axis] = 1

start = vertices[bits]

end = vertices[tuple(other)]

ax.plot(

[start[0], end[0]],

[start[1], end[1]],

[start[2], end[2]],

color="#34566f",

linewidth=1.0,

)

# 原点から出る3つのベクトルを描画

vector_specs = [

(a, "#c23b32", r"$\vec{a}$"),

(b, "#27864b", r"$\vec{b}$"),

(c, "#245b9b", r"$\vec{c}$"),

]

for vector, color, label in vector_specs:

ax.quiver(

0, 0, 0,

vector[0], vector[1], vector[2],

color=color,

linewidth=2.4,

arrow_length_ratio=0.10,

)

ax.text(*(vector * 1.05), label, color=color, fontsize=14)

# a の先端から底面に下ろした垂線を描画

ax.plot(

[a[0], foot[0]],

[a[1], foot[1]],

[a[2], foot[2]],

color="#8a4f9e",

linestyle="--",

linewidth=2.0,

)

midpoint = 0.5 * (a + foot)

ax.text(*midpoint, r"$h$", color="#8a4f9e", fontsize=14)

ax.scatter(*foot, color="#8a4f9e", s=24, depthshade=False)

# 原点を表示

ax.scatter(0, 0, 0, color="black", s=18, depthshade=False)

ax.text(0, 0, 0, " O", color="black")

# 入力条件と計算結果を図の右側に表示

info = (

"Input vectors\n"

f"a = {np.array2string(a, precision=1)}\n"

f"b = {np.array2string(b, precision=1)}\n"

f"c = {np.array2string(c, precision=1)}\n\n"

f"b × c = {np.array2string(cross_bc, precision=1)}\n"

f"a · (b × c) = {scalar_triple_product:.1f}\n"

f"Volume = |a · (b × c)| = {volume:.1f}\n"

f"Base area × height = {base_area:.2f} × {height:.2f}"

f" = {base_area * height:.1f}"

)

fig.text(

0.70, 0.82, info,

va="top",

family="monospace",

fontsize=10.5,

linespacing=1.55,

bbox={

"boxstyle": "round,pad=0.7",

"facecolor": "#f7f9fb",

"edgecolor": "0.6",

)

# 軸・視点などの設定

fig.suptitle("Scalar triple product and parallelepiped volume", y=0.97, fontsize=15)

ax.set_xlabel("x")

ax.set_ylabel("y")

ax.set_zlabel("z")

ax.set_xlim(-0.4, 4.5)

ax.set_ylim(-0.4, 4.5)

ax.set_zlim(-0.4, 4.5)

ax.set_box_aspect((1, 1, 1))

ax.view_init(elev=24, azim=-55)

ax.grid(True, alpha=0.25)

plt.show()

print(f"Scalar triple product: {scalar_triple_product:.1f}")

print(f"Parallelepiped volume: {volume:.1f}")

print(

f"Base area * height: "

f"{base_area:.2f} * {height:.2f} = {base_area * height:.1f}"

)

　こちらが実行結果です。スカラー三重積で計算された値と底面積×高さで計算された値は一致しました。

ベクトル三重積

　ベクトル三重積は、外積の結果にもう一度外積を適用する演算です。一般に、次の公式(21)で計算できます。

\[ \vec{a}\times(\vec{b}\times\vec{c}) = \vec{b}(\vec{a}\cdot\vec{c}) - \vec{c}(\vec{a}\cdot\vec{b}) \tag{21} \]

　右辺はベクトル $\vec b,\vec c$ のスカラー倍どうしの差なので、結果もベクトルになります。この公式は、ベクトル三重積の計算を内積とベクトルのスカラー倍に分けて扱えるため便利です。

まとめ

　本記事では、人工衛星の力学と制御を学ぶための数学的準備として、座標系とベクトルの基礎を整理しました。内積はベクトルの向きの関係をスカラーで表し、外積は大きさと向きを持つベクトルを与えます。また、スカラー三重積の絶対値は、3つのベクトルがつくる平行六面体の体積に等しいことを、Pythonで可視化して確認しました。

　ベクトル三重積では、外積を含む式を内積とベクトルのスカラー倍に変形できます。これらの関係は、今後、座標変換や力学・制御の式を扱う際にも使われます。

参考文献

[1] JAXA, 人工衛星の力学と制御：宇宙機ダイナミクス・姿勢制御技術ユニット, 2006, JAXA-SP-05-025.

GitHub

　本記事のコードは以下のGitHubに置きました。
・GitHub：https://github.com/watlablog/python-spacecraft-dynamics

異分野を学ぶためにJAXAの資料を読み始めました！結構しっかりと、しかも専門的な知識を学べそうです！
Xでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！

人工衛星の力学と制御に必要な数学の基礎：座標系・ベクトル

モチベーション

座標系と座標

Matplotlibで座標系を可視化してみよう！

ベクトル

例題1-1

ベクトルの内積

ベクトルの外積

直交基底ベクトル同士の関係

スカラー三重積

Pythonでスカラー三重積を体験しよう！

ベクトル三重積

まとめ

参考文献

GitHub

コメントを残すコメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

モチベーション

座標系と座標

Matplotlibで座標系を可視化してみよう！

ベクトル

例題1-1

ベクトルの内積

ベクトルの外積

直交基底ベクトル同士の関係

スカラー三重積

Pythonでスカラー三重積を体験しよう！

ベクトル三重積

まとめ

参考文献

GitHub

SNSでもご購読できます。

コメントを残す コメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

コメントを残すコメントをキャンセル