MATHEMATICAチュートリアル

剛体ソルバ

はじめに

重心が原点となっている自由剛体の運動方程式は,以下のオイラー方程式で与えられる([MR99]を参照のこと).

以下はこの系の第1積分の2つの二次方程式である.

最初の制約条件は,実質的に運動をから球面へ制限する.第2の制約条件は系の運動エネルギーを表し,最初の不変式と連動して,運動を球面の楕円に制限する.

さまざまな方法での数値的な実験は[HLW02]に挙げてある.ここでは,NDSolveのさまざまな用法を比較する.

多様体生成および効用関数

まず,便利なパッケージをいくつかロードする.
In[6]:=
Click for copyable input
剛体運動と系の不変式に対してオイラー方程式を定義する.
In[8]:=
Click for copyable input
運動方程式は単位球上で閉曲線として進化する.これにより,その単位球を表す3Dグラフィックスオブジェクトが生成される.
In[13]:=
Click for copyable input
以下の関数は,与えられた多様体の上にNDSolveの解を付加する.
In[14]:=
Click for copyable input
この関数は,さまざまな解の成分をプロットする.
In[15]:=
Click for copyable input

方法の比較

オイラー方程式は,いろいろな積分法を使って解くことができる.そのどれもが異なるコスト,動的特性を持つ.

アダムスの多段階法

ここでは運動方程式を解くのにアダムス法を使う.
In[21]:=
Click for copyable input
以下は,解の軌道を単位球に重ね合せて表示している.
In[22]:=
Click for copyable input
Out[22]=
視覚的には,解は球面に閉曲線を描いているように見える.しかし,誤差のプロットによりどちらの制約条件もあまり保持されていないことが分かる.
In[23]:=
Click for copyable input
Out[23]=

オイラー法と陰的中点法

まず,指定された固定刻み幅を使ったオイラー法で運動方程式を解く.
In[16]:=
Click for copyable input
これは,指定された固定刻み幅を使った陰的中点法で運動方程式を解く.
In[17]:=
Click for copyable input
以下により,オイラー法を使った運動方程式の数値解を単位球に重ね合わせたもの(左)と,陰的中点法を使った場合のもの(右)を表示する.
In[19]:=
Click for copyable input
Out[21]=
これはオイラー法を使った場合の数値解の要素(左)と陰的中点法を使った場合の要素(右)を表示する.
In[30]:=
Click for copyable input
Out[32]=

"OrthogonalProjection"法

ここでは法を使って方程式を解いてみる.
In[33]:=
Click for copyable input
直交の制約しか守られていないので,曲線は閉じていない.
In[34]:=
Click for copyable input
Out[34]=
時間に対する不変量の誤差をプロットすると,正投影法では2つの不変量のうちの1つしか保存されていないことが分かる.
In[35]:=
Click for copyable input
Out[35]=

"Projection"法

法は制約条件の集団を取り,各積分段階の終りにその解を多様体に投影する.

一般的に問題の不変量すべてを投影に使わなければならない.そうでなければ,数値解は投影されない解よりも実際には質的に劣る可能性がある.

以下は積分法を指定し,NDSolveが呼び出されるまで制約条件の決定を保留する.
In[36]:=
Click for copyable input

1つの制約条件の投影

以下は最初の制約条件を多様体に投影する.
In[37]:=
Click for copyable input
Out[39]=
最初の不変量だけが保存される.
In[40]:=
Click for copyable input
Out[40]=
第2の制約条件を多様体上に投影する.
In[41]:=
Click for copyable input
Out[43]=
2番目の不変量だけが保存される.
In[44]:=
Click for copyable input
Out[44]=

複数の制約条件の投影

以下で両方の制約条件を多様体上に投影する.
In[45]:=
Click for copyable input
Out[47]=
両方の不変量が保存される.
In[48]:=
Click for copyable input
Out[48]=

"Splitting"法

効率的な陽的積分法を作成する分割法は,McLachlan [M93]およびReich [R93]によって別々に導き出された.

常微分方程式 の流れを と書く.

微分方程式系は の3つの要素に分割される.これらはすべてハミルトニアンで,正確に解くことができる.

ハミルトン系が解かれ,各積分段階で以下のように再結合される.

ハミルトニアンのべクトル場への分割を定義する.
In[49]:=
Click for copyable input
これはオイラーの微分方程式系である.
In[58]:=
Click for copyable input
Out[58]=
3つに分割されたべクトル場である.
In[59]:=
Click for copyable input
Out[59]=
In[60]:=
Click for copyable input
Out[60]=
In[61]:=
Click for copyable input
Out[61]=

これは対称二次分割法を定義する.係数は方程式の構造から自動的に決定され,ストラング(Strang)の分割法の拡張となっている.
In[62]:=
Click for copyable input
系を解き,解をグラフィカルに表示する.
In[63]:=
Click for copyable input
Out[64]=
不変量の1つは丸めまで保存されているが,2つ目の不変量の誤差は一定である.
In[65]:=
Click for copyable input
Out[65]=
New to Mathematica? Find your learning path »
Have a question? Ask support »