Parameters: Transient Solver
モデル要素Param_Transientは、時間領域に基づいた非線形動解析のシミュレーション制御パラメータを定義します。
説明
- 積分器の選択。
- 積分誤差トレランス(または解の精度)。
- シミュレーション時の積分器のパフォーマンス。
動的シミュレーションは、自由度が1以上のシステムに対して実行されます。動的シミュレーションでは、すべての慣性の影響、すべての適用される力、および内部拘束が考慮されます。これによって、複雑な機構システムの正確なシステムレベルのシミュレーションが実行可能となります。
- 微分代数方程式(DAE)形式の方程式。これは、多くの場合、より正確で、堅牢で、効率的なシミュレーションにつながるMotionSolve内のデフォルト定式化です。
- 常微分方程式(ODE)形式の方程式。これは、古い定式化で、古いモデルとの下位互換性をサポートするために含まれています。
フォーマット
<Param_Transient
[ integrator_type = "string" ]
[ integr_tol = "real" ]
[ h_max = "real" ]
[ h_min = "real" ]
[ h0_max = "real" ]
[ vel_tol_factor = "real" ]
[ rel_abs_tol_ratio = "real" ]
[ max_order = "integer"]
[ dae_alg_tol_factor = "real" ]
[ dae_index = "integer"]
[ dae_constr_tol = "real" ]
[ dae_corrector_maxit = "real" ]
[ dae_corrector_minit = "real" ]
[ dae_vel_ctrl = "TRUE|FALSE" ]
[ dae_jacob_eval = "integer" ]
[ dae_eval_expiry = "integer" ]
[ dae_jacob_init = "integer" ]
[ dae_interpolation = "TRUE|FALSE" ]
[ newmark_beta = "real" ]
[ newmark_gamma = "real" ]
[ hht_alpha = "real" ]
/>
属性
- integrator_type
- ODE形式の運動方程式で使用される積分器を定義します。次のいずれかを選択します:
- integr_tol
- 積分器が変位、速度、微分方程式の状態を計算する際に許容されるステップあたりの最大絶対誤差を表します。
- h_max
- 積分器で許容される最大ステップサイズを定義します。このパラメータを使用して、結果の精度を制御することもできます。一般に、時間ステップが小さいほど、結果の精度が高まります。非連続性があるモデルの場合(接触など)、h_maxに小さい値を指定する必要があります。
- h_min
- 積分器で許容される最小ステップサイズを定義します。
- h0_max
- 最大初期ステップサイズ。
- vel_tol_factor
- 速度状態の誤差トレランスを得るため、integr_tolと掛け合わせる係数。
- rel_abs_tol_ratio
- 積分器の相対誤差トレランスを得るため、integr_tolと掛け合わせる係数。
- max_order
- 積分器で使用される最高次数を指定します。各積分器には、独自の最高次数が組み込まれています。次数が高いほど精度が高まりますが、数値的手法での安定性は低下します。上の表に、MotionSolve内のさまざまな積分器の最大次数を示します。接触などの非連続イベントを含むモデルの場合は、max_order = 2に設定することをお勧めします。
属性(DAEソルバーに固有)
- dae_alg_tol_factor
- Lagrange乗数のような代数状態の誤差トレランスを得るため、integr_tolと掛け合わせる係数。
- dae_index
- DAE定式化のインデックス。
- dae_constr_tol
- 収束時に修正子が満たす必要のあるすべての代数拘束方程式のトレランス。
- dae_corrector_maxit
- 収束を達成するために修正子が実行できる反復計算の最大回数。
- dae_corrector_minit
- 修正子が実行する反復計算の最少回数。この回数に達すると、修正子の発散が生じていないかチェックされます。
- dae_vel_ctrl
- 各ステップで局所積分誤差が生じていないかを確認するために速度状態をチェックするかどうかを制御する論理フラグ。これはSI1積分器にのみ有効です(dae_indexをご参照ください)。
- “TRUE”の場合、積分器は、各ステップの最後に実行される局所積分の有無のチェックに速度状態を含めます。
- “FALSE”の場合、積分器は、各ステップの最後に実行される局所積分の有無のチェックに速度状態を含めません。
- dae_jacob_eval
- この属性は、修正子の反復計算時にヤコビアンマトリクスの評価の頻度を制御します。
- 0
- 積分器は収束速度を調べることで、新しいヤコビアンが必要なタイミングを自動的に特定します。
- 1
- 反復計算ごとに新しいヤコビアンが計算されます。
- 2
- 1つおきの反復計算で新しいヤコビアンが計算されます(例: 反復計算1-3-5など)。
- 3
- 2つおきの反復計算で新しいヤコビアンが計算されます(例: 反復計算1-4-7など)。
- dae_eval_expiry
- これはオプションの属性です。積分ステップの数を指定します。この数の積分ステップの後は、dae_jacob_evalで定義された評価パターンが無視され、デフォルトの評価パターンが使用されます。修正子の収束速度を調べることで、新しいヤコビアンが必要になると、デフォルトの評価パターンが自動的に特定されます。
- dae_jacob_init
- この属性は、修正子の反復計算時の初期ヤコビアンマトリクスの評価数を指定します。
- 0
- 積分器は収束速度を調べることで、新しいヤコビアンが必要なタイミングを自動的に特定します。
- 1
- 最初の反復計算で新しいヤコビアンが計算され、その後は計算されません。
- 2
- 最初と2番目の反復計算で新しいヤコビアンが計算され、その後は計算されません。ここに3以上の値を設定した場合も同様です。
- dae_interpolation
- 出力ステップでの結果に対して積分器が補間を使用するかどうかを指定します。
- TRUEは、ユーザーが要求した出力ポイントで、MotionSolveが(必要に応じて)結果を補間することを意味します。
- FALSEは、ユーザーが要求した出力ポイントで、MotionSolveがDASPK積分器を使用して(必要に応じて)結果を補間することを意味します。
- newmark_beta
- Newmark-β法におけるβ係数を指定します。係数についての詳細は、コメントの参考文献をご覧ください。
- newmark_gamma
- Newmark-β法におけるγ係数を指定します。係数についての詳細は、コメントの参考文献をご覧ください。
- hht_alpha
- Hilber-Hughes-Taylor法におけるα係数を指定します。β係数とγ係数はα係数から導かれます。これらの係数についての詳細は、コメントの参考文献をご覧ください。
コメント
- 陽解法と陰解法での積分
陽解法および陰解法は、時間に依存する方程式の解の近似を求めるための数値的手法です。陽解法を使用した積分器は、tn時点の状態からtn+1時点のシステムの状態を計算します。この計算は1ステップで明示的に行うことができます。陰解法を使用した積分器は、tn+1とtnの両方の時間からの状態を含む方程式を解くことによって解を求めます。後者の場合、解を求めるための反復処理が必要となるため、各時間ステップの計算コストが高くなります。また、アプリケーションによっては、システムに数値的な減衰をもたらすこともあります。しかし、陰解法を使用した積分器は通常、より大きなステップサイズ Δt = tn+1 – tn,許容し、数値的に硬い式に適しており(コメント 2参照)、微分代数方程式を解くことができます(コメント 3参照)。
- 数値的剛性と物理的剛性
数値的剛性は、“物理的”剛性とはかなり異なります。数値的剛性は、システムの減衰特性によって決まりますが、物理的特性は、システムに内在する剛性(システム内のスプリングなど)によって決まります。数値的剛性の概念は、作用点でのシステムの固有値に密接に関係しています。
システムの固有値には、実部と虚部の2つの成分があります。実部はシステムの減衰特性を定義し、虚部は作用点でのシステムの振動周波数を定義します。- 虚部成分のない固有値は、過減衰振動モードを定義します。
- 実部成分のない固有値は、不減衰振動モードを定義します。
- 複素固有値は、不足減衰振動モードを定義します。
数値的に硬いシステムは、次のプロパティを示す、非常に離れた固有値によって表現されます。- 高周波成分は過減衰で、急速に減衰します。
- 解の低周波成分は不足減衰の挙動を示します。
数値的剛性の概念は、以下に示す簡単な例を使用して説明されます。
下の図は、2つの結合されていないスプリング質量システムを示しています。最初のシステムでは、質量M1がプロパティK1、C1、およびL01を持つスプリングで支持されています。2番目のシステムでは、質量M2がプロパティK2、C2、およびL02を持つスプリングで支持されています。全体座標系を基準にして測定された質量M1の変位は、x1と表されています。全体座標系を基準にして測定された質量M2の変位は、x2と表されています。
以下を前提とします:- M1=1Kg、C1=2N sec/m、K1=104 N/m
- M2=1Kg、C2=2x104 N sec/m、K2=108 N/m
図 1. 結合されていない2つの単純なスプリング質量システム
各システムiにおいて、以下のことが容易にわかります:- システムiの運動方程式:
- システムiの固有値:
- 最初のシステムの固有値と応答:
- 2番目のシステムの固有値と応答:
システム1は不足減衰です。初期外乱に応じて振動の減衰が示されます。0.01秒後に、初期外乱は元の値のe-0.01 = 0.99倍に減少します。システム1は、99.995ラジアン / 秒の周波数で振動します。
一方、システム2は過減衰です。0.01秒後に、初期外乱は元の値のe-100 = 3.72 x 10-44倍に減少します。どのような初期外乱であっても、システム応答はほぼ瞬間的になくなります。システム2が振動することはありません。
このような特性を表すシステムを数値的に硬いと言います。解のある成分が初期励起に対してゆっくりした減衰を示すのに対し、別の成分はほぼ瞬間的に消えます。
多くの機構システムは、数値的に硬いプロパティを有しているか、微分代数方程式で表現されます(コメント 3を参照)。一般に、これらは陰解法を使ってのみ効率的に解くことができます。
- 微分代数方程式常微分方程式(ODE)と代数拘束方程式を組み合わせた連立方程式は、微分代数方程式(DAE)と呼ばれます。
図 2. マルチボディシステム。
例えば、マルチボディシステム(MBS)で、ニュートンの第二法則は2次の常微分方程式(ODE)で表され、このシステム内の拘束(ジョイントなど)は代数方程式(AE)で表されます。したがって、通常のマルチボディシステムの数式的表現は、連立微分代数方程式(DAE)になります。上の図に示されたメカニズムは、通常、微分代数方程式で表現されます。
常微分方程式は次のように表されます:y は、積分される状態量です。
微分代数方程式は次のように表されます: - DAEのインデックス連立DAEのインデックスは、連立常微分方程式を取得するために、独立変数(通常は時間)に関してシステム内の拘束を微分する回数として定義されます。これは、下の図に示すとおりです。
図 3. 連立DAEのインデックス
- 方程式定式化MotionSolveは、次の3種類の方程式の定式化を動解析用にサポートしています:
- ODE定式化
- ODE定式化は、stiff積分器とnon-stiff積分器の両方をサポートします。MotionSolveは、coordinate partitioning法を用いてDAE形式の運動方程式をODE形式に変換し、そのうえで、ODE積分器を使って結果の方程式を解きます。stiff積分器とnon-stiff積分器の両方がサポートされています。この定式化によってサポートされるstiff積分器は、(a) VSTIFFと(b) MSTIFFです。non-stiffソリューションを積分するには、積分器ABAMが使用されます。
- DAE Index-3 (I3)定式化
- I3定式化は、DAE形式の運動方程式を、この形式の運動方程式を解くことのできる積分器に提供します。DSTIFF、CSTIFF、NSTIFF、HSTIFFはI3形式の運動方程式を解くことができます。I3定式化では、積分器は、速度の状態での積分のローカルエラーを監視しません。したがって、I3ソリューションは通常、非常に高速ですが、速度に関して多少正確さに欠ける傾向があります。
- DAE Stabilized Index-1 (SI1)定式化
- SI1定式化は、さらに“安定化させた”index-1 DAE形式の運動方程式を、この形式の運動方程式を解くことのできる積分器に提供します。“安定化させた”とは、I3セットに追加される速度と加速度の拘束の余剰と思われるセットに対して述べている言葉です。これらは、速度および加速度の、速度および加速度の拘束との一貫性を保ちます。解が正確であれば、これは自動的に行われます。DSTIFF、CSTIFF、NSTIFF、HSTIFFはSI1形式の運動方程式を解くことができます。SI1定式化では、積分器は、変位と速度の両方の状態での積分のローカルエラーを監視します。したがって、SI1ソリューションは通常、非常に正確であると言えます。SI2ソリューションのスピードはI3ソリューションのスピードに比較すると若干遅くなります。
- Index-3 DAE形式の運動方程式システムが次のように定義されているものとします:
- “N”個の変位座標q=[q1, q2, … qN]Tを使用してシステムを表します。
- 座標qは、すべて独立しているわけではありません。
- “M”個の変位拘束Φ =[ Φ1(q,t), Φ2(q,t), … ΦM(q,t)]Tが、“N”個の座標qの間に存在します。
- λ = [ λ1, λ2 … λM]Tが、拘束Φを維持するために必要なLagrange乗数を表すとします。
- F = [F1, F2 … FNa]Tが、システム上で作用する外力を表すとします。
- システムの運動エネルギーは、
、位置エネルギーは、
とします。
- Lagrange関数
は、
と表されます。
これらの仮定の下で、システムの運動方程式は
次のように表現できます:
上記の方程式は、運動方程式のEuler-Lagrange表現と呼ばれます。マトリクス形式では、式(1)は次のように記述できます:式(2.1)は、力の釣り合いの方程式を表します。最初の項は慣性力、2番目の項は拘束力、3番目の項は一般化外力を表します。式(2.1)は、したがって、ニュートンの運動の第二法則を書き換えたものとなります。式(2.1)は、速度(加速度)の時間導関数と、システムに作用する外力および内力を関連付ける微分方程式です。
式(2.2)は、システム内の運動微分方程式を表します。これらの式は、1次形式の運動方程式を表します。これを行うには、新しい一連の速度の変数uを作成します。速度は、変位の1次時間導関数として定義されます。式(2.2)も微分方程式です。
式(2.3)は、システム内の代数拘束を表します。これらが成立するのは、これらの方程式が、独立座標セットについて定式化されなかったためです。システムにはp=N-Mの自由度があります。したがって、N個の座標を、M個の代数拘束方程式により、互いに関連付けることができます。
式(2.4)は、Control_Diff、Control_SISO、およびControl_StateEqn要素を使用して指定される、ユーザー定義の微分方程式を表します。ユーザー定義の状態は、xで表されます。
式(2)は、一連の陰的な微分代数方程式(DAE)です。Z = [u, q, λ, x]Tとします。式(2)は、圧縮表記で次のように表すことができます:式(2.1)を1回、式(2.3)を3回微分することで、式(3)を次の形式の一連の常微分方程式(ODE)に変換することができます:DAEのインデックスは、DAEを一連のODEに変換するために拘束を微分する回数として定義されます。式(4)で表されるODE形式に達するためには拘束(2.3)の3回微分が必要であるため、DAEのインデックスは3となります。Index-3 DAEは、以下のような特性もあります。式(2)または(3)の数値解法には、専用のDAE積分器が必要です。これらの手法は、{u}、{q}、および{λ}の精度はさまざまなレベルに制御されるという基本概念を共有しています。一般的なアプローチの1つ(Index-3)が、{u}および{λ}での局所積分誤差を無視し、{q}での局所積分誤差のみを監視するというものです。したがって、Index-3の解法は精度は低下しますが、非常に高速になる傾向があります。非常に正確な解法が必要な場合には、代替アプローチが推奨されます。
- 陰解法を使用した積分器を使用した運動方程式の数値積分それぞれの陰解法積分器は、各時間ステップで3つの重要な操作を実行する必要があります:
- 予測予測子の目的は、修正子に優れた初期推測値を提供することです。予測子は、以前計算した解析ポイントを介して多項式をフィットさせます。その後、その多項式を時間の新しい値に外挿し、以下のように解を“予測”します。
図 4. 予測子
- 修正
修正子の目的は、運動方程式を満たすよう、予測値を改善することです。この反復プロセスには、しばしばNewton-Raphson法などの手法が用いられます。
- 局所切り捨て誤差の推定
修正子の主な役割は、解が運動方程式を満たすことを強制することですが、全体的な精度を保証するものではありません。次に、ステップ内の積分誤差を計算し、精度を推定します。これは、補正された解が必ずしも"真の"解を表しているとは限らないからです。例えば、初期予測値が真の解から大きく外れている場合、修正子はもはやモデルの初期条件と一致しない解に収束している可能性があります。つまり、得られた結果は、運動方程式に忠実ではあるけれども、解の多様体の中で別の解へと導かれた可能性があります。
- 予測
- 積分器の失敗によりモデルが失敗します。この問題をどのように診断し、修正することができるでしょうか?
DSTIFFやその他のマルチステップ積分器は、各ステップで、解の予測値と修正された値の差を計算し、積分器で計算される“誤差係数”でそれをスケーリングすることで、局所的な誤差を推定します。したがって、不適切な予測により、積分誤差が生じる可能性があります。
予測が不正確になる理由として次の2つが考えられます:- システム内に不連続なモーション入力がある。
- システム内にほぼ不連続な力がある。
このようなモデリングの効果を以下に示します。図 5. 不連続性が積分の失敗の根本原因
上記のような不連続性は、IF-THEN-ELSEロジックを不用意に使用した関数表現やユーザーサブルーチンによって生じることがあります。不連続性や急な変化がないかを調べる必要のあるモデリング要素には、次のようなものがあります:- Forces
- 変数
- ユーザーの微分方程式
- 一般状態方程式
IF-THEN-ELSEロジックを滑らかなSTEP関数に置き換えることをお勧めします。
また、MotionSolveのXMLステートメント:
を使用して、詳細なニュートン-ラプソン反復の出力を有効にしてください。この出力を確認することで、通常は失敗の原因になりそうなモデリング要素まで戻ることができます。<DebugOutput> SwitchOn="True" </DebugOutput>
典型的な積分問題に対処するための標準的なアプローチには、以下のステップがあります:- 剛性と減衰係数が現実的な値に設定されていることを確認してください。
- 質量と慣性に物理的に意味のある値が割り当てられていることを確認します。
- モデリング要素の不連続性を特定し、修正します。
- 極端に急なステップ関数に注意してください。
- 一時的な措置として、積分誤差(INTEGR_TOL)を緩和することを検討し、根本的な問題の理解を深めてください。
- 結果は収束していますか?つまり、誤差トレランスと共に変化しなくなりましたか?収束スタディを実施して、積分器から取得した結果が変化しなくなっている(収束している)ことを確認することをお勧めします。積分器は正確な解を近似するため、積分器内で許容される誤差が多すぎる場合、誤差トレランスが厳しい(小さい)場合とは異なる結果を積分器が提供することがあり得ます。必要に応じて積分器を調整できるようにする設定はある程度ありますが、最初に調査を開始するために使用できるいくつかの設定を以下に示します:
- integr_tol
- h_max
- DSTIFF用のvel_tol_factor。これには、dae_indexパラメータを経由してSI1が必要です。
収束した結果はさまざまな方法で取得できますが、ユーザーに有効な一例を以下に示します:- モデル内の注目するデータのチャネルを選択します(変位、力など)。
- 積分器を選択します。 - その積分器で何の誤差制御を行いますか?(変位、速度、微分状態)
- 積分器が収束を制御する状態まで、誤差を低減します。収束結果を得ることができる最大誤差を取得します。
- まだ積分器が誤差トレランスを使用して制御できない状態に注目している場合は、そのような状態が収束するまで積分器の最大ステップサイズ(h_max)を下げ、収束結果を得ることができる最大のh_maxを使用します。
- シミュレーションが失敗した場合、誤差トレランスを緩和(大きく)して積分器の最大ステップサイズ(h_max)の精度を制御し、問題のシミュレーションを成功させます。これは可変ステップの積分器より時間がかかるため、最後の手段として使用します。
- その他多くの積分器設定も使用可能です。自分のモデルに適用できるその他のオプションも使用してみてください。
- さまざまな定式化アプローチの相対比較以下の表は、さまざまなDAEアプローチの相対的な利点をまとめたものです。さまざまなアプローチを特徴付けるにあたり、次の基準を使用しています:
- 精度
- 解の配列で累積される積分誤差を測定します。
- 速度
- さまざまなアプローチで共通する相対CPU時間を測定します。これはモデルに大きく依存し、さまざまなアプローチをテストする際のこれまでの経験が反映されます。
- 不連続性の処理
- 定式化がシステム入力(モーションまたは力)の尖った角や不連続性を処理する能力を測定します。不連続性や尖った角(非微分可能性)は、モデルで回避する必要があります。
- 一貫性
- 解が拘束Φ、その1次および2次時間導関数をいかに満たしているかを測定します。
- 高周波の追跡
- システム内の高周波がいかに追跡されているかを測定します。
- 数値減衰
- 特定のインデックス定式化の実装により適用された数値減衰を測定します。一般に、数値減衰が大きい場合、精度は低下しますが、高速でより堅牢な解が得られます。
- 堅牢性
- デフォルトのDAE定式化はIndex-3です。最初にIndex-3定式化を使用し、上記の基準を満たさない解が見つかった場合にのみ別の定式化への切り替えを検討することをお勧めします。
注: 次の表は、主観的な情報を含んでおり、絶対的な推奨事項ではなく、一般的なガイドラインとして使用してください。(*ABAMの場合、硬くない問題に対しては高速、硬い問題に対しては非常に低速になります。)表 1. さまざまな方程式定式化の相対的な利点 I-3 SI-1 ODE (SI-0) ABAM 精度 中 高 高 高 速度 高 中 低 高 / 非常に低い*
不連続性の処理 高 低 低 高 一貫性 低 高 高 高 高周波の追跡 低 高 高 高 数値減衰 高 低 低 低 堅牢性 高 中 中 高 - DAE_CONSTR_TOLを0.001 mmより大きくすることは、一般的にお勧めできません。このトレランスが、モデルの長さ単位で指定されたシステム内の運動拘束の最低限の解析精度となります。DAE_CONSTR_TOLが大きい(0.001mmより大きい)場合、システム拘束は緩く満たされるだけになります。不正確な応答が作成されるだけでなく、拘束充足における誤差により数値的な解が不安点になり、失敗につながることもあります。
- dae_jacob_evalとdae_jacob_init はどちらも、修正子の反復計算時にヤコビアンを評価するタイミングを制御します。この2つの主な違いは、ヤコビアン評価の頻度に影響を与える方法です。Mをニュートン-ラプソン反復の反復カウンターとします。このインデックスが0インデックスとして実装されていることに注意してください。これは0から開始することを意味します。ヤコビアンは、次の条件のいずれかまたは両方が満たされるたびに再評価されます。
- MOD(M, dae_jacob_eval) = 0
- M < dae_jacob_init
これを以下の表にわかりやすく示します。
修正子反復カウンター(M) dae_jacob_eval = 2 dae_jacob_init = 2 新しいヤコビアン? MOD (M, 2) = 0? M < 2? 0 ○ ○ ○ 1 × ○ ○ 2 ○ × ○ 3 × × × 4 ○ × ○ 修正子反復カウンター(M) dae_jacob_eval = 3 dae_jacob_init = 2 新しいヤコビアン? MOD (M, 3) = 0? M < 2? 0 ○ ○ ○ 1 × ○ ○ 2 × × × 3 ○ × ○ 4 × × × 5 × × × 6 ○ × ○ 7 × × × - 参考資料DSTIFF:
- “The Numerical Solution of Initial Value Problems in Differential-Algebraic Equations”、L. Petzold、K. E. Brenan、およびS. L. Campbell、Elsevier Science Publishing Co.、(1989年)、2版、SIAM Classics Series、1996年。
- “Consistent Initial Condition Calculation for Differential-Algebraic Systems”、P. N. Brown、A. C. Hindmarsh、およびL. R. Petzold、SIAM J. Sci.Comp.19(1998年)、1495-1512ページ。
NSTIFF- "A method of computation for structural dynamics", Newmark, Nathan M., Journal of the Engineering Mechanics Division, 85 (EM3): 67–94, 1959.
HSTIFF- " Improved Numerical Dissipation for Time Integration Algorithms in Structural Dynamics", Hilber, H.M, Hughes,T.J.R and Talor, R.L. , Earthquake Engineering and Structural Dynamics, 5:282-292, 1977.
VSTIFF- “VODE, a Variable-Coefficient ODE Solver”、P. N. Brown、G. D. Byrne、およびA. C. Hindmarsh、SIAM J. Sci.Stat.Comput.、10、1038-1051ページ、1989年。
- “A Poly-algorithm for the Numerical Solution of Ordinary Differential Equations”、G. D. ByrneおよびA. C. Hindmarsh、ACM Trans.Math.Software、1(1975年)、71-96ページ。
- “EPISODE: An Effective Package for the Integration of Systems of Ordinary Differential Equations," A. C. Hindmarsh and G. D. Byrne, LLNL Report UCID-30112, Rev. 1, April 1977.
- Scientific Computingの“ODEPACK, a Systematized Collection of ODE Solvers”、 A. C. Hindmarsh、R. S. Stepleman他(編)、North-Holland、Amsterdam、1983年、55-64ページ。
MSTIFF- “The Integration Of Stiff Initial Value Problems In O.D.E.S Using Modified Extended Backward Differentiation Formulae”、J. Cash、Comp. and Maths. with Applics.、9:645-657、1983年。
- “The Integration Of Stiff Initial Value Problems In O.D.E.S Using Modified Extended Backward Differentiation Formulae”、J.R.Cash、Comp.And Maths.With Applics.、9、645-657、(1983年)。
- “Stable Recursions With Applications To The Numerical Solution Of Stiff Systems”、J.R. Cash, Academic Press, (1979).
- “Ordinary Differential Equation System Solver”、A.C.Hindmarsh、Gear, C.W.、Ucid-30001 Rev. 3、Lawrence Livermore Laboratory、Livermore、Ca 94550、1974年12月。
- “An Mebdf Code For Stiff Initial Value Problems”、J.R.CashおよびS. Considine、ACM Transactions On Mathematical Software、1992年。
ABAM- “Computer Solution of Ordinary Differential Equations”、Lawrence F. ShampineおよびMarilyn K. Gordon、W.H.Freeman & Co Ltd、1975年6月。