質点系の運動

物理学、特に古典力学におけるシミュレーションの基礎となるのは、ニュートンの運動方程式 (Newton’s Equation of Motion) です。本節では、運動方程式を数値計算可能な形、すなわち常微分方程式 (ODE) の形に定式化する方法を確認します。

運動方程式の定式化

質量の質点にかかる力を、位置をとすると、運動方程式は以下のように記述されます。

ここで、は速度です。常微分方程式で学んだ標準的な数値解法（オイラー法やルンゲ＝クッタ法）を適用するためには、この2階の微分方程式を1階の連立微分方程式に書き直す必要があります。

状態ベクトルを位置と速度を並べたものとして定義します。

このとき、状態ベレトルの時間微分は次のように得られます。

これを一般化した形式とみなせば、常微分方程式のアルゴリズムがそのまま利用可能になります。

自由落下と空気抵抗

例として、空気抵抗がある場合の自由落下を考えます。重力加速度を、空気抵抗が速度に比例すると仮定します（係数）。

運動方程式を1階の連立微分方程式に分解すると：

Rustによる実装例

ndarray入門で学んだndarrayを使用して、このダイナミクスを実装してみましょう。状態ベクトルの各要素をと対応させます。

use ndarray::{Array1, arr1};

struct FallingBody {
    m: f64, // 質量
    k: f64, // 空気抵抗係数
    g: f64, // 重力加速度
}

impl FallingBody {
    fn new(m: f64, k: f64) -> Self {
        Self { m, k, g: 9.8 }
    }

    /// ODEソルバーが期待する形式: f(t, y) -> dy/dt
    /// 状態ベクトル y = [x, y, vx, vy]
    fn dynamics(&self, _t: f64, y: &Array1<f64>) -> Array1<f64> {
        // 現在の速度を取り出す
        let vx = y[2];
        let vy = y[3];

        // 加速度の計算: a = F/m = (重力 + 空気抵抗) / m
        // ax = - (k/m) * vx
        let ax = -(self.k / self.m) * vx;
        // ay = -g - (k/m) * vy
        let ay = -self.g - (self.k / self.m) * vy;

        // dy/dt = [dx/dt, dy/dt, dvx/dt, dvy/dt] = [vx, vy, ax, ay] を返す
        arr1(&[vx, vy, ax, ay])
    }
}

fn main() {
    let model = FallingBody::new(1.0, 0.1);

    // 初期状態: 原点から初速 (10.0, 15.0) で投げ出された状態
    let y = arr1(&[0.0, 0.0, 10.0, 15.0]);

    // 時刻 t=0 における変化率 (微分値) を計算
    let dy_dt = model.dynamics(0.0, &y0);

    println!("現在の状態: {:?}", y0);
    println!("微分値 (dy/dt): {:?}", dy_dt);
    // これを Runge-Kutta ソルバー等に渡すことで
    // 時間発展をシミュレーションできる
}

自由落下のシミュレーション