Mathematica >
Mathematicaチュートリアル

NDSolveの合成・分割法

はじめに

微分方程式系を部分系に分割し,それぞれを適切な積分法を用いて解くと便利なことがある.それぞれの解を合成することで,体積等の動的特性を保存することができることがある.分割と合成に関する詳細は[MQ02, HLW02]に,NDSolveに関連する特定の問題は[SS05, SS06]に記載されている.

定義

y (0)=y0nであるときの初期値問題y' (t)=f (y (t))について考える.

合成

合成は,数値積分法の次数を上げるための便利な手法である.
外挿で使われるAitken-Nevilleアルゴリズムとは対照的に,合成ではシンプレクティック法等の基礎的な積分法の幾何的特性を保存することができる.
を,指定された実数1, ..., sを持ち,刻み幅ih を取る基本的積分法とする.
すると,s 段階の合成法f, h は以下のように与えられる.
同じ基本的メソッド= (i), i=1, ..., s.を含む合成法f, h に興味があることがある.
面白い特殊ケースとして,対称合成i=s-i+1(ここでi=1, ..., s/2)がある.
最も一般的な合成には次のようなものがある.
  • 対称2次法の対称合成
  • 1次法(随伴行列*を持つメソッド)の対称合成
  • 1次法の合成

分割

s 段階の分割法は,加法的に分割されたf の合成法の一般化である.
大切なのは,fではなくfiを含む問題を解く際に,計算上の利点がある場合が多いということである.
s 段階分割法は,以下の形式の合成である.
ここでf1, ..., fs は必ずしも異なるものでなくてもよい.
これでそれぞれの基本的積分法が,問題の一部だけを解くようになったが,適切な合成により有利な特性を持つ数値スキームが生まれる可能性がまだある.
べクトル場fi が積分可能ならば,厳密解もしくは流れfi, h を数値積分法の代りに使うことができる.
分割法では,流れと数値法を混合して使うこともできる.
その例にLie-Trotterの分割[T59]がある.
f=f1+f21=2=1で分割する.すると が1次の積分法を生成する.
計算上では,集団特性を使って流れを合成すると利点が多い.

実装

分割法と合成法を実装するために,新しいNDSolveフレームワークにいくつかの変更を加えることが必要であった.
  • メソッドにより,任意数の部分メソッドが呼び出せるようにする.
  • べクトル場全体ではなく部分体を数値的に評価するための関数を順に渡していく機能.
  • 流れを計算するLocallyExactメソッド.部分系を分析的に解いて(局所的な)解を数値的に求める.
  • 初期化の繰返しを避けるために同一メソッドのデータをキャッシュに格納.同一の部分体を数値的に評価するためのデータもキャッシュに格納される.
簡略化された入力文法により,省略されたべクトル場とメソッドを循環的に挿入することができる.これらは明確に定義しなければならない.
{f1, f2, f1, f2}{f1, f2}で入力できる.
{f1, f2, f3, f2, f1}{f1, f2, f3}で入力できない.{f1, f2, f3, f2, f1}{f1, f2, f3, f1, f2}に対応するからである.

ネストされたメソッド

以下の例では合成を使って,低次の分割法から高次の分割法を構築する.

簡略化

上記の例で,流れの集合特性を使い,Splittingメソッドを直接呼び出すと,より効率的な積分法を得ることができる.

例題

次の例題は,ストラング(Strang)の分割[S68],[M68]として知られる2次の対称分割を使う.分割係数は,方程式の構造から自動的に決定される.
メソッドSymplecticPartitionedRungeKuttaについてシンプレクティック法のリープフロッグ公式として知られるメソッドを定義する.
In[2]:=
Click for copyable input
役に立つ例題を含むパッケージをロードする.
In[3]:=
Click for copyable input

シンプレクティック分割

シンプレクティックリープフロッグ

SymplecticPartitionedRungeKuttaは長期力学を使った,分離可能なハミルトン系H (p, q)=T (p)+V (q)を解くための効率的なメソッドである.
Splittingの方が汎用のメソッドであるが,分割されたシンプレクティック法を構築するために使うことができる(しかし,SymplecticPartitionedRungeKuttaよりもやや効率性に劣る).
分離可能なハミルトン系H (p, q)=p2/2+q2/2によって支配される線形微分方程式系から派生する調和振動を例に取る.
In[5]:=
Click for copyable input
Out[5]=
ハミルトンべクトル場を,運動量および位置を支配する独立の成分に分割する.これは,方程式の右辺をゼロに設定することで実行できる.
In[6]:=
Click for copyable input
この重み付き(1次)オイラー積分ステップの合成は,シンプレクティック(2次)リープフロッグ公式に当たる.
In[11]:=
Click for copyable input
Out[14]=
ExplicitEulerメソッドは1度しか指定できなかった.2回目と3回目は循環的に埋められるからである.
積分ステップの最後の結果である.
In[15]:=
Click for copyable input
Out[15]//InputForm=
シンプレクティックリープフロッグ積分法に対応する組込みの積分法を実行する.
In[16]:=
Click for copyable input
Out[16]=
積分ステップの最後の結果は,分割法の結果と同一である.
In[17]:=
Click for copyable input
Out[17]//InputForm=

シンプレクティックリープフロッグの合成

これは,基本的な積分法としてシンプレクティック法のリープフロッグ公式を使い,Ruth-Yoshidaの対称合成[Y90]を使って4次のシンプレクティック積分法を構築する.
In[18]:=
Click for copyable input
Out[20]=
積分ステップの最後の結果である.
In[21]:=
Click for copyable input
Out[21]//InputForm=
これはRuth,Yoshidaの4次メソッドの係数を使って組込みのシンプレクティック積分法を実行する.
In[22]:=
Click for copyable input
Out[23]=
積分ステップの最後の結果は,合成法の結果と同一である.
In[24]:=
Click for copyable input
Out[24]//InputForm=

複合メソッド

閉形式の解はべクトル場全体に存在しないことがよくあるが,ベクトル場の一部について微分方程式系を解析的に解くことが可能な場合がある.
DSolveで解が見付かる場合は,局所的に解を求めるのに直接数値評価を使うことができる.
この考えはメソッドLocallyExactに実装されている.

調和振動テスト例題

以下の例題は,調和振動方程式の分割成分の厳密な流れの解が,各分割成分にオイラー法を適用したものと同じであることを調べる.
In[25]:=
Click for copyable input
In[33]:=
Click for copyable input
In[34]:=
Click for copyable input
Out[34]//InputForm=
In[37]:=
Click for copyable input
Out[38]//InputForm=

数値・記号複合分割法(ABC flow)

3次元流れを体積保存するためのモデルで,広く研究されているABC flow (Arnold,Beltrami,Childress flow)を考えてみる.
In[39]:=
Click for copyable input
Out[39]=
体積保存積分法を直接適用すると,一般に対称性は保存しない.陰的中点公式等の対称保存積分法では,体積を保存しない.
右辺の要素を零に設定することで系の分割を定義する.
In[40]:=
Click for copyable input
In[45]:=
Click for copyable input
Out[45]=
In[46]:=
Click for copyable input
Out[46]=
Y1系はDSolveで厳密に解くことができるので,LocallyExactメソッドを使うことができる.
しかし,Y2は解けないので,分割法の目標の属性を得るためには,適切な数値積分法を使う必要がある.
組込みのImplicitRungeKuttaメソッドについて陰的中点公式を計算するメソッドを定義する.
In[47]:=
Click for copyable input
2次の体積保存・反転対称群積分法を定義する[MQ02].
In[48]:=
Click for copyable input
Out[48]=

ロトカ・ボルテラ(Lotka-Volterra)方程式

この系のさまざまな数値積分法が「ロトカ・ボルテラ方程式の数値解法」の中で比較してある.

オイラー方程式

オイラー方程式のさまざまな数値積分法が,「剛体ソルバ」の中で比較してある.

非自律系のべクトル場

強制平面非自律微分方程式系のDuffing振動子を考える.
In[49]:=
Click for copyable input
Out[49]=
これは系の分割を定義する.
In[50]:=
Click for copyable input
べクトル場の時間成分の分割はあいまいであるため,メソッドはエラーメッセージを発する.
In[52]:=
Click for copyable input
Out[52]=
Z[T]Tとなるような新しいダミー変数Z[T]を導入し,分割された微分方程式系でそれをどのように分配するかを指定することにより,方程式を拡張することができる.
In[53]:=
Click for copyable input
これは,あらゆる有限時間に対して1+2=-(ここで12はリアプノフ指数)を満足する幾何的分割法を定義する[MQ02].
In[59]:=
Click for copyable input
Out[59]=
解のプロットである.
In[60]:=
Click for copyable input
Out[60]=

オプションの要約

Compositionのデフォルトの係数オプションは,Methodオプションを使って指定するメソッドの特性によって,SymmetricCompositionCoefficientsSymmetricCompositionSymmetricMethodCoefficientsのどちらかを自動的に選ぶことを試みる.
オプション名
デフォルト値
CoefficientsAutomatic合成法で使用する係数を指定する
DifferenceOrderAutomaticメソッドの局所的な精度の次数を指定する
MethodNone数値積分で使用する基礎メソッドを指定する

Compositionメソッドのオプション

オプション名
デフォルト値
Coefficients{}分割法で使用する係数を指定する
DifferenceOrderAutomaticメソッドの局所的な精度の次数を指定する
Equations{}方程式を分割する方法を指定する
MethodNone数値積分で使用する基礎メソッドを指定する

Splittingメソッドのオプション