Title : Math::RungeKutta
|
Location : Home > Languages > Perl > Package Title : Math::RungeKutta |
|
Math::RungeKutta - 微分方程式の積分
use Math::RungeKutta;
sub dydt { my ($t, @y) = @_; # 導関数
my @dydt; ... ; return @dydt;
}
@y = @initial_y; $t=0; $dt=.4; # 初期条件
# 自動タイムステップ調整
while ($t < $tfinal) {
($t, $dt, @y) = &rk4_auto(\@y, \&dydt, $t, $dt, 0.00001);
&display($t, @y);
}
# または固定タイムステップ
while ($t < $tfinal) {
($t, @y) = &rk4(\@y, \&dydt, $t, $dt); # Merson の4次の方法
&display($t, @y);
}
# 他の方法。しかし正確さが劣る
($t, @y) = &rk2(\@y, \&dydt, $t, $dt); # Heun の2次の方法
# また他の方法。ただしデフォルトでエクスポートされない
import qw(:ALL);
($t, @y) = &rk4_classical(\@y, \&dydt, $t, $dt); # Runge-Kutta 4次
($t, @y) = &rk4_ralston(\@y, \&dydt, $t, $dt); # Ralston 4次
RungeKutta.pm 同次微分方程式を数値積分するアルゴリズムを提供する。
dY/dt = F(t,Y)
ただし Y は既知の初期値 Y(0) を持つ変数の配列、F は問題から示される関数とする。
Runge-Kutta 法はタイムステップ内の様々な点で導関数 F(t,Y) を評価し、結果を結合して正確な値 Y(t+dt) を導出する。本モジュールは陽的 Runge-Kutta 法のみを利用している。陰的な方法は各タイムステップで Y(t) 及び F(t,Y) をともに含む同次方程式の集合を説く必要があり、一般に処理が難しい。
３つの主なアルゴリズムが提供されている。
rk2 は Heun の2次 Runge-Kutta アルゴリズムで、ややあいまいな面があるが、広い範囲での安定性という点では有用である。rk4 は Merson の4次 Runge-Kutta アルゴリズムで、指定されたステップサイズで通常の選択を行うものである。rk4_auto は自動的に指定された精度を確保するために rk4 のステップサイズを調整する step-doubling メソッドを利用している。これは適切なステップサイズを探索する手間を省いてくれ、解がゆっくり変わるときにはステップサイズを自動的に大きくして CPU タイムを節約する。
Perl は全地球気候シミュレーションや銀河の衝突と言った大規模な数値シミュレーションには不適な言語である（もしそれが必要なら xmds をチェックしよう）。しかし Gear が言うように「デジタルコンピュータで解かれる多くの問題は不十分な手法を与えられているという小さなことで区分される。コンピュータでの時間消費は少ない。プログラムを準備する人間の時間を最小化する方法が最良の方法だ。」
本モジュールは頑強で使いやすいように設計されており、Perl の文脈で起こりうる、経済・金融・人口学や生態学のモデリング・機械やプロセスでの力学問題で見られる微分方程式を解くことができる。
rk2( \@y, \&dydt, $t, $dt )
ただし引数は以下のとおり
\@y：変数の初期値の配列への参照
\&dydt：微分を計算する関数への参照
$t：初期時刻
$dt：タイムステップ
使用しているアルゴリズムは Ralston により導出されたもので、微分には Lotkin の境界条件を用いている。解の誤差を最小化する (gamma=3/4)。これは Heun メソッドとして知られているが、不幸にも同じメソッド名が他の分野で用いられている。各タイムステップで２回の関数評価が必要で、残余誤差は3次でタームごとに大きくなる。
rk2 は ($t, @y) を返す。ただし $t と @y はタイムステップ終了時には新しい値となっている。
rk4( \@y, \&dydt, $t, $dt )
引数は rk2 と同じである。
使用しているアルゴリズムは Merson により開発されたもので、タイムステップごとに５回の関数評価を行う。同一のタイムステップでは rk4 のほうが rk4_classical より正確である。余分な関数の評価を行っているおかげである。
rk4 は ($t, @y) を返す。ただし $t と @y はタイムステップ終了時には新しい値となっている。
rk4_auto( \@y, \&dydt, $t, $dt, $epsilon )
rk4_auto( \@y, \&dydt, $t, $dt, \@errors )
引数の１つめの形態は以下のとおり。
\@y：変数の初期値の配列への参照
\&dydt：微分を計算する関数への参照
$t：初期時刻
$dt：タイムステップ
$epsilon：$epsilon*$ymax の周りのステップごとの誤差
２つめの形態は
\@errors：最大許容誤差の配列への参照
１つめの $epsilon を呼び出す形態は @y の全要素が同一単位で同一の型（例えば y[10] が10歳から11歳の人口で y[25] は25歳から26歳の人口 など）の場合は有用である。4番目の引数のデフォルト値では $epsilon = 0.00001 である。
２つめの \@errors を呼び出す形態はその他の場合（例えば y[1] は国民総生産、y[2] は利子率 など）に有用である。この場合、許容誤差はそれぞれの変数に対して絶対値で特定され、スケールであわすことができない。
rk4_auto は要求された精度を保持するために自動的にタイムステップを調整する。これはあるタイムステップで計算し、その半分のタイムステップで計算し、その結果を比較することで行う。(rk4 で用いる Merson 法では残余局所誤差を記録のために保持することができる。これは各タームで 0.2 * ($ynp1[i] - $eta4[i]) となる。
rk4_auto では Ay + bt という形の線形の dydt 関数用にしか用いないのでこの性質を利用しない。
rk4_auto はダブルタイムステップごとに14回の関数評価が必要でタイムステップを調整するために13回の再計算が必要。
rk4_auto は ($t, $dt, @y) を返す。ただし $t, $dt, @y はタイムステップ終了時には新しい値となっている。
rk4_auto_midpoint()
rk4_auto は $dt でダブルタイムステップ計算をし、最終値を返す。通常は中間値を返さない。しかしなっがら滑らかなグラフを描くためや滑らかに動くディスプレイを示すためには中間点での値は有用である。それゆえ rk4_auto_midpointはそれらを取得する方法を提供する。
まず rk4_auto を呼び出さねばならないことに注意すること。$t+$dt における値を返し、続いて rk4_auto_midpoint を呼び出して $t+$dt/2 における値を求めるからである。言い換えれば時間順に呼び出す必要がある。例えば以下のようにする。
while ($t<$tfinal) {
($t, $dt, @y) = &rk4_auto(\@y, \&dydt, $t, $dt, $epsilon);
($t_midpoint, @y_midpoint) = &rk4_auto_midpoint();
&update_display($t_midpoint, @y_midpoint);
&update_display($t, @y);
}
rk4_auto_midpoint は ($t, @y) を返す。ただし $t と @y は rk4_auto を呼び出す前の値である。
dydt( $t, @y );
本サブルーチンは rk2, rk4, rk4_auto の2番めの引数として参照により渡される。もちろん名称は問題ではない。
以下の引数が想定されている。
$t：時刻（方程式が時間依存の場合）
@y：変数の値の配列
時間に従った変数の微分の配列を返す。
The following routines are not exported by default, but are exported under the I
import Math::RungeKutta qw(:ALL);
rk4_classical( \@y, \&dydt, $t, $dt )
引数と返り値は rk2 及び rk4 と同じである。
使用しているアルゴリズムは古典的でエレガントな4次の Runge-Kutta 法でタイムステップごとに4回の関数評価を行っている。
k0 = dt * F(y(n)) k1 = dt * F(y(n) + 0.5*k0) k2 = dt * F(y(n) + 0.5*k1) k3 = dt * F(y(n) + k2) y(n+1) = y(n) + (k0 + 2*k1 + 2*k2 + k3) / 6
rk4_ralston( \@y, \&dydt, $t, $dt )
引数と返り値は rk2 及び rk4 と同じである。
使用しているアルゴリズムは Ralston により開発されたもので、各タイムステップで誤差境界を最小化するために rk4_classical を最適化する。本モジュールはデフォルトの4次の rk4 を用いていない。それは Merson のアルゴリズムは精度が高い代わりに余分な関数評価があるたｍである。
build ディレクトリのサブディレクトリである examples/ に例示ようスクリプトがある。まずはそのスクリプトから始めるとよい。
sine-cosine
本スクリプトはアルゴリズムの選択のために Term::Clui を（矢印キー・Returnキー・終了のための q キーを使うため）用いている。積分のタイムステップ及び誤差基準として円の周りの単純な正弦波・余弦波を利用している。これは開発の間中使っていたスクリプトである。
three-body
本スクリプトは３体重力問題を表示するために vt100または xterm 'moveto' 及び 'reverse' シークウェンスを使っている。これは rk4_auto を用いている。互いに接近した２体の問題には短いタイムステップでよいからである。表示を滑らかにするために rk4_auto_midpoint も利用している。初期条件を変えることで出力がどれくらい依存しているかを経験することができる。
哀しいことに、数値シミュレーションにはマズイことがある。
最も基本的なものは不安定性である。もしタイムステップ $dt を導関数 @dydt が要求する時定数よりも大きく設定すると数値解は激しく振動し、元の方程式とは無関係な挙動を示す。このようなことが起こればさらに短い $dt を選択すべきである。
難問のいくつかにはいわゆる「不自然な」導関数を持つものがある。これは @dydt は広い時定数の幅を導出するときに起こる。不安定性を回避するために最小の時定数に対応する $dt を設定すると、長時間にわたる方程式の解を追跡することが困難になる。均衡点を見出し、長周期部分と短周期部分を分離すべきである。
同様に、数値積分は、例えばボールが壁に衝突するときのような突然の時定数の変化を望まない。各タイムステップで直接 @y 配列に手をいれ介入すべきである。例えば $y[17] が床にボールが衝突する高さとして、$y[20] が速度の垂直成分だとすると
if ($y[17] < 0.0) { $y[17] *= -0.9; $y[20] *= -0.9; }
このようにして数値積分が滑らかになるように調節する。
Peter J Billam, http://www.pjb.com.au/comp/contact.html
examples/sine-cosine 及び examples/three-body のスクリプト, http://www.pjb.com.au/, http://www.pjb.com.au/comp/, Math::WalshTransform, Math::Evol, Term::Clui, Crypt::Tea, http://www.xmds.org/ を見よ。
|
|
Updated : 2008/10/31 |