線形代数学と工学

主なトピックと工学への応用例

1. 行列とベクトル
   • 行列とベクトルの基本演算、行列の積、逆行列、行列式など。
   • 応用例：電気回路の解析、機械の動作のモデリング、画像処理におけるフィルタリング。

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter
from skimage import data

# サンプル画像の読み込み (カメラマン画像)
image = data.camera()

# ガウシアンフィルタ適用
filtered_image = gaussian_filter(image, sigma=3)

# 元画像とフィルタリング後の画像を並べて表示
fig, ax = plt.subplots(1, 2, figsize=(12, 6))
ax[0].imshow(image, cmap='gray')
ax[0].set_title("Original Image")
ax[0].axis('off')

ax[1].imshow(filtered_image, cmap='gray')
ax[1].set_title("Gaussian Filtered Image")
ax[1].axis('off')

plt.show()

import numpy as np
import matplotlib.pyplot as plt

def gaussian_kernel(size: int, sigma: float) -> np.ndarray:
    """ガウシアンフィルタのカーネル行列を作成する関数。
    
    Args:
        size (int): カーネルのサイズ（奇数であることが望ましい）。
        sigma (float): ガウス関数の標準偏差。

    Returns:
        np.ndarray: ガウシアンカーネル行列。
    """
    # カーネルの中心を基準に座標系を設定
    ax = np.linspace(-(size // 2), size // 2, size)
    xx, yy = np.meshgrid(ax, ax)

    # ガウス関数を使用して各座標の値を計算
    kernel = np.exp(-0.5 * (np.square(xx) + np.square(yy)) / np.square(sigma))

    # 正規化（カーネル行列の値の総和を1にする）
    return kernel / np.sum(kernel)

# カーネル行列を作成
kernel_size = 5  # カーネルサイズ
sigma_value = 1.0  # 標準偏差

# ガウシアンカーネルの行列を作成
gaussian_kernel_matrix = gaussian_kernel(kernel_size, sigma_value)

# カーネル行列を表示
print("Gaussian Kernel Matrix:")
print(gaussian_kernel_matrix)

# カーネル行列の可視化
plt.imshow(gaussian_kernel_matrix, interpolation='nearest', cmap='viridis')
plt.colorbar()
plt.title(f"Gaussian Kernel (size={kernel_size}, sigma={sigma_value})")
plt.show()

1. 線形方程式系
   • 連立一次方程式の解法（ガウス消去法、クラメルの公式など）。
   • 応用例：構造工学における力の平衡解析、交通流解析、電力供給網の解析。

import numpy as np
import matplotlib.pyplot as plt

# 節点 (ノード) の位置 [x, y]
nodes = np.array([[0, 0], [2, 0], [1, 1.5]])

# 要素（部材）の接続（始点と終点のノード番号）
elements = np.array([[0, 1], [1, 2], [2, 0]])

# 各要素の剛性行列を設定（簡単のため、すべて同じ値に設定）
element_stiffness = 1.0

# グローバル剛性行列を作成
K = np.zeros((6, 6))  # 2D平面なので各ノードは2つの自由度を持つ (3ノード x 2自由度)

# 各要素の剛性行列をグローバル剛性行列に組み込む
for element in elements:
    node1, node2 = element
    x1, y1 = nodes[node1]
    x2, y2 = nodes[node2]
    length = np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
    c = (x2 - x1) / length
    s = (y2 - y1) / length
    k_local = element_stiffness * np.array([[ c**2,  c*s, -c**2, -c*s],
                                            [ c*s,  s**2, -c*s, -s**2],
                                            [-c**2, -c*s,  c**2,  c*s],
                                            [-c*s, -s**2,  c*s,  s**2]])

    # グローバル剛性行列への組み込み
    dof = [2*node1, 2*node1+1, 2*node2, 2*node2+1]
    for i in range(4):
        for j in range(4):
            K[dof[i], dof[j]] += k_local[i, j]

# 外力ベクトル（ノードごとの外力） [Fx, Fy] (例として、節点2に外力を加える)
F = np.array([0, 0, 0, 0, 0, -10])

# 支点条件（拘束条件）を設定 (節点0と節点1のx, y方向を固定)
fixed_dof = [0, 1, 2, 3]  # 拘束された自由度

# 剛性行列と外力ベクトルの削減
K_reduced = np.delete(np.delete(K, fixed_dof, axis=0), fixed_dof, axis=1)
F_reduced = np.delete(F, fixed_dof)

# 変位の計算
U_reduced = np.linalg.solve(K_reduced, F_reduced)

# 拘束条件を考慮した全変位ベクトルを組み立て
U = np.zeros(K.shape[0])
U[np.setdiff1d(np.arange(K.shape[0]), fixed_dof)] = U_reduced

# 変位の表示
print("Nodal Displacements (x, y):")
for i in range(len(U)//2):
    print(f"Node {i}: {U[2*i]:.4f}, {U[2*i+1]:.4f}")

# 交通流解析のための係数行列と外部交通流量のベクトルを設定
A = np.array([[ 1, -1,  0,  0],
              [ 0,  1, -1,  0],
              [ 0,  0,  1, -1],
              [ 1,  0,  0, -1]])

# 各リンクの交通流量（外部流入量、正の値は流入、負の値は流出）
b = np.array([100, 0, 0, -100])

# リンク内の交通流量を求める
x = np.linalg.solve(A, b)

# 各リンクの交通流量を表示
print("Traffic Flow in Each Link:")
for i, flow in enumerate(x):
    print(f"Link {i+1}: {flow} vehicles/hour")

# リンク間の交通流量を図示
plt.figure(figsize=(8, 4))
plt.bar(range(len(x)), x)
plt.xlabel("Link")
plt.ylabel("Traffic Flow (vehicles/hour)")
plt.title("Traffic Flow in Each Link")
plt.show()

1. 固有値・固有ベクトル
   • 固有値問題の定式化とその解法。
   • 応用例：振動系の解析（モード解析）、動力学系の安定性解析、画像認識における主成分分析（PCA）。

import numpy as np
import matplotlib.pyplot as plt

# 剛性行列 K と 質量行列 M の定義（例：2自由度系）
K = np.array([[4, -2],
              [-2, 2]])  # 剛性行列 (Stiffness matrix)
M = np.array([[2, 0],
              [0, 1]])  # 質量行列 (Mass matrix)

# 一般化固有値問題を解く (K*u = ω^2 * M*u)
eigvals, eigvecs = np.linalg.eig(np.linalg.inv(M) @ K)

# 固有振動数の取得
natural_frequencies = np.sqrt(eigvals)

# 固有ベクトルの正規化
eigvecs = eigvecs / np.linalg.norm(eigvecs, axis=0)

# 結果の表示
print("固有振動数 (Natural Frequencies):", natural_frequencies)
print("振動モード (Mode Shapes):\n", eigvecs)

# モードのプロット設定
modes = eigvecs.T

# プロット
fig, ax = plt.subplots(1, 2, figsize=(12, 5))

for i, mode in enumerate(modes):
    ax[i].plot([0, 1], [0, mode[0]], marker='o', label=f'Mass 1')
    ax[i].plot([0, 1], [0, mode[1]], marker='o', label=f'Mass 2')
    ax[i].axhline(0, color='black', linewidth=0.5)
    ax[i].set_title(f'Mode {i+1}: Frequency = {natural_frequencies[i]:.2f}')
    ax[i].set_xlabel('Degree of Freedom')
    ax[i].set_ylabel('Amplitude')
    ax[i].legend()
    ax[i].grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.datasets import load_digits

# データの読み込み（例: 手書き数字データ）
digits = load_digits()
X = digits.data

# データの平均をゼロにする
X_centered = X - np.mean(X, axis=0)

# PCA を適用
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_centered)

# 固有値と固有ベクトルの取得
eigenvalues = pca.explained_variance_
eigenvectors = pca.components_

# 結果のプロット
plt.figure(figsize=(8, 6))
plt.scatter(X_pca[:, 0], X_pca[:, 1], c=digits.target, cmap='viridis', edgecolor='k', s=50)
plt.colorbar()
plt.title("PCA of Digits Dataset")
plt.xlabel("Principal Component 1")
plt.ylabel("Principal Component 2")

# 固有ベクトルを矢印でプロット
for i, v in enumerate(eigenvectors):
    plt.arrow(0, 0, v[0]*5, v[1]*5, head_width=0.3, head_length=0.3, linewidth=2, color='red')
    plt.text(v[0]*5, v[1]*5, f"Eigenvector {i+1}", color='red')

plt.grid()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 状態行列 A を定義（例: 2次元系）
A = np.array([[0, 1],
              [-2, -3]])

# 固有値と固有ベクトルの計算
eigenvalues, eigenvectors = np.linalg.eig(A)

# 固有値の表示
print(f"Eigenvalues: {eigenvalues}")

# 固有ベクトルの表示
print(f"Eigenvectors:\n{eigenvectors}")

# 固有値とベクトルのプロット
plt.figure(figsize=(8, 8))

# x, y 軸の設定
plt.axhline(0, color='gray', linewidth=0.5)
plt.axvline(0, color='gray', linewidth=0.5)

# 固有ベクトルのプロット
origin = np.array([[0, 0], [0, 0]])  # 原点
plt.quiver(*origin, eigenvectors[0, :], eigenvectors[1, :], color=['r', 'b'], scale=5)

# 安定性の判定を可視化
for eig in eigenvalues:
    plt.scatter(np.real(eig), np.imag(eig), color='green', s=100, label=f'Eigenvalue: {eig:.2f}')
    
# グラフの設定
plt.title('Stability Analysis: Eigenvalues and Eigenvectors')
plt.xlabel('Real Part')
plt.ylabel('Imaginary Part')
plt.grid()
plt.axvline(0, color='black', linestyle='--')
plt.axhline(0, color='black', linestyle='--')

# 凡例を表示
plt.legend(loc='best')
plt.show()

1. ベクトル空間と基底
   • ベクトル空間、部分空間、基底と次元の概念。
   • 応用例：信号処理におけるフーリエ変換、圧縮センシング、ロボティクスにおける姿勢制御。

import numpy as np
import matplotlib.pyplot as plt

# ロボットアームのパラメータ設定
length = 1.0  # ロボットアームの長さ（メートル）
theta = np.pi / 4  # 初期角度（45度）

# 目標角度
target_angle = np.pi / 6  # 目標角度（30度）

# 制御パラメータ
Kp = 2.0  # 比例制御ゲイン（Pゲイン）

# 時間設定
dt = 0.01  # 時間刻み（秒）
time = np.arange(0, 5, dt)  # 0秒から5秒までの時間配列

# 角度の履歴を保存するリスト
theta_history = []

# 姿勢制御のシミュレーション
for t in time:
    # 現在の角度誤差
    error = target_angle - theta

    # 比例制御（P制御）による角度変化量
    delta_theta = Kp * error

    # 角度を更新
    theta += delta_theta * dt

    # 角度の履歴を保存
    theta_history.append(theta)

# シミュレーション結果をプロット
plt.figure(figsize=(12, 6))

# 角度の時間変化をプロット
plt.subplot(1, 2, 1)
plt.plot(time, np.degrees(theta_history), label='Current Angle')
plt.axhline(np.degrees(target_angle), color='r', linestyle='--', label='Target Angle')
plt.title("Angle Control Over Time")
plt.xlabel("Time [s]")
plt.ylabel("Angle [degrees]")
plt.legend()
plt.grid()

# 最終姿勢のプロット
plt.subplot(1, 2, 2)
x = [0, length * np.cos(theta_history[-1])]  # アームのx座標
y = [0, length * np.sin(theta_history[-1])]  # アームのy座標
plt.plot(x, y, marker='o')
plt.xlim(-length, length)
plt.ylim(-length, length)
plt.title("Final Robot Arm Position")
plt.xlabel("X position [m]")
plt.ylabel("Y position [m]")
plt.axhline(0, color='k', linestyle='--')
plt.axvline(0, color='k', linestyle='--')
plt.grid()

# グラフを表示
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linprog

# スパース信号の生成
n = 100  # 元信号の次元
m = 30   # 観測値の次元
k = 5    # スパース信号の非ゼロ要素数

# 元信号 x の生成（ランダムに k 個の非ゼロ要素を持つ）
x_true = np.zeros(n)
x_true[np.random.choice(n, k, replace=False)] = np.random.randn(k)

# 観測行列 A の生成（ランダムな正規分布行列）
A = np.random.randn(m, n)

# 観測値 y の生成
y = A.dot(x_true)

# 最小 L1 ノルムを解くための線形計画問題の設定
c = np.ones(2 * n)  # コスト関数 (L1 ノルム)
A_eq = np.hstack([A, -A])
b_eq = y

# 線形計画法で解を求める
result = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=(0, None), method='highs')

# 求めた解の信号
x_recovered = result.x[:n] - result.x[n:]

# プロット
plt.figure(figsize=(12, 6))

# 元信号
plt.subplot(1, 2, 1)
plt.stem(x_true, use_line_collection=True)
plt.title('Original Sparse Signal')
plt.xlabel('Index')
plt.ylabel('Amplitude')

# 復元信号
plt.subplot(1, 2, 2)
plt.stem(x_recovered, use_line_collection=True)
plt.title('Recovered Signal from Measurements')
plt.xlabel('Index')
plt.ylabel('Amplitude')

plt.tight_layout()
plt.show()

1. 線形写像と変換

   • 線形写像の定義、核と像、行列表示。
   • 応用例：ロボットアームの位置制御、画像のアフィン変換や投影変換。

2. 直交性と内積空間

   • ベクトルの直交性、内積、グラム・シュミットの正規直交化法。
   • 応用例：信号の直交分解（フーリエ級数展開）、機械学習のサポートベクトルマシン（SVM）。

1. 特異値分解 (SVD) と行列分解

   • 特異値分解、QR分解、LU分解の定義と計算方法。
   • 応用例：データ圧縮、レコメンダシステム、画像復元。

2. 数値的手法とアルゴリズム

   • 数値解法、反復法（共役勾配法など）、条件数と数値安定性の概念。
   • 応用例：大規模システムの数値解析、有限要素法（FEM）による構造解析。