WOLFRAM言語チュートリアル

NDSolveの"Extrapolation"メソッド

はじめに

外挿法とは,次数と刻み幅が自動制御される任意次数のメソッドのひとつである.誤差推定は,同じメソッドをさまざまなステップ数で使い,計算された解全体にフィットする多項式の外挿を使って,任意区間上の解を計算することで生じ,複合的な高次法を与える[BS64].また,この多項式は誤差推定の方法を提供する.

一般に低精度では,外挿法はルンゲ・クッタ型のメソッドに及ぶものではなかった.しかし,高精度で任意次数法を使うということは,非常に正確な許容誤差の計算では固定次数法を使うよりも格段に速くなり得るということである.

次数および刻み幅の制御は,[HNW93]と[HW96]に記述されているコードodex.fseulex.fに基づいている.

以下のようにして,ステップシーケンスのプロットのためのユーティリティ関数およびすでに定義された問題パッケージをロードする.
In[1]:=
Click for copyable input

"Extrapolation"

法は各1ステップの積分法にリチャードソン(Richardson)の補外を1度実行する.これは「NDSolveのDoubleStep法」内に記述されている.

は,リチャードソンの補外を一連の微調整に一般化するものである.

次の微分方程式系を考えてみる.

を基本的な刻み幅とし,次のように正の整数の単調増加列を選ぶ.

これに対応する刻み幅

で定義する.次に次数 の数値法を選び,刻み幅 ステップ実行することにより初期値問題の解を計算すると,以下が得られる.

値の表を構築することにより,AitkenNevilleアルゴリズムを使って外挿が実行される.

ここで は,基本メソッドが補外のもとで対称かどうかにより1か2のいずれかになる.

(1)の値の依存性グラフにより,相互関係が明らかになる.

と考えることは,リチャードソンの補外に等しい.

硬くない問題では,(2)の の次数は である.硬い問題では,分析はもっと複雑であり,特異な振動問題に現れる振動項の調査が必要になる[HNW93, HW96].

外挿の数列

外挿の数列は,実装で指定することができる.以下は,一般的に使われる数列の例である.

Romberg数列.
In[3]:=
Click for copyable input
Out[3]=
Bulirsch法.
In[4]:=
Click for copyable input
Out[4]=
調和数列.
In[5]:=
Click for copyable input
Out[5]=

を満足する数列には, 次の基本積分法の丸め誤差を最小化する効果がある.

二次の基本メソッドでは,数列の最初の項目は以下のように与えられる.
In[6]:=
Click for copyable input
Out[6]=
以下は,メソッドの次数がオプションである調和数列を定義する関数を加える例である.
In[7]:=
Click for copyable input

最低コストの数列は調和数列である.しかし丸め誤差が捨てられないので,問題がないわけではない.

丸め誤差の累積

高次の外挿では,AitkenNevilleアルゴリズム(3)における丸め誤差の累積を考慮することが大切である.

[HNW93]のセクションII.9のExercise 5を例に取る.

項目 が丸め誤差 を持ち,これらの誤差の外挿表の中への伝播を計算するとする.

外挿過程(4)は線形なので, はゼロに等しく,を取ると仮定する.

以下は,調和数列と対称二次基本積分法 を使用して,初期データ上のAitkenNevilleアルゴリズム(5)の進化を示している.

16次のメソッドでは,丸め誤差の累積のためにおよそ2桁が失われてしまう.

このモデルは幾分大雑把なものである.というのも,「丸め誤差の軽減」で説明するように,丸め誤差は のとき ではなく で生じることの方が多いからである.

丸め誤差の軽減

ここで,高次の外挿における丸め誤差による効果が軽減できるアプローチを探すのも意味のあることである.

ひとつのアプローチとして,丸め誤差を小さくするために異なるステップシーケンスを選ぶというものがある.しかしこれには,外挿表の最初の列に を形成するのに必要な積分ステップ数が多くなるという欠点がある.

STEP等,コードの中には,厳しい許容誤差に対する丸め誤差の効果を軽減するために積極的な方法を取るものもある [SG75].

別の方法として,外挿ということにおいてはそれほどの注目を集めてはいないが,浮動小数点操作における丸め誤差の大きさを減少させるために,基本積分法を修正するというものがある.このアプローチは[G51]の考えに基づいており,2体問題[F96b]で使われている.(背景は[K65],[M65a],[M65b],[V79]も参照のこと)次のセクションでは,これについて説明する.

基本メソッド

以下のメソッドは,外挿の基本積分法として一般的に選ばれるものである.

  • (平滑化ステップを含むグラッグ(Gragg)法(6))
  • (平滑化ステップを含まないBaderDeuflhard公式(7))
    • (平滑化ステップを含むBaderDeuflhard公式(8))

    効率性を考えて,これらのメソッドはNDSolveに組み込まれており,Methodオプションで別々のメソッドとして呼び出すことができる.

    上記のメソッドの実装には,および内の複数のサブステップに対する特殊な解釈が含まれる.

    1ステップメソッドに対するNDSolveフレームワークでは,増分あるいは解のアップデートを返す公式が使われる.これは,数値誤差が長期に渡る積分で取り除かれない幾何的数値積分に有利である.またこれは,補償される総和などの効率的補修ストラテジーの適用を許可する.この式も補外のコンテキストで便利である.

    これらのメソッドは,累積される丸め誤差を軽減するために使われる,増分の再定式化と共に記述されている.

多段階のオイラー(Euler)法

があるとき,刻み幅 回の連続した積分ステップがオイラー法で実行されるとする.

ここで である.

陽的ルンゲ・クッタ法との一致

ある基本積分スキームでは,(9)から生成された外挿表の項目 は,陽的ルンゲ・クッタ法と一致するということがよく知られている([HNW93]のセクションII.9のExercise 1を参照).

例えば,(10)は 段の陽的ルンゲ・クッタと等しい.

係数は以下の表の通りである.

再定式化

とすると,積分(11)は陽的ルンゲ・クッタ法(12, 13)との一致を反映するために以下のように書き換えることができる.

ここで(14)の右辺の項は,同値の から逸脱するものと考えられる.

(15)の は(16)の に相当する.

とすると,必要な結果は以下のように回復できる.

数学的に言うと,公式(17)と(18, 19)は等しい.しかし のとき,(20)の計算の方が有限精度の浮動小数点演算の丸め誤差累積を軽減させる の数量(増分)の累積和が小さいという利点がある.

多段階法の陽的中点公式

リチャードソンの補外を効果的に実装するためには, の偶数乗での拡大が非常に重要である[S70].

一段階のオイラー法の後に多段階法の陽的中点公式を使って実行される刻み幅 の連続した積分ステップ を考える.

(21)を 段階の中点公式で計算すると,メソッドには誤差の対称拡大が含まれる([G65], [S70]).

再定式化

増分に関する(22)の再定式化は,以下のように求められる.

グラッグ法(平滑法)

グラッグ法は,陽的中点公式が安定ではないということに歴史的起源がある.

(23)を利用するためには,公式(24)を 段階の中点公式で計算する.これには安定領域を増大し,基本ステップの最後で関数を評価するという利点がある[HNW93].

構造上の理由で,増分の和と2つの連続した増分がアルゴリズムの最後で利用可能になる点に注目のこと.これにより,以下の公式が導かれる.

さらに(25)は,有限精度の計算では(26)よりも有利である.それは,通常増分 よりも大きい値 が計算に関与しないためである.

グラッグ法は,メソッドの後に外挿を使う場合はそれほど重要ではない.これよりもShampineの提案する平滑法の方がわずかに効率的である[SB83].

メソッドは段階法を使い,メソッドは 段階法使った後で平滑法(27)を使う.

安定領域

以下は,平滑法が線形安定領域に及ぼす影響を示したものである(FunctionApproximations.mパッケージを使って実行).

明示的な中点法(左)と平滑化を使った明示的な中点法(右)での における線形安定領域である.
Click for copyable input

任意の基本メソッドに対する正確な安定境界は計算が複雑なので,簡単な近似が用いられる.次数 の補間メソッドでは,負の実軸の交点がその点と考えられる.ここでは以下が成り立つ:

安定領域は負の半平面の場合,この半径で原点が(0,0)の円板として近似される.

陰的微分方程式

微分方程式系の一般化(28)は放物型偏微分方程式の空間離散化等多くの状況で発生する.

ここで は広く「質量行列」を言われている定数行列である.

線形方程式系の解を含む補外の基本メソッドは,形式(29)の問題を解くよう簡単に変更することができる.

線形多段階の後退オイラー法

多くの半陰的および陰的手法の記述では,増分は自然に発生する. および である系(30)に対する線形後退オイラー法を使って実行される連続した積分ステップを考えてみる.

ここで は質量行列を, のヤコビアンを表す.

(31)の増分に対する方程式の解は,行列 のLU分解を1度,それに続き各右辺の三角線形系の解を使うことにより得られる.

所望の結果は(32)から以下のように得られる.

再定式化

からの出発した増分の再定式化は,以下のように求めることができる.

(33)を使った の結果は(34)から得られる.

のとき,(35)と(36)は等価であることに注目のこと.

線形多段階の陰的中点公式

BaderとDeuflhardの公式[BD83]を使って,線形陰的オイラー法に続く および である線形多段階の陰的中点公式の1ステップを考える

線形段階の陰的中点法で(37)を計算すると,メソッドには誤差の対称拡大が含まれる[BD83].

再定式化

増分に関する(38)の再定式化は,以下のように求めることができる.

平滑法

線形の陰的中点公式に適した平滑法は[BD83]である.

増分に関して書き直されたBaderの平滑法(39)は,以下のようになる.

必要量は,(40)が ステップで実行されたときに得られる.

線形の陰的中点公式の平滑法は,陽的中点公式のグラッグの平滑法とは異なる役割を持つ([BD83]と[SB83])を参照).前者には排除の必要がある安定ではない項がないので,その目的は漸近安定性を向上させることである.

メソッドは段階法を使い,メソッドは 段階法の後で平滑法(41)を使う.

増分による多項式の外挿

増分によって外挿表の第1列の項目である を修正する方法は,前述の通りである.

しかし,ある基本的積分法では, のそれぞれが陽的ルンゲ・クッタ法に対応する.

従って,その関係がまだ完全に利用されておらず,さらなる微調整が可能なようである.

AitkenNevilleアルゴリズム(42)には線形差分方程式が含まれるので,増分を使って外挿プロセス全体が実行できる.

これにより,以下のようなAitken-Nevilleアルゴリズムの修正が導かれる.

(43)の数量は,修正された基本積分法から得られた初期量 から始めて,の影響を受けることなく繰り返し計算できる.

最終目標値 と回復できる.

この利点は,外挿表が小さい数量で構築されているため,減算の桁落ちによる丸め誤差の効果が,減少する.

実装の問題

実装の問題として考慮すべきものが多くあるが,そのいくつかをここに挙げる.

ヤコビアンの再利用

ヤコビアンは,それ自身とそれが評価される関連時間とを共に保存することにより,各時間刻みにおけるすべての項目 について1度だけ評価される.これには,拒絶されたステップについてはヤコビアンを再計算する必要がないという利点がある.

密な線形代数

密な系では,LU分解にはLAPACKルーチンのxyyTRFが,その結果の三角方程式系にはxyyTRSが使える[LAPACK99].

次数の適応法と作業推定

積分の過程で外挿の次数を適宜変更するためには,基本スキームと外挿の列に必要な作業量を測定することが大切である.

関数評価の相対コストを測定すると有利である.

線形の陰的スキームで各線形系を解くコストを考慮に入れるためには,系の次元(できれば構造に従った重み付きで)を組み入れる必要がある.

安定性のテスト

外挿法は,問題を引き起す可能性のある大きな基本刻み幅を使う.

「どちらのコードもオーバーフローを起すので,ファン・デル・ポル(Van der Pol)方程式を解くことはできない」[S87].

線形の陰的基本スキームには2つの形式の安定性テストが使われる(追加説明は[HW96]を参照のこと).

そのうちのひとつは,外挿の過程で実行される.とする.

Deuflhardは,の計算を中断するためには,完全に陰的なスキームに適用されたニュートン反復が収束するかどうかを判定すればよいと提案している.

線形後退オイラー法では,この提案により以下が導かれる.

(44) は(45)と2つ目の方程式だけが異なる.異なる右辺の解を見付けなければならないが,関数評価を追加する必要はない.

陰的中点公式では,単に基本演算操作を数回必要とするだけである.

どちらの安定性テストの形式の実装でも,増分公式の方が正確である.

例題

作業誤差比較

異なる外挿スキームを比較するために,[HW96]の例題を見てみる.
In[9]:=
Click for copyable input

厳密解は で与えられる.

増分公式

この例には,調和数列を持つ八次のが含まれている.外挿の過程で増分ベースの公式を使うことにより,精度がおよそ2桁向上している.

  • 標準公式(46)の結果は,緑で描画されている.
  • 増分公式(47)の後に標準外挿(48)を施した結果は,青で描画されている.
    • 増分公式(49)に,(50)で実行された増分の外挿を施した結果は,赤で描画されている.

    外挿の過程で増分ベースの公式を使うことにより,精度がおよそ2桁向上している.

    以下は前出の例題の外挿表の最初の列を形成する積分データの相対誤差を比較したものである.

    参照値は32桁のソフトウェア演算を使って計算され,最近傍のIEEE倍精度浮動小数点数に変換された.ここでULPとはUnit in the Last PlaceあるいはUnit in the Last Position(最下位桁単位)を意味する.

    実際には,の誤差が1ULPを超える可能性があるので,丸め誤差の増大の研究の動機として使われた丸め誤差モデルには限界があるという点に注目する.

    外挿の過程で使われた増分公式では,には1ULPよりも小さい丸め誤差が生じる.

メソッドの比較

ここでは(赤),(青),(緑)のそれぞれに基づく外挿に必要な作業を比較する.

計算はすべて32桁のソフトウェア演算を使って実行される.

次数の選択

解く問題を選ぶ.
In[15]:=
Click for copyable input
評価の次数と回数を保存する監視関数を定義する.
In[16]:=
Click for copyable input
積分が進むに従って,監視関数を使ってデータを収集する.
In[17]:=
Click for copyable input
積分中にどのように次数が変化するかを表示する.
In[18]:=
Click for copyable input
Out[18]=

メソッドの比較

解く問題を選ぶ.
In[19]:=
Click for copyable input
問題が硬いかどうかにより,1対のメソッドの間で切り替わるメソッドを使って参照解を計算する.
In[20]:=
Click for copyable input
比較するメソッドのリストを定義する.
In[22]:=
Click for copyable input
角度と作業を比較するデータは,を使って許容誤差範囲について計算される.
In[23]:=
Click for copyable input
このメソッドについての作業と誤差の比較データは以下の対数プロットで表示される.ここで大域的誤差は縦軸で,関数評価の回数は横軸で表されている.最終的に補外メソッドの次数が高いほど効率的であるということになる.また,非常に厳しい許容誤差でさえも,増分公式はよい結果を出し続けている点に注目する.
In[24]:=
Click for copyable input
Out[24]=

硬い系

回路を記述する最も簡単な非線形方程式のひとつに,ファン・デル・ポル方程式がある.
In[25]:=
Click for copyable input
基本メソッドでデフォルトの二重調和列2, 4, 6, にしたで方程式を解く.硬さ検出デバイスが積分を終了させ,別のメソッドを提案する.
In[28]:=
Click for copyable input
Out[28]=
以下ではベースメソッドでデフォルトの部分調和列2, 3, 4, にしたを使って方程式を解く.
In[29]:=
Click for copyable input
Out[29]=

ヤコビ行列が自動的に計算され(数値差分あるいは記号微分を使ってユーザが指定できる),ランタイム時に適切な線形代数ルーチンが選ばれ実行される.

一定時間上の最初の解の成分をプロットする.
In[30]:=
Click for copyable input
Out[30]=
解を計算するときに使われた刻み幅のプロットである.
In[31]:=
Click for copyable input
Out[31]=

高精度の比較

ローレンツ(Lorenz)方程式を選ぶ.
In[32]:=
Click for copyable input
次数九(8)の陽的ルンゲ・クッタ法に埋め込まれたbigfloat (ソフトウェアの浮動小数点数)を使う[V78].
In[33]:=
Click for copyable input
Out[33]=
LSODAのbigfloatバージョンを使って法を実行する.これらのメソッドの最大次数は12である.
In[34]:=
Click for copyable input
Out[34]=
を基本積分スキームとしたメソッドを実行する.
In[35]:=
Click for copyable input
Out[35]=
以下はさまざまなメソッドで使われる刻み幅である.外挿で使われる高次数は大きい刻み幅が取れることを示している.
In[36]:=
Click for copyable input
Out[38]=

質量行列 - fem2ex

以下の偏微分方程式を考える:

整数 を定義し,Galerkinの離散化法を使って において近似する:

ここで において においての区分線形関数である.

(51)に適用された離散化(52)は(53)のような定数質量行列式を持つ常微分方程式系になる.常微分方程式系は[SR97]におけるfem2ex問題であり,IMSLライブラリでも見られる.

この問題は,行列の疎配列を使うよう設定される.その疎配列は小さい次数を考慮する必要はないが,離散化点の数が増えるとうまくスケールされる.初期条件にはベクトル値の変数が使われる.この系は区間上で解かれる.
In[39]:=
Click for copyable input
でベースメソッドのを使って常微分方程式系を解く.系が質量行列形式であることを指定するためにオプションが使われている.
In[48]:=
Click for copyable input
このプロットは,メソッドで取られる比較的大きい刻み幅を示している.
In[49]:=
Click for copyable input
Out[49]=
この種の問題のデフォルトメソッドは,汎用の微分代数方程式ソルバ[HT99]のである.このメソッドは非常に汎用であるため,この例題には幾分過ぎている.
In[50]:=
Click for copyable input
以下のプロットで,微分代数ソルバではずっと多くのステップが取られることがはっきり分かる.
In[51]:=
Click for copyable input
Out[51]=
格子上に解をプロットするために使用する関数を定義する.
In[52]:=
Click for copyable input
格子上に解を表示する.
In[53]:=
Click for copyable input
Out[53]=

微調整

"StepSizeSafetyFactors"

ほとんどのメソッドでそうであるように,あまりに小さいステップと,通常は拒否されるほどの大きいステップとの間でバランスを取りたいと思うであろう.オプションは以下のようにして,刻み幅の選択を制約する.次数 のメソッドによって選ばれた刻み幅は,以下を満足する.

これには次数依存要因と次数非依存要因の両方が含まれる.

"StepSizeRatioBounds"

オプションは以下のように,次の刻み幅の範囲を指定する.

"OrderSafetyFactors"

の重要な一面に,次数の選択というものがある.

次数 の各外挿には,関連する作業推定 がある.

陽的な基本メソッドの作業推定は,関数評価の回数と使われるステップシーケンスに基づいている.

線形陰的な基本メソッドの作業推定は,ヤコビアンの評価のコスト,LU分解のコスト,線形方程式の逆計算のコストも含んでいる.

単位ステップあたりの作業推定は,作業推定 と次数 のメソッドで取る新しく使う推定刻み幅 ((54)から計算される)から形成される.

は連続した推定値を比較して,別の次数のメソッドの方が効率的となる時期を決定する.

オプションは,推定値 の比較に含まれる安全率を指定する.

次数は のとき減少する.

のときは増加する.

この他にも,いつステップあたりの次数の最大増加が1(対称法では2)となるか,ステップが拒絶された直後に,いつ次数の増加を禁止するか等の制約がある.

ノンスティッフな基本メソッドでは,デフォルト値はであり,スティッフな基本メソッドではである.

オプションの要約

オプション名
デフォルト値
"ExtrapolationSequence"Automatic外挿で使う数列を指定する
"MaxDifferenceOrder"Automatic使用する最大次数を指定する
Method"ExplicitModifiedMidpoint"使用する基本積分法を指定する
"MinDifferenceOrder"Automatic使用する最小次数を指定する
"OrderSafetyFactors"Automatic適応型の次数選択の推定で使う安全率を指定する
"StartingDifferenceOrder"Automatic使用する初期次数を指定する
"StepSizeRatioBounds"Automatic下限上限のときに,現在の刻み幅 から新しい刻み幅 への相対変化の境界を指定する
"StepSizeSafetyFactors"Automatic適応型刻み幅に使う誤差推定に組み込む安全率を指定する
"StiffnessTest"Automatic硬さ検出機能を使うかどうかを指定する

メソッドのオプション

オプションのデフォルト設定Automaticでは,基本メソッドの硬さと対称性に基づいて列が選ばれる.

オプションのデフォルト設定Automaticは,最大次数を十進の作業精度の2倍に制限する.

オプションのデフォルト設定Automaticは,基本メソッドの次数から始めた2つの外挿の最小を選ぶ.

オプションのデフォルト設定Automaticは,硬い基本メソッドには値を,硬くない基本メソッドには値を使う.

オプションのデフォルト値Automaticは, の設定により異なる.これは,基本メソッドが対称であるかどうかにより,に設定される.

オプションのデフォルト設定Automaticでは,硬い基本メソッドには値が,硬くない基本メソッドには値が使われる.

オプションのデフォルト設定Automaticは,硬い基本メソッドには値を,硬くない基本メソッドには値を使う.

オプションの基本設定Automaticは,硬くない基本メソッドが使われるならば硬さ検定が有効になることを意味する.

オプション名
デフォルト値
"StabilityTest"True連続した陰的な解(例(55)を参照)に安定性テストを実行するかどうかを指定する

の各メソッドのオプション