人工衛星の力学と制御に必要な数学の基礎：ダイアディック

　ベクトルから特定方向の成分を取り出す射影を出発点に、ダイアディックの考え方を基礎から解説します。数式の意味を一つずつ確認しながら、NumPyによる計算とMatplotlibでの可視化を経て、人工衛星の力学で扱う3次元ベクトルへのつながりまで学びます。文章とコードの両面から理解を深めましょう。

こんにちは。wat（@watlablog）です。ここでは数学のダイアディックについて理解を深めてみます！

　本記事は参考文献^[1]に沿って個人的な勉強のためにお試し計算をしているという内容です。専門家による記事ではないのでご注意ください。
・前回記事：人工衛星の力学と制御に必要な数学の基礎：座標系・ベクトル

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

ダイアディックとは？

2階のテンソル量

　ダイアディックは、ベクトルに作用して、別のベクトルをつくる2階のテンソルです。「2階」は作用を二度行うという意味ではなく、成分を表すときに二つの添字が必要な量であることを示します。まず、量の階数を簡単に整理してみましょう。

階数	量の例	特徴
0階	質量、温度などのスカラー	大きさだけを持ち、方向を持たない
1階	位置、力、角速度などのベクトル	大きさと方向を持つ
2階	慣性テンソルなど	二つの方向の関係を表す

　2階テンソルの成分は、たとえば $D_{ij}$ のように二つの添字で表されます。$D_{ij}$ は、入力ベクトルの第 $j$ 方向の成分が、出力ベクトルの第 $i$ 方向にどう関係するかを示します。そのため、二階テンソルを使うと、方向どうしの関係を一つの式にまとめて表せます。

　ここでいう「階数」は、行列のランクとは別の概念です。2階というのは、テンソルの成分を表すのに二つの添字が必要だという意味です。

慣性テンソルに見る2階テンソルの働き

　2階テンソルの具体例として、剛体の慣性テンソルを考えます。剛体の角速度を $\boldsymbol{\omega}$、角運動量を $\mathbf H$、慣性テンソルを $\mathbf I$ とすると、両者の関係は

\[ \mathbf H=\mathbf I\boldsymbol{\omega} \]

と表されます。

　慣性テンソル $\mathbf I$ は、剛体の質量分布を反映し、角速度ベクトルを角運動量ベクトルへ結び付けます。剛体の形や質量分布によっては、角速度と角運動量が同じ方向を向かないことがあります。つまり、ある方向の回転に対する剛体の応答を、回転軸ごとの独立した一つの数だけでなく、方向どうしの関係を含めて表す必要があります。

　このように、慣性テンソルを使えば、回転を表す角速度と角運動量の方向的な関係を一つの式で扱えます。この関係は人工衛星だけのものではありません。航空機やロボットなど、剛体の姿勢運動や姿勢制御を扱う場面でも基本になります。

　なお、慣性テンソルの成分行列は、どの座標系・基底で表すかによって数値が変わります。一方で、表している物理的な作用は同じです。本記事ではこの違いを踏まえ、ダイアディックそのものと、その成分行列を区別しながら見ていきます。

この記事で学ぶこと

　本記事では、まずベクトルを特定の方向へ射影する例から始めます。次に、入力ベクトルの各方向成分を組み合わせて出力ベクトルをつくる考え方へ広げ、参考文献^[1]の記法に沿ってダイアディック積と二階テンソルを説明します。

　最後に、NumPyで三次元のダイアディックを計算し、Matplotlibで入力ベクトル、出力ベクトル、各入力成分からの寄与を可視化します。数式と図、コードを対応させながら、ダイアディックが方向どうしの関係をどのように表しているのかを確かめていきましょう。

ベクトルを別のベクトルに変換したいとき

　ダイアディックの定義式に入る前に、まず「ベクトルに何らかの操作を加えて、別のベクトルを得る」とはどういうことかを考えてみましょう。

たとえば、あるベクトルから特定の方向の成分だけを取り出したいとします。入力するベクトルを $\vec v$、操作の結果得られるベクトルを $\vec w$ とすると、この関係は

\[ \vec w = \text{ベクトルに加える操作}(\vec v) \]

と表せます。この「操作」を行列 $\mathbf A$ で表すと、

\[ \vec w = \mathbf A\vec v \]

と書けます。行列 $\mathbf A$ が入力ベクトル $\vec v$ に作用し、その結果として出力ベクトル $\vec w$ が得られる、という関係です。

　このような操作は、人工衛星の力学や制御でも使われます。たとえば、ベクトルを特定の方向へ射影したり、ベクトルの各成分を組み合わせたりする計算です。重要なのは、ベクトルの成分をただ別の座標系で表し直すこととは区別することです。ここで考えるのは、入力されたベクトルに操作を加え、新しいベクトルを作ることです。

　次に、もっとも簡単な例として、2次元ベクトルから $x$ 方向の成分だけを取り出してみます。この例を通して、「ベクトルに作用するもの」を、まず具体的に見ていきましょう。

例題：ベクトルから特定方向の成分を取り出す

　ここでは、2次元ベクトルから $x$ 方向の成分だけを取り出します。入力するベクトルを

\[ \vec v= \begin{pmatrix} 3\\ 4 \end{pmatrix} = 3\vec e_x+4\vec e_y \]

とします。$\vec e_x=(1,0)^T$、$\vec e_y=(0,1)^T$ は、それぞれ $x$ 方向、$y$ 方向の単位基底ベクトルです。このベクトルには、$x$ 方向に $3$、$y$ 方向に $4$ の成分があります。したがって、$x$ 方向の成分だけを取り出すと、その成分の大きさは $3$ です。

　$x$ 方向の成分の大きさは、単位ベクトル $\vec e_x$ との内積で求められます。

\[ \vec e_x\cdot\vec v = \begin{pmatrix} 1\\ 0 \end{pmatrix} \cdot \begin{pmatrix} 3\\ 4 \end{pmatrix} =3 \]

　内積の結果 $3$ はスカラーです。これを $x$ 方向のベクトルとして表すには、$\vec e_x$ を掛けます。

\[ \vec v_x = (\vec e_x\cdot\vec v)\vec e_x = 3 \begin{pmatrix} 1\\ 0 \end{pmatrix} = \begin{pmatrix} 3\\ 0 \end{pmatrix} \]

　したがって、$x$ 方向の成分の大きさは $3$、その成分をベクトルで表したものは $(3,0)^T$ です。ここでは、内積で方向成分の大きさを取り出し、単位ベクトルを使ってその方向のベクトルに戻すという二段階で射影を求めました。次はこの計算を図にして、入力ベクトルと $x$ 方向への射影の関係を確認します。

Pythonで射影を確認

　先ほど計算した入力ベクトル $\vec v=(3,4)^T$ と、$x$ 方向への射影ベクトル $\vec v_x=(3,0)^T$ を図にします。入力ベクトルの終点から $x$ 軸へ垂線を引くと、その足が射影ベクトルの終点になります。

Pythonコードを見る

import matplotlib.pyplot as plt

# Times New Romanを基本フォントに設定
plt.rcParams["font.family"] = "Times New Roman"

# ベクトルの始点と終点
origin = (0, 0)
v = (3, 4)       # 入力ベクトル
v_x = (3, 0)     # x方向への射影ベクトル

fig, ax = plt.subplots(figsize=(6, 5))

# 入力ベクトル v = (3, 4)
ax.annotate(
    "",
    xy=v,
    xytext=origin,
    arrowprops=dict(arrowstyle="->", color="C0", linewidth=2),
)
ax.text(1.25, 2.35, r"$\vec{v}=(3,4)^T$", color="C0", fontsize=13)

# x方向への射影 v_x = (3, 0)
ax.annotate(
    "",
    xy=v_x,
    xytext=origin,
    arrowprops=dict(arrowstyle="->", color="C1", linewidth=2),
)
ax.text(1.35, -0.45, r"$\vec{v}_x=(3,0)^T$", color="C1", fontsize=13)

# 入力ベクトルの終点からx軸への垂線
ax.plot(
    [v_x[0], v[0]],
    [v_x[1], v[1]],
    linestyle="--",
    color="0.4",
    linewidth=1.2,
)

# 射影点を表示
ax.plot(v_x[0], v_x[1], "o", color="C1", markersize=5)

# x軸とy軸を原点に配置
ax.spines["left"].set_position("zero")
ax.spines["bottom"].set_position("zero")
ax.spines["right"].set_color("none")
ax.spines["top"].set_color("none")

# 上下左右に目盛を表示し、内向きにする
ax.tick_params(
    axis="both",
    which="both",
    direction="inout",
    top=False,
    right=False,
)

ax.set_xlim(-0.5, 4.5)
ax.set_ylim(-0.8, 5.0)
ax.set_aspect("equal", adjustable="box")
ax.set_xlabel(
    "x",
    loc="right",
    fontsize=13,
    fontname="Times New Roman",
    fontstyle="italic",
)
ax.set_ylabel(
    "y",
    loc="top",
    fontsize=13,
    fontname="Times New Roman",
    fontstyle="italic",
)
ax.grid(True, linestyle=":", linewidth=0.7, alpha=0.6)

fig.tight_layout()
plt.show()

import matplotlib.pyplot as plt

# Times New Romanを基本フォントに設定

plt.rcParams["font.family"] = "Times New Roman"

# ベクトルの始点と終点

origin = (0, 0)

v = (3, 4) # 入力ベクトル

v_x = (3, 0) # x方向への射影ベクトル

fig, ax = plt.subplots(figsize=(6, 5))

# 入力ベクトル v = (3, 4)

ax.annotate(

"",

xy=v,

xytext=origin,

arrowprops=dict(arrowstyle="->", color="C0", linewidth=2),

)

ax.text(1.25, 2.35, r"$\vec{v}=(3,4)^T$", color="C0", fontsize=13)

# x方向への射影 v_x = (3, 0)

ax.annotate(

"",

xy=v_x,

xytext=origin,

arrowprops=dict(arrowstyle="->", color="C1", linewidth=2),

)

ax.text(1.35, -0.45, r"$\vec{v}_x=(3,0)^T$", color="C1", fontsize=13)

# 入力ベクトルの終点からx軸への垂線

ax.plot(

[v_x[0], v[0]],

[v_x[1], v[1]],

linestyle="--",

color="0.4",

linewidth=1.2,

)

# 射影点を表示

ax.plot(v_x[0], v_x[1], "o", color="C1", markersize=5)

# x軸とy軸を原点に配置

ax.spines["left"].set_position("zero")

ax.spines["bottom"].set_position("zero")

ax.spines["right"].set_color("none")

ax.spines["top"].set_color("none")

# 上下左右に目盛を表示し、内向きにする

ax.tick_params(

axis="both",

which="both",

direction="inout",

top=False,

right=False,

)

ax.set_xlim(-0.5, 4.5)

ax.set_ylim(-0.8, 5.0)

ax.set_aspect("equal", adjustable="box")

ax.set_xlabel(

"x",

loc="right",

fontsize=13,

fontname="Times New Roman",

fontstyle="italic",

)

ax.set_ylabel(

"y",

loc="top",

fontsize=13,

fontname="Times New Roman",

fontstyle="italic",

)

ax.grid(True, linestyle=":", linewidth=0.7, alpha=0.6)

fig.tight_layout()

plt.show()

　青い入力ベクトルの終点は $(3,4)$ です。そこから $x$ 軸へ下ろした破線の足は $(3,0)$ となり、これが $x$ 方向への射影ベクトル $\vec v_x$ の終点です。図からも、入力ベクトルの $x$ 成分が $3$ であり、射影ベクトルが $(3,0)^T$ になることを確認できます。

射影の考え方を一般化する

入力ベクトルの成分から出力ベクトルを作る

　前の例では、入力ベクトルから $x$ 方向の成分だけを取り出しました。これは、入力ベクトルの各方向成分のうち、$x$ 方向の成分だけを出力に使う操作となります。ここではこの考え方を広げ、入力ベクトルの各方向成分を取り出し、それらを組み合わせて出力ベクトルを作ることを考えます。

直交規格化された基底 $\vec e_1,\vec e_2$ で表される2次元の入力ベクトルを

\[ \vec v=v_1\vec e_1+v_2\vec e_2 \]

とします。このベクトルから、基底 $\vec e_1,\vec e_2$ の方向成分を取り出すと、それぞれ

\[ \vec e_1\cdot\vec v=v_1,\qquad \vec e_2\cdot\vec v=v_2 \]

です。
これらの成分に係数を掛け、出力ベクトルの各方向に割り当てることにします。たとえば、出力ベクトル $\vec w$ の成分を

\[ \begin{aligned} w_1 &= D_{11}v_1+D_{12}v_2,\\ w_2 &= D_{21}v_1+D_{22}v_2 \end{aligned} \]

と定めます。係数 $D_{ij}$ は、入力の第 $j$ 方向の成分 $v_j$ が、出力の第 $i$ 方向の成分 $w_i$ にどれだけ寄与するかを表します。

この関係は、成分行列を使うと

\[ \begin{pmatrix} w_1\\ w_2 \end{pmatrix} = \begin{pmatrix} D_{11}&D_{12}\\ D_{21}&D_{22} \end{pmatrix} \begin{pmatrix} v_1\\ v_2 \end{pmatrix} \]

と書けます。行列の掛け算を成分ごとに書き下すと、先ほどの $w_1,w_2$ の式になります。

基底ベクトルを使って書く

　出力ベクトルは、出力の各成分を基底ベクトルに沿って足し合わせたものです。

\[ \vec w=w_1\vec e_1+w_2\vec e_2 \]

出力成分 $w_1,w_2$ の定義を代入すると、

\[ \begin{aligned} \vec w &=(D_{11}v_1+D_{12}v_2)\vec e_1 +(D_{21}v_1+D_{22}v_2)\vec e_2\\ &=D_{11}v_1\vec e_1 +D_{12}v_2\vec e_1 +D_{21}v_1\vec e_2 +D_{22}v_2\vec e_2 \end{aligned} \]

となります。

ここで、基底 $\vec e_1,\vec e_2$ は直交規格化されているため、入力ベクトルの各成分は、対応する基底ベクトルとの内積で取り出せます。

\[ v_1=\vec e_1\cdot\vec v,\qquad v_2=\vec e_2\cdot\vec v \]

したがって、$v_1,v_2$ は成分のまま書いても構いませんが、ここでは式を基底ベクトルと入力ベクトルを使った形に書き換えるため、あえて内積で表します。たとえば、

\[ D_{11}v_1\vec e_1 = D_{11}(\vec e_1\cdot\vec v)\vec e_1 \]

と書き換えられます。このように書くと、入力ベクトルからどの方向の成分を取り出し、それをどの方向の出力に使うのかを、基底ベクトルを使って表せます。

すべての成分を同じように書き換えると、

\[ \begin{aligned} \vec w ={}&D_{11}(\vec e_1\cdot\vec v)\vec e_1 +D_{12}(\vec e_2\cdot\vec v)\vec e_1\\ &+D_{21}(\vec e_1\cdot\vec v)\vec e_2 +D_{22}(\vec e_2\cdot\vec v)\vec e_2 \end{aligned} \]

となります。たとえば、第2項の

\[ D_{12}(\vec e_2\cdot\vec v)\vec e_1 \]

では、まず内積 $\vec e_2\cdot\vec v$ で、入力ベクトルの $\vec e_2$ 方向の成分を取り出します。その成分に係数 $D_{12}$ を掛け、$\vec e_1$ 方向の出力に加えています。

つまり、係数 $D_{ij}$ は、入力の第 $j$ 方向の成分を出力の第 $i$ 方向にどれだけ反映させるかを表しています。ここでは入力と出力を、同じ基底 $\vec e_1,\vec e_2$ で表しています。

3次元のベクトルに拡張する

　人工衛星の力学では、位置、速度、力、角速度などを3次元ベクトルとして扱います。そこで、同じ考え方を3次元に広げます。直交規格化された基底 $\vec e_1,\vec e_2,\vec e_3$ に対して、入力ベクトルを

\[ \vec v=v_1\vec e_1+v_2\vec e_2+v_3\vec e_3 \]

と表します。

このとき、出力ベクトルの各成分 $w_i$ は、入力ベクトルの各成分 $v_j$ を係数 $D_{ij}$ に従って組み合わせたものになります。

\[ w_i=\sum_{j=1}^{3}D_{ij}v_j \qquad (i=1,2,3) \]

成分をすべて書き下すと、

\[ \begin{aligned} w_1&=D_{11}v_1+D_{12}v_2+D_{13}v_3,\\ w_2&=D_{21}v_1+D_{22}v_2+D_{23}v_3,\\ w_3&=D_{31}v_1+D_{32}v_2+D_{33}v_3 \end{aligned} \]

です。これを成分行列で書けば、

\[ \begin{pmatrix} w_1\\ w_2\\ w_3 \end{pmatrix} = \begin{pmatrix} D_{11}&D_{12}&D_{13}\\ D_{21}&D_{22}&D_{23}\\ D_{31}&D_{32}&D_{33} \end{pmatrix} \begin{pmatrix} v_1\\ v_2\\ v_3 \end{pmatrix} \]

となります。

2次元の例と同じく、入力の各方向の成分を取り出し、係数 $D_{ij}$ に従って出力の各方向に割り当てています。次は、この関係を参考文献^[1]の記法に合わせて、基底ベクトルを使ったダイアディックの式で表します。

ダイアディック

ダイアディックの定義

　前の節では入力ベクトルの各方向成分を取り出し、係数 $D_{ij}$ に従って出力ベクトルの各方向に割り当てることを考えました。ここでは、その関係を参考文献^[1]の記法に合わせて表します。

　3次元の基底ベクトルを $\vec e_1,\vec e_2,\vec e_3$、係数を $D_{ij}$ とすると、参考文献^[1]の式(1.3.1)は次の形で表せます。

\[ \vec{\mathbf D} = \sum_{i=1}^{3}\sum_{j=1}^{3} D_{ij}\,\vec e_i\vec e_j \]

　ここで、$\vec{\mathbf D}$ はダイアディック、$D_{ij}$ はその成分を表す数値です。矢印付きの太字 $\vec{\mathbf D}$ と、矢印のない太字 $\mathbf D$ を区別します。後者は、係数 $D_{ij}$ を並べた行列です。

\[ \mathbf D= \begin{pmatrix} D_{11}&D_{12}&D_{13}\\ D_{21}&D_{22}&D_{23}\\ D_{31}&D_{32}&D_{33} \end{pmatrix} \]

　式にある $\vec e_i\vec e_j$ は、通常の数値の掛け算でも、内積や外積でもありません。ダイアディック積と呼ばれる、二つのベクトルから作る量です。この積をベクトル $\vec v$ に作用させると、次のようになります。

\[ (\vec e_i\vec e_j)\cdot\vec v = \vec e_i(\vec e_j\cdot\vec v) \]

右側から順に見ると、まず内積 $\vec e_j\cdot\vec v$ で入力ベクトルの $\vec e_j$ 方向の成分を取り出します。その成分を $\vec e_i$ 方向のベクトルにします。つまり、$\vec e_i\vec e_j$ は、入力の第 $j$ 方向成分を取り出して、第 $i$ 方向の出力にする項です。

Pythonでダイアディック計算をする

　それではここまで学んだ内容をPythonでコーディングしてみましょう。次のPythonコードを説明します。

Pythonコードを見る

import numpy as np
import matplotlib.pyplot as plt

# Times New Romanを基本フォントに設定
plt.rcParams["font.family"] = "Times New Roman"

# 標準基底ベクトル e_1, e_2, e_3
basis = np.eye(3)

# ダイアディックの係数 D_ij
D = np.array([
    [1.0,  2.0, 0.0],
    [0.0, -1.0, 1.0],
    [2.0,  0.0, 1.0],
])

# 入力ベクトル
v = np.array([1.0, 2.0, -1.0])

# 定義式の各項 D_ij e_i e_j を
# 標準基底上の行列として足し合わせる
D_dyadic = np.zeros((3, 3))

for i in range(3):
    for j in range(3):
        D_dyadic += D[i, j] * np.outer(basis[i], basis[j])

# ダイアディックを入力ベクトルに作用させる
w = D_dyadic @ v

# 入力の各成分から出力への寄与ベクトル c_j
contributions = np.array([
    D[:, j] * v[j]
    for j in range(3)
])

print("入力ベクトル v:", v)
print("出力ベクトル w:", w)
print("各寄与ベクトル c_j:\n", contributions)

fig = plt.figure(figsize=(8, 7))
ax = fig.add_subplot(111, projection="3d")


def draw_vector(start, vector, color, label, linestyle="-", linewidth=2):
    """startを始点とする3次元ベクトルを描く。"""
    start = np.asarray(start, dtype=float)
    vector = np.asarray(vector, dtype=float)

    ax.quiver(
        start[0], start[1], start[2],
        vector[0], vector[1], vector[2],
        color=color,
        linewidth=linewidth,
        linestyle=linestyle,
        arrow_length_ratio=0.12,
    )

    # ラベルを矢印の中ほどに配置
    label_position = start + 0.55 * vector
    ax.text(
        label_position[0],
        label_position[1],
        label_position[2],
        label,
        color=color,
        fontsize=11,
    )


# 入力ベクトル v
draw_vector(
    start=(0, 0, 0),
    vector=v,
    color="C0",
    label=r"$\vec{v}$",
    linewidth=2.5,
)

# 寄与ベクトル c_1, c_2, c_3 を先端から順につなぐ
start = np.zeros(3)
contribution_colors = ["C1", "C2", "C3"]

for j, contribution in enumerate(contributions):
    draw_vector(
        start=start,
        vector=contribution,
        color=contribution_colors[j],
        label=rf"$\vec{{c}}_{j + 1}$",
    )
    start += contribution

# 出力ベクトル w を原点から破線で描く
draw_vector(
    start=(0, 0, 0),
    vector=w,
    color="C4",
    label=r"$\vec{w}$",
    linestyle="--",
    linewidth=2.5,
)

# 入力条件と計算結果を図中に表示
ax.text2D(
    0.03,
    0.97,
    "D = [[1, 2, 0], [0, -1, 1], [2, 0, 1]]\n"
    "v = (1, 2, -1)\n"
    "w = (5, -3, 1)",
    transform=ax.transAxes,
    va="top",
    fontsize=10,
)

ax.set_title("Dyadic operation in 3D", fontsize=14)
ax.set_xlabel("x", fontstyle="italic")
ax.set_ylabel("y", fontstyle="italic")
ax.set_zlabel("z", fontstyle="italic")

# 入力・出力・寄与ベクトルがすべて見える範囲に設定
ax.set_xlim(-1, 6)
ax.set_ylim(-4, 3)
ax.set_zlim(-2, 4)
ax.set_box_aspect((7, 7, 6))
ax.view_init(elev=22, azim=-55)

ax.tick_params(axis="both", which="both", direction="inout")
ax.grid(True, linestyle=":", linewidth=0.7, alpha=0.6)

fig.tight_layout()
plt.show()

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

import numpy as np

import matplotlib.pyplot as plt

# Times New Romanを基本フォントに設定

plt.rcParams["font.family"] = "Times New Roman"

# 標準基底ベクトル e_1, e_2, e_3

basis = np.eye(3)

# ダイアディックの係数 D_ij

D = np.array([

[1.0, 2.0, 0.0],

[0.0, -1.0, 1.0],

[2.0, 0.0, 1.0],

])

# 入力ベクトル

v = np.array([1.0, 2.0, -1.0])

# 定義式の各項 D_ij e_i e_j を

# 標準基底上の行列として足し合わせる

D_dyadic = np.zeros((3, 3))

for i in range(3):

for j in range(3):

D_dyadic += D[i, j] * np.outer(basis[i], basis[j])

# ダイアディックを入力ベクトルに作用させる

w = D_dyadic @ v

# 入力の各成分から出力への寄与ベクトル c_j

contributions = np.array([

D[:, j] * v[j]

for j in range(3)

])

print("入力ベクトル v:", v)

print("出力ベクトル w:", w)

print("各寄与ベクトル c_j:\n", contributions)

fig = plt.figure(figsize=(8, 7))

ax = fig.add_subplot(111, projection="3d")

def draw_vector(start, vector, color, label, linestyle="-", linewidth=2):

"""startを始点とする3次元ベクトルを描く。"""

start = np.asarray(start, dtype=float)

vector = np.asarray(vector, dtype=float)

ax.quiver(

start[0], start[1], start[2],

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

color=color,

linewidth=linewidth,

linestyle=linestyle,

arrow_length_ratio=0.12,

)

# ラベルを矢印の中ほどに配置

label_position = start + 0.55 * vector

ax.text(

label_position[0],

label_position[1],

label_position[2],

label,

color=color,

fontsize=11,

)

# 入力ベクトル v

draw_vector(

start=(0, 0, 0),

vector=v,

color="C0",

label=r"$\vec{v}$",

linewidth=2.5,

)

# 寄与ベクトル c_1, c_2, c_3 を先端から順につなぐ

start = np.zeros(3)

contribution_colors = ["C1", "C2", "C3"]

for j, contribution in enumerate(contributions):

draw_vector(

start=start,

vector=contribution,

color=contribution_colors[j],

label=rf"$\vec{{c}}_{j + 1}$",

)

start += contribution

# 出力ベクトル w を原点から破線で描く

draw_vector(

start=(0, 0, 0),

vector=w,

color="C4",

label=r"$\vec{w}$",

linestyle="--",

linewidth=2.5,

)

# 入力条件と計算結果を図中に表示

ax.text2D(

0.03,

0.97,

"D = [[1, 2, 0], [0, -1, 1], [2, 0, 1]]\n"

"v = (1, 2, -1)\n"

"w = (5, -3, 1)",

transform=ax.transAxes,

va="top",

fontsize=10,

)

ax.set_title("Dyadic operation in 3D", fontsize=14)

ax.set_xlabel("x", fontstyle="italic")

ax.set_ylabel("y", fontstyle="italic")

ax.set_zlabel("z", fontstyle="italic")

# 入力・出力・寄与ベクトルがすべて見える範囲に設定

ax.set_xlim(-1, 6)

ax.set_ylim(-4, 3)

ax.set_zlim(-2, 4)

ax.set_box_aspect((7, 7, 6))

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

ax.tick_params(axis="both", which="both", direction="inout")

ax.grid(True, linestyle=":", linewidth=0.7, alpha=0.6)

fig.tight_layout()

plt.show()

　このコードを実行すると、次の図が表示されます。Matplotlibのインタラクティブプロットなので、自分の環境で是非ぐりぐりと回してみてください。

係数行列と入力ベクトルを設定する

　このコードでは、まずダイアディックの係数行列 $\mathbf D$ と入力ベクトル $\vec v$ を設定し、ダイアディックを入力ベクトルに作用させた結果 $\vec w$ を計算しています。その後、入力ベクトルの各成分が出力に与える寄与を求め、Matplotlibで3次元表示します。

D = np.array([
    [1.0,  2.0, 0.0],
    [0.0, -1.0, 1.0],
    [2.0,  0.0, 1.0],
])

v = np.array([1.0, 2.0, -1.0])

D = np.array([

[1.0, 2.0, 0.0],

[0.0, -1.0, 1.0],

[2.0, 0.0, 1.0],

])

v = np.array([1.0, 2.0, -1.0])

ここでは、係数行列 $\mathbf D$ と入力ベクトル $\vec v$ を設定しています。今回使う数値は、ダイアディックの計算と図の見方を確認するための例です。標準基底で表したダイアディックの成分行列と入力ベクトルの成分を掛けると、

\[ \vec w=\mathbf D\vec v = \begin{pmatrix} 1&2&0\\ 0&-1&1\\ 2&0&1 \end{pmatrix} \begin{pmatrix} 1\\2\\-1 \end{pmatrix} = \begin{pmatrix} 5\\-3\\1 \end{pmatrix} \]

となります。

定義式からダイアディックを組み立てる

basis = np.eye(3)
D_dyadic = np.zeros((3, 3))

for i in range(3):
    for j in range(3):
        D_dyadic += D[i, j] * np.outer(basis[i], basis[j])

basis = np.eye(3)

D_dyadic = np.zeros((3, 3))

for i in range(3):

for j in range(3):

D_dyadic += D[i, j] * np.outer(basis[i], basis[j])

basisには、3次元の標準基底ベクトル $\vec e_1,\vec e_2,\vec e_3$ を格納しています。二重ループの中では、定義式

\[ \vec{\mathbf D} = \sum_{i=1}^{3}\sum_{j=1}^{3} D_{ij}\,\vec e_i\vec e_j \]

の各項を一つずつ計算して加えています。 np.outer(basis[i], basis[j]) は、ダイアディック積 $\vec e_i\vec e_j$ を標準基底上の数値行列として表す処理です。

　この例では標準基底を使っているため、組み立てた D_dyadicは係数行列 Dと同じ数値になります。ただし、 D_dyadic は定義式から組み立てたダイアディックの行列表現、 D はその係数を並べた行列として、コード上では役割を分けています。

ダイアディックを入力ベクトルに作用させる

w = D_dyadic @ v

1	w = D_dyadic @ v

　ここで、組み立てたダイアディックを入力ベクトル $\vec v$ に作用させ、出力ベクトル $\vec w$ を求めています。計算結果は

\[ \vec w= \begin{pmatrix} 5\\-3\\1 \end{pmatrix} \]

です。 @ はNumPyで行列とベクトルを掛ける演算子です。

入力成分ごとの寄与を計算する

contributions = np.array([
    D[:, j] * v[j]
    for j in range(3)
])

contributions = np.array([

D[:, j] * v[j]

for j in range(3)

])

　 D[:, j] は係数行列の第 $j$ 列です。その列に入力ベクトルの第 $j$ 成分 $v_j$ を掛け、入力の第 $j$ 成分が出力に与える寄与ベクトル $\vec c_j$ を計算しています。

今回の各寄与は、

\[ \vec c_1= \begin{pmatrix}1\\0\\2\end{pmatrix}, \qquad \vec c_2= \begin{pmatrix}4\\-2\\0\end{pmatrix}, \qquad \vec c_3= \begin{pmatrix}0\\-1\\-1\end{pmatrix} \]

です。これらを足すと、

\[ \vec c_1+\vec c_2+\vec c_3 = \begin{pmatrix}5\\-3\\1\end{pmatrix} =\vec w \]

となります。つまり、寄与ベクトルは出力ベクトルがどの入力成分から作られたのかを確認するための分解です。

ダイアディックを計算するメリット

　ダイアディックを使うと、入力ベクトルの各方向成分が、出力の各方向にどう影響するかを一つの対象にまとめて表せます。今回の例では、出力ベクトルを求めるだけでなく、各入力成分からの寄与に分けて確認できました。

　そのため、ダイアディックの計算からは、入力ベクトルがどのように変換されたかに加えて、どの入力成分が出力のどの部分に寄与したかも読み取れます。係数を含む成分行列だけを見ると式としてまとまりますが、寄与ベクトルに分けて図示すると、その計算過程を視覚的に追えます。

※なお、この例の係数行列は計算の仕組みを示すために選んだ値です。実際の人工衛星の物理量を表す場合には、対象とする現象に応じて係数の意味や単位を定めます。

参考文献

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

GitHub

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

ためになるJAXAの資料を読んでスキルアップをします！
Xでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！

人工衛星の力学と制御に必要な数学の基礎：ダイアディック

ダイアディックとは？

2階のテンソル量

慣性テンソルに見る2階テンソルの働き

この記事で学ぶこと

ベクトルを別のベクトルに変換したいとき

例題：ベクトルから特定方向の成分を取り出す

Pythonで射影を確認

射影の考え方を一般化する

入力ベクトルの成分から出力ベクトルを作る

基底ベクトルを使って書く

3次元のベクトルに拡張する

ダイアディック

ダイアディックの定義

Pythonでダイアディック計算をする

係数行列と入力ベクトルを設定する

定義式からダイアディックを組み立てる

ダイアディックを入力ベクトルに作用させる

入力成分ごとの寄与を計算する

ダイアディックを計算するメリット

参考文献

GitHub

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

最近の投稿

アーカイブ

カテゴリー

ダイアディックとは？

2階のテンソル量

慣性テンソルに見る2階テンソルの働き

この記事で学ぶこと

ベクトルを別のベクトルに変換したいとき

例題：ベクトルから特定方向の成分を取り出す

Pythonで射影を確認

射影の考え方を一般化する

入力ベクトルの成分から出力ベクトルを作る

基底ベクトルを使って書く

3次元のベクトルに拡張する

ダイアディック

ダイアディックの定義

Pythonでダイアディック計算をする

係数行列と入力ベクトルを設定する

定義式からダイアディックを組み立てる

ダイアディックを入力ベクトルに作用させる

入力成分ごとの寄与を計算する

ダイアディックを計算するメリット

参考文献

GitHub

SNSでもご購読できます。

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

最近の投稿

アーカイブ

カテゴリー

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