多変数ニュートン法

前節では1変数の非線形方程式を扱いましたが、物理シミュレーションでは複数の変数が相互に依存する連立方程式を解く必要がしばしば生じます。

これを解くための標準的な手法が、ニュートン法を多変数に拡張した多変数ニュートン法 (Multivariable Newton’s Method) です。

アルゴリズム

1変数の場合の更新式をベクトルと行列に拡張します。

多変数のテイラー展開（1次近似）は以下のようになります。

ここではヤコビ行列 (Jacobian Matrix) です。

実用上は、このヤコビ行列の扱いが大きな問題になります。未知数が個ならは行列なので、要素数は個あります。各要素は偏微分であり、解析的に導出する場合も、差分で近似する場合も、自動微分を使う場合も、「残差ベクトルを評価する」だけの場合より設計が難しくなります。

したがって、修正量は以下の連立一次方程式を解くことで求められます。

これを解いてを更新します。

Rustによる実装

円と放物線の交点を求めてみましょう。

連立方程式は以下のようになります。

ヤコビ行列は：

計算の具体例 (1ステップ目)

初期値をとして、1回目の更新を手計算で追ってみましょう。

残差ベクトルの計算

ヤコビ行列の計算

修正量の計算 連立方程式を解きます。

これを解くと（クラメルの公式やガウス消去法などで）、

解の更新

次の推定値はとなります。これを繰り返すことで、真の解（この場合は）に近づいていきます。

実装には ndarray と ndarray-linalg を使用します。

Note

実行には Cargo.toml に ndarray と ndarray-linalg、およびバックエンド（例: openblas-src）の依存関係が必要です。

use ndarray::{arr1, arr2, Array1, Array2};
use ndarray_linalg::Solve;

fn residual(x: &Array1<f64>) -> Array1<f64> {
    let curr_x: f64 = x[0];
    let curr_y: f64 = x[1];

    arr1(&[
        curr_x.powi(2) + curr_y.powi(2) - 1.0,
        curr_y - curr_x.powi(2),
    ])
}

fn jacobian(x: &Array1<f64>) -> Array2<f64> {
    let curr_x: f64 = x[0];
    let curr_y: f64 = x[1];

    arr2(&[
        [2.0 * curr_x, 2.0 * curr_y],
        [-2.0 * curr_x, 1.0       ]
    ])
}

fn residual_norm(f_vec: &Array1<f64>) -> f64 {
    f_vec.iter().map(|v| v.powi(2)).sum::<f64>().sqrt()
}

fn newton_solve(mut x: Array1<f64>, tolerance: f64, max_iter: usize) -> Option<(Array1<f64>, usize)> {
    for i in 0..max_iter {
        // 残差ベクトル F(x)
        let f_vec = residual(&x);

        // 収束判定 (ノルムが十分小さいか)
        if residual_norm(&f_vec) < tolerance {
            return Some((x, i));
        }

        // ヤコビ行列 J(x)
        let j = jacobian(&x);

        // 連立一次方程式 J * delta_x = -F を解く
        // solve() は ndarray-linalg の機能
        let delta = j.solve(&(-f_vec)).expect("Singular Jacobian");

        // 更新
        x = x + delta;
    }

    None
}

fn main() {
    // 初期値 (x, y) = (1.0, 2.0)
    // 解に近い適切な初期値を選ぶ必要があります
    let x0 = arr1(&[1.0, 2.0]);

    let tolerance = 1e-8;
    let max_iter = 100;

    match newton_solve(x0, tolerance, max_iter) {
        Some((x, iter)) => println!("解が見つかりました: x={:.6}, y={:.6} (反復: {})", x[0], x[1], iter),
        None => println!("収束しませんでした"),
    }
}

解説

ループ: 1変数のときと同様に、収束するまでループします。
ヤコビ行列の評価: jacobian(&x) は今回のようなの問題では簡単ですが、未知数が多い問題では個の偏微分を評価・保存する必要があります。PDEを離散化して得られる問題ではヤコビ行列が疎になることが多く、その構造を使わずに密行列として扱うとすぐに破綻します。
連立一次方程式: j.solve(&(-f_vec)) の部分で、線形代数のソルバーが活躍します。密行列として解くと一般にの計算量が必要です。大規模な物理シミュレーションでは、ヤコビ行列の疎性を使うソルバーや、ヤコビ行列を明示的に作らない Newton-Krylov 法などを検討します。
逆行列: 数式上はと書けますが、数値計算では逆行列を明示的に求めてはいけません。連立方程式を解く方が、計算量的にも精度的にも有利だからです。

応用例

多変数ニュートン法は以下のような場面で使用されます。

陰的解法 (Implicit Methods): 時間発展の方程式（微分方程式）を解く際、次の時刻の状態を求めるために非線形方程式を解く必要がある場合（後退オイラー法など）。
構造解析: 大きな変形を伴う物体の釣り合い位置を求める問題。

まとめ

多変数ニュートン法は、1変数のニュートン法をベクトルと行列に一般化した手法である。
各ステップでヤコビ行列 を計算し、連立一次方程式を解くことで修正量を求める。
未知数が個ならヤコビ行列は個の要素を持つため、大規模問題では評価方法、疎性、保存形式を含めて設計する必要がある。
逆行列を計算するのではなく、LU分解などの線形ソルバーを用いるのが数値計算の鉄則である。
この手法は、非線形微分方程式の陰的解法など、物理シミュレーションの多くの場面で利用される。

次節では、方程式を解くことと密接に関連する「最適化問題」について扱います。

Keyboard shortcuts

Rust計算物理学チュートリアル

多変数ニュートン法

アルゴリズム

Rustによる実装

計算の具体例 (1ステップ目)

解説

応用例

まとめ