数学 と 数値解析 において 、 常微分方程式の数値解法 （ 数値積分 の特殊なケースを含む） の一部の手法では、 手法の誤差を制御し、 A安定性 などの 安定性特性 を保証するために、適応ステップサイズが使用されています。導関数のサイズに大きな変動がある場合、適応ステップサイズの使用は特に重要です。たとえば、地球の周りの衛星の動きを標準的な ケプラー軌道としてモデル化する場合、 オイラー法 などの固定時間ステップ法で 十分な場合があります。しかし、三体問題のように地球と月の両方を考慮して宇宙船の動きをモデル化する場合は、状況が難しくなります 。その場合、宇宙船が地球と月から遠いときには大きな時間ステップをとれるが、宇宙船が惑星の1つに衝突しそうになると、小さな時間ステップが必要になるというシナリオが生まれます。 ロンバーグ 法 と ルンゲ・クッタ・フェールベルグ法 は、適応ステップサイズを使用する数値積分法の例です。
例 簡単にするために、次の例では最も単純な積分法である オイラー法を使用しますが、実際には、優れた収束性と安定性の特性を持つ ルンゲ・クッタ 法などの高次法 が好まれます。
初期値問題を考える
ええ ′
（ t ） = ふ （ t 、 ええ （ t ） ） 、 ええ （ 1つの ） =
ええ
1つの
{\displaystyle y'(t)=f(t,y(t)),\qquad y(a)=y_{a}}
ここで、 y と f は ベクトルを表す場合があります (その場合、この方程式は複数の変数の結合された常微分方程式のシステムを表します)。
関数 f ( t , y ) と初期条件 ( a , ya ) が与えられており、 t = b における解を求めます 。 y ( b ) が b における正確な解を表し 、 y b が計算した解を表すものとします。 と書きます。 ここで は 数値解の誤差です。
ええ
b
+ ε = ええ （ b ）
{\displaystyle y_{b}+\varepsilon =y(b)}
ε
{\displaystyle \epsilon }
t の値のシーケンス ( t n )に対して 、 t n = a + nh 、オイラー法は対応する y ( t n ) の値を次のように
近似します。
ええ
ん + 1
（ 0 ）
=
ええ
ん
+ h ふ （
t
ん
、
ええ
ん
）
{\displaystyle y_{n+1}^{(0)}=y_{n}+hf(t_{n},y_{n})}
この近似の局所的打ち切り誤差は次のように定義される。
τ
ん + 1
（ 0 ）
= ええ （
t
ん + 1
） −
ええ
ん + 1
（ 0 ）
{\displaystyle \tau _{n+1}^{(0)}=y(t_{n+1})-y_{n+1}^{(0)}}
テイラーの定理 によれば、（ f が十分に滑らかであれば）局所的な打ち切り誤差はステップサイズの2乗に比例すること
が示される。
τ
ん + 1
（ 0 ）
= c
h
2
{\displaystyle \tau_{n+1}^{(0)}=ch^{2}}
ここで c は 比例定数です。
この解決策とそのエラーを でマークしました 。
（ 0 ）
{\displaystyle (0)}
c の値は 不明です。ここで、異なるステップ サイズでオイラー法を再度適用し、 y ( t n +1 ) の 2 番目の近似値を生成してみましょう。2 番目の解が得られ、これを a とラベル付けし ます。新しいステップ サイズを元のステップ サイズの半分にして、オイラー法を 2 ステップ適用します。この 2 番目の解は、おそらくより正確です。オイラー法を 2 回適用する必要があるため、局所誤差は (最悪の場合) 元の誤差の 2 倍になります。
（ 1 ）
{\displaystyle (1)}
ええ
ん +
1 2
=
ええ
ん
+
h 2
ふ （
t
ん
、
ええ
ん
）
{\displaystyle y_{n+{\frac {1}{2}}}=y_{n}+{\frac {h}{2}}f(t_{n},y_{n})}
ええ
ん + 1
（ 1 ）
=
ええ
ん +
1 2
+
h 2
ふ （
t
ん +
1 2
、
ええ
ん +
1 2
）
{\displaystyle y_{n+1}^{(1)}=y_{n+{\frac {1}{2}}}+{\frac {h}{2}}f(t_{n+{\frac {1}{2}}},y_{n+{\frac {1}{2}}})}
τ
ん + 1
（ 1 ）
= c
（
h 2
）
2
+ c
（
h 2
）
2
= 2 c
（
h 2
）
2
=
1 2
c
h
2
=
1 2
τ
ん + 1
（ 0 ）
{\displaystyle \tau _{n+1}^{(1)}=c\left({\frac {h}{2}}\right)^{2}+c\left({\frac {h}{2}}\right)^{2}=2c\left({\frac {h}{2}}\right)^{2}={\frac {1}{2}}ch^{2}={\frac {1}{2}}\tau _{n+1}^{(0)}}
ええ
ん + 1
（ 1 ）
+
τ
ん + 1
（ 1 ）
= ええ （ t + h ）
{\displaystyle y_{n+1}^{(1)}+\tau _{n+1}^{(1)}=y(t+h)}
ここで、誤差係数は 区間 にわたって一定であると仮定します 。実際には、その変化率は に比例します 。解を減算すると、誤差の推定値が得られます。
c
{\displaystyle c}
[ t 、 t + h ]
{\displaystyle [t,t+h]}
ええ
（ 3 ）
（ t ）
{\displaystyle y^{(3)}(t)}
ええ
ん + 1
（ 1 ）
−
ええ
ん + 1
（ 0 ）
=
τ
ん + 1
（ 1 ）
{\displaystyle y_{n+1}^{(1)}-y_{n+1}^{(0)}=\tau _{n+1}^{(1)}}
このローカルエラー推定は 3 次精度です。
ローカルエラー推定は、望ましい精度を達成するためにステップサイズをどのように変更するかを決定するために使用できます 。たとえば、ローカル許容値が許容される場合、 h を 次のように変化させることができます 。
h
{\displaystyle h}
トル
{\displaystyle {\text{tol}}}
h → 0.9 × h × 分
（
最大
（
（
トル
2
|
τ
ん + 1
（ 1 ）
|
）
1
/
2
、 0.3
）
、 2
）
{\displaystyle h\rightarrow 0.9\times h\times \min \left(\max \left(\left({\frac {\text{tol}}{2\left|\tau _{n+1}^{(1)}\right|}}\right)^{1/2},0.3\right),2\right)}
は 、次回の試行で成功を保証するための安全係数です。最小値と最大値は、前回のステップサイズからの極端な変化を防ぐためのものです。これにより、原則として、 次回の試行で約 の誤差が生じるはずです。 の場合 、ステップは成功したとみなされ、誤差推定値を使用して解が改善されます。
0.9
{\displaystyle 0.9}
0.9 ×
トル
{\displaystyle 0.9\times {\text{tol}}}
|
τ
ん + 1
（ 1 ）
|
<
トル
{\displaystyle |\tau _{n+1}^{(1)}|<{\text{tol}}}
ええ
ん + 1
（ 2 ）
=
ええ
ん + 1
（ 1 ）
+
τ
ん + 1
（ 1 ）
{\displaystyle y_{n+1}^{(2)}=y_{n+1}^{(1)}+\tau _{n+1}^{(1)}}
このソリューションは、実際にはローカル スコープでは 3 次 精度 (グローバル スコープでは 2 次精度) ですが、エラー推定がないため、ステップ数を減らすのに役立ちません。この手法は、 リチャードソン外挿 と呼ばれます。
この理論では、 初期ステップサイズ から始めて、局所的な誤差許容度を与えられた最適なステップ数を使用して、点 から までの常微分方程式の制御可能な積分を容易に 行うことができます。欠点は、特に低次 オイラー法 を 使用する場合に、ステップサイズが法外に小さくなる可能性があることです。
h = b − 1つの
{\displaystyle h=ba}
1つの
{\displaystyle a}
b
{\displaystyle b}
同様の方法は、4 次ルンゲ・クッタ法などの高次法にも適用できます。また、ローカル エラーをグローバル スコープにスケーリングすることで、グローバル エラー許容度を実現できます。
埋め込みエラー推定 いわゆる「埋め込み」誤差推定を使用する適応ステップサイズ法には、 Bogacki–Shampine 法 、 Runge–Kutta–Fehlberg 法 、 Cash–Karp 法、および Dormand–Prince 法があります。これらの方法は計算効率が高いと考えられていますが、誤差推定の精度は低くなります。
埋め込み法の考え方を説明するために、 を更新する次のスキームを考えます 。
ええ
ん
{\displaystyle y_{n}}
ええ
ん + 1
=
ええ
ん
+
h
ん
ψ （
t
ん
、
ええ
ん
、
h
ん
）
{\displaystyle y_{n+1}=y_{n}+h_{n}\psi (t_{n},y_{n},h_{n})}
t
ん + 1
=
t
ん
+
h
ん
{\displaystyle t_{n+1}=t_{n}+h_{n}}
次のステップは 前回の情報から予測されます 。
h
ん
{\displaystyle h_{n}}
h
ん
= グ （
t
ん
、
ええ
ん
、
h
ん − 1
）
{\displaystyle h_{n}=g(t_{n},y_{n},h_{n-1})}
埋め込みRK法の場合、の計算には 低次のRK法が含まれます 。その場合、誤差は次のように簡単に記述できます。
ψ
{\displaystyle \psi}
ψ 〜
{\displaystyle {\チルダ {\psi }}}
間違い
ん
（ h ） =
ええ 〜
ん + 1
−
ええ
ん + 1
= h （
ψ 〜
（
t
ん
、
ええ
ん
、
h
ん
） − ψ （
t
ん
、
ええ
ん
、
h
ん
） ）
{\displaystyle {\textrm {err}}_{n}(h)={\チルダ {y}}_{n+1}-y_{n+1}=h({\チルダ {\psi }}( t_{n},y_{n},h_{n})-\psi (t_{n},y_{n},h_{n}))}
間違い
ん
{\displaystyle {\textrm {err}}_{n}}
は正規化されていない誤差です。これを正規化するには、絶対許容値と相対許容値で構成されるユーザー定義の許容値と比較します。
トル
ん
=
アトル
+
ルトル
⋅ 最大 （
|
ええ
ん
|
、
|
ええ
ん − 1
|
）
{\displaystyle {\textrm {tol}}_{n}={\textrm {Atol}}+{\textrm {Rtol}}\cdot \max(|y_{n}|,|y_{n-1}| )}
え
ん
=
規範
（
間違い
ん
/
トル
ん
）
{\displaystyle E_{n}={\textrm {norm}}({\textrm {err}}_{n}/{\textrm {tol}}_{n})}
次に、正規化された誤差を 1 と比較して予測値を取得します 。
え
ん
{\displaystyle E_{n}}
h
ん
{\displaystyle h_{n}}
h
ん
=
h
ん − 1
（ 1
/
え
ん
）
1
/
（ q + 1 ）
{\displaystyle h_{n}=h_{n-1}(1/E_{n})^{1/(q+1)}}
パラメータ q は、RK 法に対応する次数であり 、RK 法の方が次数が低い。上記の予測式は、推定された局所誤差が許容値よりも小さい場合はステップを拡大し、そうでない場合はステップを縮小するという意味で妥当である。
ψ 〜
{\displaystyle {\チルダ {\psi }}}
上記の説明は、明示的RKソルバーのステップサイズ制御で使用される簡略化された手順です。より詳細な説明は、Hairerの教科書に記載されています。 [1] 多くのプログラミング言語のODEソルバーは、この手順を適応ステップサイズ制御のデフォルト戦略として使用し、システムの安定性を高めるために他のエンジニアリングパラメータを追加します。
参照
参考文献
^ E. Hairer、SP Norsett G. Wanner、「常微分方程式の解法 I: 非剛性問題」、第 II 節。
さらに読む William H. Press、Saul A. Teukolsky、William T. Vetterling、Brian P. Flannery、 『Numerical Recipes in C』 、第 2 版、CAMBRIDGE UNIVERSITY PRESS、1992 年 。ISBN 0-521-43108-5 ケンドール・E・アトキンソン 『数値解析』第2 版 、ジョン・ワイリー・アンド・サンズ、1989年。ISBN 0-471-62489-6 
![{\displaystyle [t,t+h]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fe3196e4e40380aeccb623b8fec709179e935da6)
