目次

概要

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

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

曲面上のラグランジアン

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

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

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

( x(q1, q2), y(q1, q2), z(q1, q2))

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

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

φq(q1, q2)= φ ( x(q1, q2), y(q1, q2), z(q1, q2))

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

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

x' =
∂x
q1
q1' +
∂x
q2
q2'
y' =
∂y
q1
q1' +
∂y
q2
q2'
z' =
∂z
q1
q1' +
∂z
q2
q2'

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

q =(q1, q2)T
g =[
g11g12
g21g22
]
gij=
∂x
qi
∂x
qj
+
∂y
qi
∂y
qj
+
∂z
qi
∂z
qj

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

T =
1
2
m q'T g(q) q'

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

L =
1
2
m q'T g(q) q' φq(q)

曲面上のハミルトン形式

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

p =
∂L
q'
= m g(q) q'

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

H = pT q' L =
1
2
m q'T g(q) q + φq(q)=
1
2m
pT g−1(q) p + φq(q)

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

g−1=[
g11g12
g21g22
]
d = g11g22 g122
g11= g22/ d
g22= g11/ d
g12= g21=g12/ d

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

d
dt
q =
∂p
H =
1
2m
g−1(q) p
d
dt
p =
∂q
H =
1
2m
pT(
∂q
g−1(q)) p
∂q
φ(q)

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

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

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

ハミルトン形式の導出

単位球面は、

x(q1, q2)=cos q1sin q2
y(q1, q2)=sin q1sin q2
z(q1, q2)=cos q2

という式で表すことができます。 (要するに、q1 が水平角で、 q2 が仰角。 q2=0 付近が一番下で、 q2= π 付近が一番上。)

微分すると、

[
x'
y'
z'
] = [
sin q1sin q2 cos q1cos q2
cos q1sin q2 sin q1cos q2
0 sin q2
] [
q1'
q2'
]

で、計量 g は、

g =[
sin2 q20
01
] , g−1=[
1/sin2 q20
01
]

となります。 また、 g−1, φ =mG cos q2q で微分すると、

q1
g−1= 0 ,
q2
g−1=[
2cos q2/sin3 q20
00
]
q1
φ = 0 ,
q2
φ = mG sin q2

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

d
dt
q1=
p1
m sin2 q2
d
dt
q1=
p2
m
d
dt
p1=0
d
dt
p2=
p12cos q2
m sin3 q2
mG sin q2

数値計算

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

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

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

x(t + Δt)= x(t)+ f(x(t)) Δt

この方法はオイラー法と呼ばれているもので、 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」。

更新履歴

ブログ