概要

「ラグランジュ形式」や「ハミルトン形式」で解説したような、 座標系の取り方によらない運動方程式に関して、 曲面上に拘束された物体の運動について考えてみます。

また、具体例として、単位球面上の運動を考えて、 数値計算により物体の軌跡を求めてみます。

曲面上のラグランジアン

正規直交座標系 x, y, z を用いると、 ラグランジアンは L = T − V =

1

2

m ( x'²+ y'²+ z'²)− φ (x, y, z) と表されます。 （ただし、変数に付いた ' は時間微分を表すものとする。） これが、曲面上ではどう表されるかを考えてみましょう。

3次元空間上の曲面は、 2つの媒介変数 q₁, q₂ を用いて、

と表すことができます。 で、物体がこの曲面上に垂直抗力で拘束されているものとします。

垂直抗力は仕事をしない（＝ ポテンシャルには影響しない）ので、 まず、φ の方は、 単純に x = x(q₁, q₂) 等を代入して、

と表すことができます。 （座標 q で表したポテンシャルであることを明示するために φ_q と書きましたが、 文脈で分かる場合には、 φ(q₁, q₂)= φ ( x(q₁, q₂), y(q₁, q₂), z(q₁, q₂)) と書いたりもします。）

一方、運動エネルギーの方は、 微分の座標変換法則が必要になります。 全微分公式から、

となるので、 運動エネルギーはこれらの二乗和になるんですが、 二乗和をいちいち書くのも面倒なので、 以下のような記号を用意します。

ちなみに、g は対称行列になっています。 この行列 g を曲面上の計量といいます。 また、g の各要素は q の関数になっているので、 これを明示するために、ここでは g(q) と書きましょう。 すると、運動エネルギー T は以下のように表すことができます。

したがって、ラグランジアンは以下のようになります。

曲面上のハミルトン形式

ラグランジアンの形が分かれば、 ラグランジュ形式やハミルトン形式の運動方程式が立てられます。 まずは、 q =(q₁, q₂)^T の共役運動量 p =(p₁, p₂)^T を求めると、以下の通り。

したがって、ハミルトニアンは、

となります。 g⁻¹(q) は g(q) の逆行列なんですが、 一応、要素を書き下しておくと、以下の通りです。

で、ハミルトン形式の運動方程式を立てると、以下の通り。

具体例: 単位球面上の運動

前節までで、曲面上のハミルトン形式を導出したわけですが、 まだ「よく分からない式」だと思うので、 具体例を挙げてみましょう。

簡単な例ということで、 高さに比例したポテンシャル φ = mGz （m は質量で、G は重力定数） が働いているときの、 単位球面上の運動を考えてみます。

ハミルトン形式の導出

単位球面は、

という式で表すことができます。 （要するに、q₁ が水平角で、 q₂ が仰角。 q₂=0 付近が一番下で、 q₂= π 付近が一番上。）

微分すると、

で、計量 g は、

となります。 また、 g⁻¹, φ =−mG cos q₂ をq で微分すると、

が得られます。 これらを前節の結果に代入すると、以下の微分方程式が得られます。

数値計算

微分方程式を導出したからといって、 厳密解を解析的に求められるわけではないんですが、 ハミルトン形式にしてしまえば、 数値計算は簡単にできます。

ハミルトン形式のように、

d

dt

x = f(x) という形になっている微分方程式は、 例えば、以下のような反復計算で近似解が得られます。

この方法はオイラー法と呼ばれているもので、 f のテイラー展開に関して1次程度の近似度の数値計算法です。 （式自体は簡単だけど、精度はあまり高くない。）

では、オイラー法を用いて、球面上の物体の運動を数値計算してみましょう。 プログラム例を C# 3.0 で示すと、以下のようになります。 （参考： 「C# によるプログラミング入門」「C# 3.0 の新機能」）

using Func2 = System.Linq.Func<double, double, double>;
using Func4 = System.Linq.Func<double, double, double, double, double>;

static void Simulate()
{
  const double M = 0.1;
  const double G = 10;

  Func2 x = (q1_, q2_) => (Math.Cos(q1_) * Math.Sin(q2_));
  Func2 y = (q1_, q2_) => (Math.Sin(q1_) * Math.Sin(q2_));
  Func2 z = (q1_, q2_) => (-Math.Cos(q2_));

  Func4 fq1 =
    (q1_, q2_, p1_, p2_) => (p1_ / (M * Math.Sin(q2_)));
  Func4 fq2 =
    (q1_, q2_, p1_, p2_) => (p2_ / M);
  Func4 fp1 =
    (q_1, q2_, p1_, p2_) => (0);
  Func4 fp2 =
    (q1_, q2_, p1_, p2_) => (
      (p1_ * p1_ * Math.Cos(q2_))
        / (M * Math.Sin(q2_) * Math.Sin(q2_) * Math.Sin(q2_))
      - M * G * Math.Sin(q2_)
      );

  double q1 = 0;
  double q2 = Math.PI / 2;
  double p1 = 0.1;
  double p2 = 0;

  const double dt = 0.01;
  const double t_end = 10;
  const int DISPLAY_INTERVAL = 5;

  Console.Write("t,x,y,z\n");

  int n = 0;
  for (double t = 0; t < t_end; t += dt, ++n)
  {
    q1 += dt * fq1(q1, q2, p1, p2);
    q2 += dt * fq2(q1, q2, p1, p2);
    p1 += dt * fp1(q1, q2, p1, p2);
    p2 += dt * fp2(q1, q2, p1, p2);

    if (n == DISPLAY_INTERVAL)
    {
      Console.Write("{0},{1},{2},{3}\n",
        t,
        x(q1, q2), y(q1, q2), z(q1, q2)));
    }
  }
}

ちなみに、 もう少し込み入ったことをしてる（精度の高い数値計算法を使ったり）サンプルプログラムも作ったので、 置いておきます →

surface.cs 。 （GUI 版、.NET Framework 3.0 必須 → ソースファイル一式。）

別ページにてバージョンアップ → その1：「曲面上の物体の運動シミュレーション」、 その2：「Expression Tree ＋ CodeDom ＋ WPF」。