はじめに

解析力学において，系の運動方程式はラグランジアン $L = T - V$（運動エネルギー $T$ と位置エネルギー $V$ の差）から，オイラー=ラグランジュ方程式

$$
\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}}\right) - \frac{\partial L}{\partial q} = 0
$$

を通じて導出される．単振り子のように系が単純な場合，この方程式は解析的に整理でき，角加速度 $\ddot{\theta}$ を閉じた形で得ることができる．

しかし，より複雑な系（多自由度・非線形拘束・速度に依存する散逸項を含む系など）では，オイラー=ラグランジュ方程式を $\ddot{q}$ について解析的に解くことが困難な場合がある．本稿では，そのような場合にも適用可能な一般的アプローチとして，オイラー=ラグランジュ方程式の残差を定義し，これを最急降下法（勾配降下法）によってゼロに近づけることで角加速度を数値的に求める手法を示す．題材として単振り子を取り上げ，解析解と比較することで手法の妥当性を検証する．

1. 理論

1.1 単振り子のラグランジアン

質量 $m$，長さ $l$ の単振り子を考える．鉛直下方を基準とした振れ角を $\theta$，重力加速度を $g$ とすると，運動エネルギー $T$ と位置エネルギー $V$ は

$$
T = \frac{1}{2} m l^2 \dot{\theta}^2, \qquad V = -m g l \cos\theta
$$

と書ける．したがってラグランジアンは

$$
L(\theta, \dot{\theta}) = T - V = \frac{1}{2} m l^2 \dot{\theta}^2 + m g l \cos\theta
$$

である．

1.2 オイラー=ラグランジュ方程式と解析解

上記の $L$ をオイラー=ラグランジュ方程式に代入する．

$$
\frac{\partial L}{\partial \dot{\theta}} = m l^2 \dot{\theta}, \qquad
\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\theta}}\right) = m l^2 \ddot{\theta}, \qquad
\frac{\partial L}{\partial \theta} = -m g l \sin\theta
$$

これらを方程式に代入すると，

$$
m l^2 \ddot{\theta} + m g l \sin\theta = 0
\quad\Longrightarrow\quad
\boxed{\ddot{\theta} = -\frac{g}{l}\sin\theta}
$$

という周知の単振り子の運動方程式が得られる．本稿では，この解析解を数値解法の検証用の真値として用いる．

1.3 残差の定式化

本稿の主眼は，上記のような解析的な整理を行わずに，オイラー=ラグランジュ方程式をそのまま満たす $\ddot{\theta}$ を数値的に探索する点にある．そこで，時刻 $t$ における状態 $(\theta, \dot{\theta})$ を既知として，試行角加速度 $\ddot{\theta}$ の関数として残差を次のように定義する．

$$
R(\ddot{\theta}) = \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\theta}}\right) - \frac{\partial L}{\partial \theta}
$$

真の運動方程式を満たす $\ddot{\theta}$ に対しては $R(\ddot{\theta}) = 0$ となる．したがって，運動方程式を解く問題は

\ddot{\theta}^{*} = \operatorname*{arg\,min}_{\ddot{\theta}} \; R(\ddot{\theta})^2

という最小化問題に帰着できる．

1.4 数値微分による偏導関数の近似

$\partial L / \partial \theta$ および $\partial L / \partial \dot{\theta}$ は，微小量 $h$ を用いた中心差分によって近似する．

$$
\frac{\partial L}{\partial \theta} \approx \frac{L(\theta+h,\dot{\theta}) - L(\theta-h,\dot{\theta})}{2h}, \qquad
\frac{\partial L}{\partial \dot{\theta}} \approx \frac{L(\theta,\dot{\theta}+h) - L(\theta,\dot{\theta}-h)}{2h}
$$

一般化運動量 $p \equiv \partial L / \partial \dot{\theta}$ の時間微分 $dp/dt$ は，試行角加速度 $\ddot{\theta}$ を用いて1ステップ先（$\Delta t$ 後）の状態

$$
\theta_{\text{next}} = \theta + \dot{\theta},\Delta t, \qquad
\dot{\theta}_{\text{next}} = \dot{\theta} + \ddot{\theta},\Delta t
$$

を予測し，運動量の前進差分

$$
\frac{dp}{dt} \approx \frac{p(\theta_{\text{next}}, \dot{\theta}_{\text{next}}) - p(\theta, \dot{\theta})}{\Delta t}
$$

によって近似する．これにより，残差 $R(\ddot{\theta})$ は数値微分のみを用いて計算可能となる．

1.5 最急降下法による解法

残差の2乗 $R(\ddot{\theta})^2$ を目的関数とし，これを最急降下法で最小化する．目的関数の勾配は連鎖律より

\frac{d}{d\ddot{\theta}}\bigl[R(\ddot{\theta})^2\bigr] = 2\,R(\ddot{\theta})\,\frac{dR}{d\ddot{\theta}}

であり，$dR/d\ddot{\theta}$ もまた数値微分（前進差分）で近似する．学習率を $\alpha$ とすると，更新則は

\ddot{\theta} \leftarrow \ddot{\theta} - \alpha \cdot 2\,R(\ddot{\theta})\,\frac{dR}{d\ddot{\theta}}

となる．なお，単振り子の場合 $R(\ddot{\theta}) = m l^2,\ddot{\theta} + m g l \sin\theta$ は $\ddot{\theta}$ について線形であるため，誤差 $x = \ddot{\theta} - \ddot{\theta}^{*}$ は

x_{k+1} = (1 - 2\alpha\, m l^2)\, x_k

という単純な等比数列に従う．本稿の実装（$m=l=1$）では $|1-2\alpha| < 1$，すなわち $0 < \alpha < 1$ を満たす限り収束するが，収束速度は $\alpha$ に強く依存する．後述するように，$\alpha$ を小さく取りすぎると既定の反復回数内で収束せず，見かけ上「解けたことにされた」不正確な近似解が時間発展に用いられてしまう点に注意が必要である．

1.6 時間発展：シンプレクティック・オイラー法

各時刻で $\ddot{\theta}$ が求まれば，状態 $(\theta, \dot{\theta})$ を時間発展させて軌道を得る．単純な陽的オイラー法

\theta \leftarrow \theta + \dot{\theta}\,\Delta t, \qquad \dot{\theta} \leftarrow \dot{\theta} + \ddot{\theta}\,\Delta t

を用いると，振動系ではエネルギーが単調に増大し，振幅が時間とともに発散的に増加してしまう．そこで本稿では，速度を先に更新してから位置の更新に反映するシンプレクティック・オイラー法

\dot{\theta} \leftarrow \dot{\theta} + \ddot{\theta}\,\Delta t, \qquad \theta \leftarrow \theta + \dot{\theta}\,\Delta t

を採用する．この方法はエネルギーを長時間にわたり有界に保つ性質（シンプレクティック性）を持ち，振動系のシミュレーションに適している．

2. 数値実装

2.1 ラグランジアンと偏導関数

def lagrangian(theta, d_theta):
    """L = T - V"""
    T = 0.5 * m * l**2 * d_theta**2
    V = -m * g * l * np.cos(theta)
    return T - V


def dL_dtheta(theta, d_theta):
    """∂L/∂θ を中心差分で計算"""
    return (lagrangian(theta + h, d_theta) - lagrangian(theta - h, d_theta)) / (2 * h)


def dL_ddtheta(theta, d_theta):
    """∂L/∂θ̇ を中心差分で計算（一般化運動量 p）"""
    return (lagrangian(theta, d_theta + h) - lagrangian(theta, d_theta - h)) / (2 * h)

2.2 残差と最急降下法

def residual(dd_theta, theta, d_theta, dt):
    theta_next = theta + d_theta * dt
    d_theta_next = d_theta + dd_theta * dt

    p_now = dL_ddtheta(theta, d_theta)
    p_next = dL_ddtheta(theta_next, d_theta_next)
    dp_dt = (p_next - p_now) / dt

    return dp_dt - dL_dtheta(theta, d_theta)


def solve_dd_theta(theta, d_theta, dt, dd_theta0=0.0):
    dd_theta = dd_theta0
    for _ in range(max_iter):
        f = residual(dd_theta, theta, d_theta, dt)
        if f**2 < tol:
            break
        f_plus = residual(dd_theta + h, theta, d_theta, dt)
        df = (f_plus - f) / h          # dR/d(θ̈)
        grad = 2 * f * df              # d(R^2)/d(θ̈)
        dd_theta -= alpha * grad
    return dd_theta

各時刻において solve_dd_theta を呼び出す際，直前の時刻で得られた $\ddot{\theta}$ を初期値として与える（ウォームスタート）．状態 $(\theta,\dot{\theta})$ は時間刻み $\Delta t$ ごとにわずかしか変化しないため，これにより収束までの反復回数を大幅に削減できる．

2.3 時間発展ループ

def simulate(theta0, d_theta0, t_end, dt):
    n = int(t_end / dt)
    ts = np.zeros(n)
    thetas = np.zeros(n)
    d_thetas = np.zeros(n)

    theta, d_theta = theta0, d_theta0
    dd_theta_guess = 0.0

    for i in range(n):
        ts[i] = i * dt
        thetas[i] = theta
        d_thetas[i] = d_theta

        dd_theta = solve_dd_theta(theta, d_theta, dt, dd_theta_guess)
        dd_theta_guess = dd_theta

        d_theta = d_theta + dd_theta * dt
        theta = theta + d_theta * dt

    return ts, thetas, d_thetas

3. 実験結果と考察

3.1 学習率の選定

最急降下法の学習率 $\alpha$ は，収束の速度と安定性を左右する重要なハイパーパラメータである．本実装の初期検討では $\alpha = 10^{-3}$，最大反復回数 $500$ で試験したところ，残差は既定回数内でゼロに収束せず，得られた $\ddot{\theta}$ は解析解の約 $63.2%$ にとどまるという系統的な誤差が観測された．

1.5節で述べた線形収束モデル $x_{k+1} = (1-2\alpha),x_k$（$m=l=1$）に基づけば，初期値を $\ddot{\theta}_0 = 0$ とした場合の $n$ 回反復後の誤差比は $(1-2\alpha)^n$ で与えられる．$\alpha=10^{-3}$，$n=500$ を代入すると

(1-2\times10^{-3})^{500} \approx e^{-1} \approx 0.368

となり，解に対する残存誤差比は理論値と実測値（$1 - 0.632 = 0.368$）が良く一致する．この結果は，最急降下法を用いた数値解法において，学習率と反復回数の設計を誤ると，収束前の値が「見かけ上の解」として時間発展に使われてしまう危険性を定量的に裏付けている．

そこで本実装では $\alpha = 0.3$，最大反復回数 $2000$，収束判定 $R^2 < 10^{-14}$ を採用した．この設定では，代表的な $\theta$ の値に対して以下のように，20回程度の反復で解析解に一致する結果が得られた．

$\theta$ [rad]	数値解 $\ddot{\theta}$	解析解 $-,(g/l)\sin\theta$	反復回数
0.1	-0.978367	-0.978367	18
0.5	-4.698370	-4.698370	20
1.0	-8.246416	-8.246416	21
1.5	-9.775451	-9.775451	21

3.2 時間発展のシミュレーション

初期条件 $\theta_0 = \pi/2$，$\dot{\theta}_0 = 0$，時間刻み $\Delta t = 0.01$，計算時間 $10,\mathrm{s}$ でシミュレーションを実行した結果を図1に示す。

図1：最急降下法による単振り子の数値解。上段：角度 $\theta(t)$，下段：角速度 $\dot{\theta}(t)$。

得られた $\theta(t)$ は振幅 $|\theta| \approx \pi/2$ を保ったまま安定に振動しており，初期条件から決まるエネルギー保存則と整合する結果となった。実際に力学的エネルギー

E = \frac{1}{2}ml^2\dot\theta^2 - mgl\cos\theta

を全時間ステップにわたって評価したところ，理論初期値 $E_0 \approx 0$ に対し，シミュレーション全体を通じて $E \in [-0.133,\ 0.136]$ の範囲に収まっており，時間とともに単調増加も単調減少もしないことを確認した。この有界な振動は，1.6節で述べたシンプレクティック・オイラー法に特有の挙動であり，長時間シミュレーションにおいてもエネルギーが発散しないことを示している。

比較として，シンプレクティック化を行わない単純な陽的オイラー法（同一の $\Delta t=0.01$，$t=10,\mathrm{s}$）では，振幅が $\theta_0=\pi/2,(\approx1.571)$ から $|\theta|\approx 2.88$ まで増大し，エネルギーも系統的に増加することを別途確認している。この比較からも，最急降下法によって求めた $\ddot{\theta}$ の精度自体は妥当であり，時間発展の積分スキームの選択が長時間安定性を左右する主要因であることがわかる。

4. まとめ

本稿では，単振り子のラグランジュの運動方程式を題材に，オイラー=ラグランジュ方程式の残差を最急降下法によって最小化することで角加速度を数値的に求める手法を示した。得られた知見を以下に整理する。

オイラー=ラグランジュ方程式は，残差 $R(\ddot{\theta})$ の2乗を目的関数とする最小化問題として定式化でき，解析的に $\ddot{q}$ について解けない系にも適用可能な一般的枠組みを与える。
最急降下法の学習率と反復回数の設定は収束精度に直結し，過小な学習率は既定反復回数内での未収束という見えにくい誤差を生む。本稿では線形収束モデルによってこの誤差を定量的に説明した。
求めた $\ddot{\theta}$ を用いた時間発展には積分スキームの選択が重要であり，陽的オイラー法はエネルギーの増大を招く一方，シンプレクティック・オイラー法はエネルギーを有界に保ち，長時間の振動シミュレーションに適する。
単振り子について得られた数値解は，解析解 $\ddot{\theta} = -(g/l)\sin\theta$ および理論的な振幅・エネルギー保存と良く一致し，手法の妥当性が確認された。

本手法は解析的取り扱いが困難な多自由度系や非線形拘束系への拡張が可能であり，今後の課題としたい。

参考文献

江沢洋，解析力学（物理学One Point），共立出版，1990．
H. Goldstein, C. Poole, J. Safko, Classical Mechanics, 3rd ed., Addison-Wesley, 2002.
森正武，数値解析（第2版），共立出版，2002．

付録：プログラム全文

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib

# ============================================================
# 単振り子のラグランジュの運動方程式を、各時刻ごとに
# オイラー=ラグランジュ方程式の残差を最急降下法で最小化して
# 角加速度を求める（数値的に）ことで解くプログラム
# ============================================================

# 物理パラメータ
m = 1.0      # 質点の質量 [kg]
l = 1.0      # 振り子の長さ [m]
g = 9.8      # 重力加速度 [m/s^2]

# 数値微分用の微小量（積分刻みdtとは独立に小さく取る）
h = 1e-6

# 最急降下法のパラメータ
alpha = 0.3        # 学習率
tol = 1e-14        # 収束判定（残差の2乗）
max_iter = 2000    # 最大反復回数


def lagrangian(theta, d_theta):
    """L = T - V"""
    T = 0.5 * m * l**2 * d_theta**2
    V = -m * g * l * np.cos(theta)
    return T - V


def dL_dtheta(theta, d_theta):
    """∂L/∂θ を中心差分で計算"""
    return (lagrangian(theta + h, d_theta) - lagrangian(theta - h, d_theta)) / (2 * h)


def dL_ddtheta(theta, d_theta):
    """∂L/∂θ̇ を中心差分で計算（一般化運動量 p）"""
    return (lagrangian(theta, d_theta + h) - lagrangian(theta, d_theta - h)) / (2 * h)


def residual(dd_theta, theta, d_theta, dt):
    """
    オイラー=ラグランジュ方程式の残差
        R(θ̈) = d/dt(∂L/∂θ̇) - ∂L/∂θ
    d/dt(∂L/∂θ̇) は θ̈ を使って1ステップ先の状態を予測し、
    運動量 p=∂L/∂θ̇ の前進差分で近似する。
    """
    theta_next = theta + d_theta * dt
    d_theta_next = d_theta + dd_theta * dt

    p_now = dL_ddtheta(theta, d_theta)
    p_next = dL_ddtheta(theta_next, d_theta_next)
    dp_dt = (p_next - p_now) / dt

    return dp_dt - dL_dtheta(theta, d_theta)


def solve_dd_theta(theta, d_theta, dt, dd_theta0=0.0):
    """
    最急降下法で R(θ̈)^2 を最小化し、θ̈ を求める。
    直前ステップの解を初期値にする（ウォームスタート）ことで収束を速める。
    """
    dd_theta = dd_theta0
    for _ in range(max_iter):
        f = residual(dd_theta, theta, d_theta, dt)
        if f**2 < tol:
            break
        f_plus = residual(dd_theta + h, theta, d_theta, dt)
        df = (f_plus - f) / h          # dR/d(θ̈)
        grad = 2 * f * df              # d(R^2)/d(θ̈)
        dd_theta -= alpha * grad
    return dd_theta


def simulate(theta0, d_theta0, t_end, dt):
    """θ(t), θ̇(t), θ̈(t) を時間発展させて求める"""
    n = int(t_end / dt)
    ts = np.zeros(n)
    thetas = np.zeros(n)
    d_thetas = np.zeros(n)

    theta, d_theta = theta0, d_theta0
    dd_theta_guess = 0.0

    for i in range(n):
        ts[i] = i * dt
        thetas[i] = theta
        d_thetas[i] = d_theta

        dd_theta = solve_dd_theta(theta, d_theta, dt, dd_theta_guess)
        dd_theta_guess = dd_theta  # 次ステップのウォームスタート

        # シンプレクティック・オイラー法で状態を更新
        # （速度を先に更新してから位置に反映することでエネルギーの
        #   増大を抑え、振動系でも安定した軌道が得られる）
        d_theta = d_theta + dd_theta * dt
        theta = theta + d_theta * dt

    return ts, thetas, d_thetas


if __name__ == "__main__":
    theta0 = np.pi / 2   # 初期角度
    d_theta0 = 0.0        # 初期角速度
    t_end = 10.0
    dt = 0.01

    ts, thetas, d_thetas = simulate(theta0, d_theta0, t_end, dt)

    fig, axes = plt.subplots(2, 1, figsize=(8, 6), sharex=True)

    axes[0].plot(ts, thetas)
    axes[0].set_ylabel("θ [rad]")
    axes[0].set_title("最急降下法によるラグランジュの運動方程式の数値解（単振り子）")
    axes[0].grid(True)

    axes[1].plot(ts, d_thetas, color="tab:orange")
    axes[1].set_xlabel("t [s]")
    axes[1].set_ylabel("dθ/dt [rad/s]")
    axes[1].grid(True)

    plt.tight_layout()
    plt.savefig("lagrange_result.png", dpi=150)
    plt.show()