変分問題を関数プログラミングしてみる

変分問題は 汎関数 （関数を変数とする関数）を対象に、その最小化や最大化などを扱います。

関数を変数とする関数、ということで JavaScript で関数プログラミングするとどんな感じになるか試してみました。

  • なお小生、いまだに JavaScript に不自由な方、です。。。（自慢すんな！）
  • 一部だけですが、恐る恐る ECMAScript 6 の機能を使ってみました

例題

ここでは古典力学のラグランジアンと最小作用の原理を例にします。また簡単のため１粒子の１次元運動のみを扱います。

ラグランジアン

エネルギーの次元をもつ以下の量を表す関数をラグランジアンといいます。ここで x は粒子の位置座標、 v は速度。

$$
\begin{eqnarray}
L(x, v) &=& \mathbf{運動エネルギー - ポテンシャルエネルギー} \\
&=& \frac{1}{2} M v^2 - U(x)
\end{eqnarray}
$$

例えば、質量 $M$ バネ定数 $K$ の振動子のラグランジアンは以下になります。

$$
L_{oscillator}(x, v) = \frac{1}{2} M v^2 - \frac{1}{2} K x^2
$$

作用積分

あるラグランジアンに対して、以下の形の時間積分を作用あるいは作用積分といいます。

$$
S[x(t)] = \int_{t0}^{t1} L(x(t), \dot{x}(t)) dt
\ \ ただし \ \dot{x}(t) = \dfrac{d x(t)}{dt}
$$

ここで、t0, t1 は任意に固定した時刻、$x(t)$ は粒子の運動を表す 任意の関数 です。
$S$ は任意の関数 $x(t)$ を変数とする関数、すなわち汎関数で、汎関数であることを示すために、$S[x(t)]$ と書きます。

最小作用の原理

上の作用積分 $S[x(t)]$ の変数 x(t) は任意の関数でよいのですが、
$S[x(t)]$ を最小にするような関数が 実際に起こる運動 となることが知られています。（正しくは「最小」でなく「停留」です）

このための必要条件は、関数 $x(t)$ を少し変えた関数 $x(t) + \delta x(t)$ を考えたとき、$S$ が変化しないこと、すなわち以下で表せます。

$$
\delta S[x(t)] = S[x(t) + \delta x(t)] - S[x(t)] = 0 \hspace{5em} (A)
$$

ただし 境界条件 として、積分の両端では値が変わらないとします。すなわち

$$\delta x(t0) = \delta x(t1) = 0 \hspace{5em} (B)
$$

プログラム

通常は作用積分の停留条件から得られるオイラー・ラグランジュの方程式（２階の微分方程式）を解いて実際の運動を求めるのですが、
ここでは上の作用積分 $S[x(t)]$ を数値的に求めるプログラムを作成し、$x(t)$ を変化させて $S[x(t)]$ の変化を調べてみます。

振動子のラグランジアン

以下は、質量 M, バネ定数 K の振動子のラグランジアン $L(x,v)$ を作成します。

1
function oscillator(M,K) {
    return function(x,v) {
	return 0.5*(M*v*v - K*x*x);
    };
}

数値微分演算子

以下はきざみ幅 dt で、与えられた関数を数値微分する微分演算子を作成します。

1
function D_t(dt) {
    return function(func_t) {
	return function(t) {
	    return (func_t(t + 0.5*dt) - func_t(t - 0.5*dt)) / dt;
	};
    };
}

D_t(dt)(func) は関数 func の微分を表す関数（すなわち導関数）を返します。
（これも関数を変数にとる関数だが、通常、値・数値を返す関数を 汎関数、
　別の関数を返す関数は 演算子,作用素 というようです）

テストをかねて、以下の f, df, ddf を google chart でプロットしてみます。

1
var d_t = D_t(0.0001); // 数値微分演算子
var f = Math.cos;      // 微分する関数
var df  = d_t(f);      // その導関数
var ddf = d_t(d_t(f)); // 2次の導関数

数値定積分演算子

以下は、区間 [t0, t1] を分割数 ndiv で数値積分する数値定積分演算子を作成します。

1
function Int_t(t0, t1, ndiv) {
    var dt = (t1 - t0)/ndiv;
    return function(func_t) {
	var t = t0 + 0.5*dt;
	var s = 0.0;
	for (var i = 0; i < ndiv; i++) {
            s += func_t(t);
	    t += dt;
	}
	return s * dt;
    };
}

作用積分

以下は、与えられたラグランジアンに対して、作用積分汎関数 $S[x(t)]$ を作成します。
int_t は上の定積分演算子、d_t は微分演算子です。

1
function ActionIntegral(int_t, d_t, lagrangian) {
    return function(xt) {
	var dxdt = d_t(xt);
        var func_t = function(t) {
	    return lagrangian(xt(t), dxdt(t));
	};
	return int_t(func_t);
    };
}

トライアル

準備ができたので、試行関数 $x(t)$ をいろいろ変えて $S[x(t)]$ の変化を調べてみます。
$x(t)$ は 任意の関数 ですが、もちろんそんなことはできないので、いくつかのパラメータで $x(t)$ を表し、パラメータを変化させてみることになります。

ここではごくごく単純に １つのパラメータ k を持つ以下を試行関数にして、以下の条件で調べてみます。

  • ラグランジアン：　上の振動子のものを使用
    （解は $x(t) = Asin(\omega t) + B cos(\omega t)$ ($\omega = \sqrt{K/M}, A, B は定数$）
  • 積分範囲: $[t_0, t_1]=[0,1]$
  • 境界条件：$x(t0)=0、x(t1)=1$
  • 試行関数：$x(t;k) = sin(kt)/sin(k t_1)$

また積分範囲は とし、境界条件は とします。

1
function trialFunc(t, w) {
    return Math.sin(w*t)/Math.sin(w*t1);
}

Comment

2016-02-18
▲