硬さ検出

概要"NonstiffTest"メソッドオプション
はじめに例題
線形安定オプションの要約
"StiffnessTest"メソッドオプション

概要

多くの微分方程式には，刻み幅を制限し，それにより陽的な解法の効率をも制限するある種の硬さというものが存在する．

この問題を回避するために，多数の陰的メソッドが長年に渡り開発されてきた．

同じ刻み幅の場合，陰的メソッドは内的線形代数に関連するオーバーヘッドのために，陽的メソッドよりも効率がはるかに劣ることがある．

このコストは，ある特定の領域において陰的メソッドはかなり大きい刻み幅を取ることができるということにより埋め合せられる．

ランタイム時に硬さを自動的に検出し，必要に応じて適切なメソッドに切り替えるユーザフレンドリーなコードを提供するため，数々の試みが行われてきた．

ここでは，コードに硬さ検出を自動的に装備するために提案されてきた多数のストラテジーを紹介する．

硬さ検出をどのようにNDSolveに実装するかを記述するために，行列の優勢固有値の推定問題に特に注目してみる．

ストラテジーの効率は数値例題で例示する．

初期化

定義済みの例題とユーティリティ関数を含むパッケージをロードする：

はじめに

初期値問題の数値解を考える：

硬さとは問題，解法，初期条件，局所的誤差許容度の組合せである．

硬さはメソッドが取ることのできる刻み幅を制約するため，陽的解法の効率を制約する．

硬さは線の方法による偏微分方程式の数値解法だけでなく，多数の実践的な系で発生する．

例題

ファン・デル・ポル(van der Pol)振動子は，非線形減衰を受けた非保存系の振動子であり，常微分方程式の硬い系の例である：

このとき ;

初期条件

を考慮し，区間 上で解く．

メソッド"StiffnessSwitching"ではデフォルトで1ペアの外挿メソッドが使われる．

解

以下のようにしてパッケージから問題をロードする：
硬くないメソッドを使って系を数値的に解く：
硬さが発生したら切り替わるメソッドを使って系を解く：
次は2つの解の成分のプロットである．尖った頂点（青）は大きさほぼ450くらいまで広がっているので切り取られている：

硬さは急速な過渡解に従う領域で起ることが多い．

以下は時間に対する刻み幅をプロットしたものである：

問題は，解が急速に変化するとき，局所的確度が優勢な問題であるため硬いソルバを使う点はほとんどないということである．

効率を考えると，局所的確度（安定性ではない）が重要である領域を，メソッドが自動的に検出できたら便利であろう．

線形安定

線形安定理論は，ダールキスト(Dahlquist）のスカラー線形テスト方程式の研究から派生する：

これは初期問題(1)を研究するための簡約化されたモデルである．

安定とは下の

（ここで であり は有理安定関数）を得るために(2)に適用されたメソッドを分析することである．

絶対安定の境界は，領域を考慮することで得られる：

前進オイラー法

前進オイラー法：

を(3)に適用すると，以下が求まる：

陰影の付いた領域は不安定を表し，ここで >1である．

線形安定境界(LSB)は負の実軸との交点と取られることが多い．

前進オイラー法ではである

固有値 では，線形安定の必要条件は，刻み幅が非常にゆるい制約の を満足する必要があるということを意味する．

しかし固有値 のとき，線形安定の必要条件は，刻み幅が非常に厳しい制約である を満足する必要があるということを意味する．

例題

以下の例題は，硬い系を解くのに陽的ルンゲ・クッタ法を使ったときの刻み幅シーケンスに対する硬さの効果を示している．

以下の系は化学反応をモデル化する：
この系は，組込みの硬さ検出を無効にすることで解かれる：
安定領域に達すると，刻み幅シーケンスは振動し始める：

後退オイラー法

後退オイラー法

を(4)に適用すると

が得られる．このメソッドは左半分の平面には無条件で安定である．

つまり，安定を保つために刻み幅には制約が必要ないということである．

欠点は，方程式の陰的な系をそれぞれの積分ステップで解かなければならないということである．

タイプの不感性

タイプ不感性ソルバは各ステップで効率よく硬さを認識し反応するため，問題のタイプ（変化）に敏感でない．

このクラスで最もよく構築されたソルバの一つがLSODA [H83]，[P83]である．

CVODE等，後の世代のLSODAはもはや固さ検出デバイスを装備しない．それは，後で説明するが，LSODAは優勢固有値を推測するためにノルム限界を使うが，その限界が非常に不正確であり得るためである．

A(α)安定のBDFメソッドの低次数は，高確度の系あるいは優勢固有値に大きい虚部のある系をとくのにLSODAおよびCVODEがあまり適していないことを意味する．線形の陰的スキームの外挿に基づくメソッド等，別のメソッドではそのようなことは問題にならない．

硬さ検出については1980年代，1990年代にスタンドアロンのFORTRANコードを使って多数の研究がされた．

それ以降新しい線形代数技術や効率的なソフトウェアが利用できるようになり，それがWolfram言語ですぐにアクセスできるようになっている．

硬さは一時的な現象であり得るため，硬さのないことを検出することも同様に重要である[S77]，[B90]．

"StiffnessTest"メソッドオプション

硬さのないソルバから硬い祖父場に切り替えるために使用できるアプローチがいくつかある．

直接の推定

硬さを検出する便利な方法に，問題のヤコビアン の優勢固有値を直接推定するというものがある（[S77]，[P83]，[S83]，[S84a]，[S84c]，[R87]，[HW96]を参照）．

そのような推定は数値積分の二次的結果として利用できることが多いので，比較的安価である．

がヤコビアンの優勢固有値に対応する固有ベクトルの近似を表し，が十分小さい場合，平均値定理により，主固有値のよい近似は以下のようになる：

リチャードソン(Richardson)の補外法は，ある種の陽的ルンゲ・クッタ(Runge–Kutta)法で行うように，この公式の量を算出する一連の微調整を提供する．

コストはせいぜい2回の関数評価ですむが，少なくともそのうちの1回は数値積分の二次的結果として利用できるので比較的安価である．

が線形安定領域を表すとする—線形安定領域と負の実軸との交点である．

積 は硬さを検出するためにメソッドの線形安定領域と比較することのできる推定値を与える：

ここで は安全因子である．

解説

メソッド"DoubleStep"，"Extrapolation"，"ExplicitRungeKutta"にはオプション"StiffnessTest"があり，これは与えられた問題に，指定されたAccuracyGoalおよびPrecisionGoalで適用されたメソッドが硬いかどうかを見分けるために使われる．

メソッドオプションの"StiffnessTest"自体は(5)の弱形式を実装する多数のオプションを受け入れる．ここでテストは決まった回数だけ失敗することが許される．

その理由は，問題の中にはある領域だけがほどほどに硬いだけということがあり，明示的な積分法はまだ有効であるためである．

"NonstiffTest"メソッドオプション

"StiffnessSwitching"メソッドにはオプション"NonstiffTest"がある．これは硬いメソッドから硬くないメソッドへと切り替えるために使われる．

次の設定がオプション"NonstiffTest" に許される．

硬くないソルバへの切換え

硬いメソッドと独立のアプローチが使われる．

ヤコビアンJ（あるいは近似）があるとき，以下のうちの一つを計算する：

ノルム限界：

スペクトル半径： ρ(J)max

優勢固有値

線形代数の手法の多くは，1つの問題を高確度までとくことに焦点を当てている．

硬さ検出では，1桁か2桁までの解を持つ問題の連続が妥当である．

数値的離散化

で，部分区間のいくつかにある行列のシーケンス を考えてみる：

連続行列のスペクトルはステップからステップでゆっくりと変化することが多い．

目的は下のような行列の連続 の優勢固有値（の限界）を推定する方法を探すことである：

NormBound

優勢固有値の限界を得るための簡単で効率的な方法は，ヤコビアン のノルム（ここで通常 あるいは ）を使うことである．

この方法には，硬いソルバで実施される作業よりも小さい複雑性 がある．

これはLSODAで使われるアプローチである．

例題

以下のヤコビ行列は，硬いソルバを使ったファン・デル・ポルの数値解で現れる：
ノルムに基づく限界は，スペクトル半径を大きさの次数より多く見積もる：

固有値の直接計算

小さい問題（例えば ）では，優勢固有値を直接計算した方が効率がよいことがある．

ベキ乗法

Shampineはヤコビアンの優勢固有値を推定するためにベキ乗法を使うことを提唱した[S91]．

ベキ乗法はあまり高く評価されている方法ではないが，Google社のページランキングに使用されているため，非常に関心が高まっている．

ベキ乗法は次の場合に使うことができる．

解説

出発ベクトルを与えると は以下を計算する．

優勢固有値の近似を計算するために，レイリー商が使われる：

実際には，近似固有ベクトルは各ステップでスケールされる：

プロパティ

ベキ乗法は以下の割合で線形に収束する：

これは遅い．

特に，このメソッドは優勢固有値の複素共役ペアを持つ行列に適用されると収束しない．

一般化

ベキ乗法は，同一モジュール固有値の問題を克服するために適用することができる（例：NAPACK）．

しかし，個の修正は一般にクラスタ化された固有値の遅い収束率という問題は回避しない．

ベキ乗法を一般化するための主なアプローチが2つある：

部分空間反復

部分空間（同時）反復は，各ステップで m 個のベクトルに作用することによりベキ乗法の考えを一般化する．

ベクトル の正規直交集合から始める．ここで通常 である．

行列 との積から：

すべてのベクトルが の同じ優勢固有ベクトル の倍数へと収束するのを妨げるため，正規直交化される：

正規直交化のステップは，行列の積と比較すると高価である．

レイリー・リッツ(Rayleigh-Ritz)の投影法

入力：行列 A とベクトルの正規直交集合 V

収束

部分空間（同時）反復は，各ステップで 個のベクトルに作用することによりベキ乗法の考えを一般化する．

SRRITは割合

で線形に収束する．特に，優勢固有値の割合は

となる．したがって，優勢固有値にのみ関心があるとしても例えば 以上を取るとよいことがある．

誤差制御

優勢固有値の連続した近似に対する相対的な誤差検定は

である．これは収束が遅いときは満足されないので十分ではない．

あるいはならば， の 番目の列は一意に決められない．

SRRITで使われる残差検定は

であり，ここで ，は

の 番目の列，は の 番目の列である．

これは同一モジュラの固有値にも使えるので便利である．

次数付きのシューア分解を使用するため，上三角行列 の最初の列の位置が検定される．

実装

部分空間反復の実装にはいくつかある．

KrylovIteration

列が，指定された部分空間 の直行基底をなすような 行列 があるとする：

レイリー・リッツ法は， の計算と関連する固有問題 を解くことからなる．

もとの問題 ， の近似固有ペアはそれぞれリッツ値およびリッツベクトルと呼ばれる および を満足する．

この方法は部分空間 が の不変量部分空間を近似するときに最も効果がある．

この方法は， が行列 および与えられた初期ベクトル に関連付けられたクリロフ部分空間に等しい場合に効果的である：

解説

アーノルディ法はクリロフ部分空間の直交基底を計算し， のときの投影された 行列 を生成するクリロフベースの投影アルゴリズムである．

入力：行列 ，ステップ数 ，ノルム1の初期ベクトル

出力：であり

アーノルディ法の場合， は簡約されないヘッセンベルグ形式（追加の非零部分対角要素を含む上三角）を持つ．

正規直交化は通常グラム・シュミット法で実行される．

このアルゴリズムで計算された数量は以下を満足する：

残差 は不変量部分空間への近さを示し，関連付けられたノルム は計算されたリッツペアの正確さを示す：

再出発

リッツペアは，初期ベクトル が所望の固有値の方向に富んでいれば迅速に収束する．

そうでない場合，作業およびメモリが過剰になるのを避けるために再出発が必要になる．

再出発には，特に以下のようないくつかのストラテジーがある：

実装

多数のソフトウェア実装が利用できる：

自動ストラテジー

"Automatic"設定では，下のようなメソッドの融合が使われる．

ステップの却下

評価時間のキャッシュは，優勢固有値推定が却下されたステップに対して再計算されないことを確実にする．

硬さ検出は以下の理由で，却下されたステップに対しても実行される．

反復メソッドオプション

"NonstiffTest"の反復メソッドには，変更することのできるオプションがある：

デフォルトの誤差許容範囲はただ1つの正しい数値を目標としているが，特に連続するステップで数回の反復が成功した後，実質的にもっと多くの正確な値を得ることがある．

反復の回数を制限するデフォルトの値は以下のようなものである：

レイテンシと切換え

"StiffnessSwitching"メソッドが継続的に硬いメソッドと硬くないメソッドとを切り換えようとするサイクルを避けるために，ある形でレイテンシを含ませることが重要である．

"StiffnessTest"および"NonstiffTest"のオプション"MaxRepetitions"，"SafetyFactor"はこの目的で使用される．

デフォルト設定では，切換えが非常に反応が早い．これは1ステップの積分メソッドに適している．

例題

ヴァンデルポール

例題の系を選ぶ：

StiffnessTest

この系は，与えれたメソッドと"StiffnessTest"のデフォルトオプション設定での積分が成功する：
積分が長くなると放棄され，硬さ検出条件が満足されないとメッセージが出力される：
単位安全因子を使い，硬さの失敗が1回だけ許されることを指定すると，厳密な検定が効果的に与えられる．この指定ではネストされたメソッドオプションのシンタックスが使われている：

NonstiffTest

小さい系の場合は，直接の固有値計算が使われる．

この例題は，全体的な硬さ切換えフレームワークが期待通りに動いているというより検定の例である．

硬いメソッドと硬くないメソッドとの切換えおよび取られるステップを監視する関数を設定する．硬いソルバと硬くないソルバのデータは，"Sow"に対して異なるタグを使うことで別々のリストに入れられる：
系を解き，メソッド切換えのデータを収集する：
明示的ソルバ（青）と暗示的ソルバ（赤）を使って，取られる刻み幅をプロットする：
取られた硬くないステップおよび硬いステップの数（却下されたステップを含む）を計算する：

CUSP

神経インパルスメカニズムに対するcusp catastrophyモデル[Z72]：

ヴァンデルポールと組み合せるとCUSP系になる[HW96]：

ここで

であり，および である．

線の方法を使った拡散項の離散化が，次元の常微分方程式系を得るために使われる．

ヴァンデルポール系とは異なり，問題の大きさのため，固有値推定には反復メソッドが使われる．

刻み幅と次数選択

解く問題を選ぶ：
使用されるメソッドのタイプと刻み幅を監視する関数を設定する．さらにメソッドの次数がツールチップとして含まれる：
積分が進むにつれて，メソッドの次数のデータを収集する：
明示的ソルバ（青）と暗示的ソルバ（赤）を使って，取られる刻み幅をプロットする．各ステップでのメソッドの次数は，ツールチップで示される：
取られる硬くないステップと硬いステップの合計数（却下されたステップを含む）を計算する：

ヤコビアンの例題

最初の2，3のヤコビアン行列を集める関数を定義する：

硬いメソッドへの切換えは0.00113425付近で起きており，硬くないことの最初の検定は次のステップ で起きている．

ヤコビアン のグラフィックスによる表示である：
計算する関数を定義し， , , … の最初のいくつかの固有値とノルム限界を計算する関数を定義する：

この例題では，ノルム限界は非常に急である．

Korteweg–deVries

Korteweg–deVriesの偏微分方程式は浅い水面の波動の数学的モデルである：

次の境界条件

を考え，区間 上で解く．

線の方法を使った離散化は192の常微分方程式系を形成するために使われる．

刻み幅

解く問題を選ぶ：
LSODAで使用されるBackward Differentiation Formulaメソッドでは，この問題を解く際に困難が生じる：
このプロットで刻み幅が急激に減少しているのが分かる：
これとは対照的に，StiffnessSwitchingはほとんど積分ステップを必要としない線形後退オイラー法を使うよう，直ちに切り替える：

積分のはじめに硬いソルバが選ばれたら，外挿法は決して硬くないソルバに切り換えない．

従って，これは硬くないことの検出で最悪の場合の例だといえる．

それでも部分空間反復を使用するコストは積分全体の時間のほんの数パーセントに過ぎない．

硬くないメソッドへの切換えを無効にして時間を計算する：

ヤコビアンの例

前に定義した監視関数を使って最初のいくつかのヤコビアン行列のデータを集める：
初期ヤコビアン のグラフによる表示である：
最初のいくつかの , , …の固有値およびノルム境界を計算し表示する：

ノルム境界はわずかに過大評価であるが，より重要なのは，実部および虚部の相対的な大きさにてついて言及していないという点である．

オプションの要約

StiffnessTest

オプション名
デフォルト値
"MaxRepetitions"{3,5}硬さ検定(6)の失敗が許される連続および合計回数を指定する
"SafetyFactor"硬さ検定(7)の右辺と左辺に使う安全係数を指定する

メソッドオプション"StiffnessTest"のオプション

NonstiffTest

オプション名
デフォルト値
"MaxRepetitions"{2,∞}硬さ検定(8)の失敗が許される連続および合計回数を指定する
"SafetyFactor"硬さ検定(9)の右辺と左辺に使う安全係数を指定する

メソッドオプション"NonstiffTest"のオプション