- 2025/12/22: プログラムを C# 14.0 に修正
こんにちは、Rafkaです。
今回は、前回紹介したオイラー法よりも誤差を抑えられる、
ホイン法（ 2 次のルンゲクッタ法とも呼ばれます）の導出と実装を行っていきます。
アルゴリズム
導出する前に、ホイン法のアルゴリズムを紹介します。
ちなみにホインの綴りは少し特殊で、Heunと書きます。
ホイン法では以下の更新式で
x(t) の値を更新していきます。
上の式の感覚的な解釈としては、
時間
t における傾き
k1 を計算し、
次の時間ステップまでをその傾きの１次関数として近似した場合の、
次のステップでの傾き
k2を計算し、
それらの平均を次の時間ステップへの値の更新に使用するといった感じです。
それでは導出の方に進みましょう。
導出
とりあえず
x(t+Δt) をテイラー展開します。
前回のオイラー法よりも精度を上げたいので、今回は 2 次の項まで展開します。
ここで厄介なのは 2 階微分の項です。
2 階微分の項は以下のように、
一度定義の式で
f(x(t),t) に置き換えた後、
f について微分を行うことで計算を進められます。
微分は合成関数の微分ですね。
これをテイラー展開した式の 2 階微分の項に代入すると、以下のようになります。
x(t+Δt)=x(t)+f(x(t),t)Δt+21(∂x∂f(x(t),t)f(x(t),t)+∂t∂f(x(t),t))Δt2+O(Δt3)
偏微分の項が出てきてしまったので、
未知定数
a,b を用いて表現された以下の式を組み合わせて、
上の式と同じ形に変形し、
それぞれを比較して、
定数を決定して最終的な更新式の導出を行いたいと思います。
後で使うので、先に
k2 を 1 次の項までテイラー展開しておきます。
k1 とテイラー展開した
k2 を 3 番目の更新式に代入して変形すると、
以下の式になります。
x(t+Δt)=x(t)+(a+b)f(x(t),t)Δt+b(∂x∂f(x(t),t)f(x(t),t)+∂t∂f(x(t),t))Δt2
同じ様な形になりましたね。
未知定数の部分の比較をすると、
a+b=1 、
b=1/2 であることが分かるので、
とすれば良いことが分かります。
これらの未知定数を更新式に代入して整理すると、
アルゴリズムの部分で示した以下の更新式が導かれます。
k1k2x(t+Δt)=f(x(t),t)=f(x(t)+k1Δt,t+Δt)=x(t)+2Δt(k1+k2)
途中示したテイラー展開の式から、ホイン法が計算精度 2 次の手法であることが分かります。
実装
オイラー法のメソッドと同様な形になるように、以下のように実装しました。
void Heun(ModelFunction f, Span<double> x, double t, double dt) {
var n = x.Length;
Span<double> k1 = stackalloc double[n];
Span<double> k2 = stackalloc double[n];
Span<double> tempX = stackalloc double[n];
f(x, t, k1);
ScaleAdd(x, k1, dt, tempX);
f(tempX, t + dt, k2);
for (var i = 0; i < n; i++) {
var slope = (k1[i] + k2[i]) * 0.5;
x[i] += slope * dt;
}
}
8 行目で
k1 を計算し、
9 行目で
x(t)+k1Δt を計算して
tempX に保存しています。
10 行目で
k2 を計算し、
最後のループで更新式に基づいて
x(t+Δt) を計算しています。
まとめ
今回はオイラー法よりも計算精度の良いホイン法の導出と実装を行いました。
次回はルンゲクッタ法の紹介を行いますが、
導出は書くと大変なので割愛させていただきます。
導出を行っている他のサイトさんを見つけたのでそちらの紹介を行い、
アルゴリズムの紹介とその実装をメインに説明したいと思います。
